UNIVERSIDAD POLITÉCNICA DE MADRID ESCUELA TÉCNICA SUPERIOR DE INGENIEROS INDUSTRIALES Grado en Ingeniería en Tecnologías Industriales Especialidad Técnicas Energéticas TRABAJO DE FIN DE GRADO MODELIZACIÓN POR ELEMENTOS FINITOS DE PARQUES EÓLICOS Inés Mateos Canals Tutores: David Portillo García Ignacio Romero Olleros Convocatoria de Julio de 2016 “Vivimos en el fondo de un océano de aire” Evangelista Torricelli, 1608-1647 Agradecimientos A mi tutor, David, por guiarme a través del proyecto y tener paciencia ante mis frustraciones con la programación. A mi familia: a los ingenieros, por compartir mi entusiasmo, y a los no ingenieros, por aguantarnos. A mis amigos, por sacarme del pozo cuando es necesario y estar dispuestos a chapotear conmigo cuando no. Inés Mateos Canals 5 Resumen Introducción La relevancia de la energía eólica ha aumentado drásticamente en los últimos años debido a los desafíos medioambientales y la demanda energética creciente. Un aspecto fundamental de esta fuente de energía es su carácter sostenible, al ser el viento consecuencia del calentamiento de la atmósfera por la radiación solar y de la rotación de la Tierra. Esta característica supone un gran atractivo con vistas al futuro. La apuesta por la instalación de parques eólicos supone una inversión en innovación tecnológica en este campo, en el que cobran especial importancia las herramientas de simulación. El diseño y la optimización de aerogeneradores y parques eólicos ha encontrado un fuerte aliado en la Mecánica de fluidos computacional, una rama de la mecánica de fluidos que emplea análisis numérico y algoritmos para resolver y estudiar problemas relacionados con el flujo de fluidos. La modelización de los aerogeneradores mediante métodos analíticos permite estudiar el comportamiento del viento en torno a un aerogenerador e incluso, en un parque eólico que integre un elevado número de turbinas. Objetivos El fin último de este trabajo es la modelización y simulación de un parque eólico para posibilitar el desarrollo de una herramienta de optimización sobre el programa de elementos finitos IRIS. Este objetivo se aborda a través de la definición e implementación de un modelo de aerogenerador y la realización de simulaciones de complejidad creciente, empezando por un aerogenerador aislado y concluyendo con la simulación de un parque eólico existente. Metodología En el desarrollo de este Trabajo de Fin de Grado (TFG) pueden diferenciarse tres bloques. El primero tiene un enfoque teórico y comprende los dos primeros capítulos del presente documento, se trata de estudio analítico de la aerodinámica de un aerogenerador y de una serie de modelos matemáticos, cuyo objeto es la modelización de un aerogenerador para su simulación por ordenador. En esta parte se pretende comprender el mecanismo de extracción de potencia que tiene lugar en el rotor de un aerogenerador, de forma que pueda crearse un modelo del mismo. El segundo y tercer bloque cuentan con un carácter práctico al tratarse de la aplicación de uno de los modelos estudiados, con ciertas simplificaciones, para la obtención de los resultados del trabajo. En el segundo bloque se han llevado a cabo varias simulaciones de un aerogenerador aislado y se han analizado las predicciones de potencia generada que resultan de ellas mediante el contraste con las especificaciones técnicas del aerogenerador y otros resultados obtenidos en estudios similaInés Mateos Canals i RESUMEN res. Por otro lado, se ha estudiado la influencia de determinados parámetros del modelo sobre las estimaciones realizadas, así como la relevancia de la elección del mallado empleado. El tercer bloque trata la simulación de varios aerogeneradores con vistas a observar las consecuencias que tiene agruparlos, en términos de interacción entre las estelas y potencia generada. Finalmente, la simulación de un parque eólico de mediano tamaño supone el nexo entre el trabajo realizado y la continuación lógica del mismo, la optimización de un parque eólico. Es en este punto en el que se han alcanzado los límites de la capacidad de computación de la que dispone el departamento. Se exponen los resultados obtenidos y se detallan las inestabilidades que derivan de las limitaciones mencionadas. Resultados y conclusiones El estudio de los modelos básicos de la aerodinámica de un aerogenerador ha permitido identificar el modelo del disco actuador (ADM) como opción adecuada para cumplir los objetivos de este trabajo. El modelo satisface el requisito de lograr un equilibrio entre simplicidad, en términos de exigencia de poder computacional, y precisión, en tanto que los resultados obtenidos permitan llevar a cabo análisis válidos de escenarios de distribución de aerogeneradores. La implementación en IRIS del modelo ADM axisimétrico ha supuesto la realización de una serie de simplificaciones. Cabe destacar que se trata de un modelo limitado al rotor (buje y palas) del aerogenerador que no considera la influencia de la topografía del terreno, el efecto del suelo y las características de la capa límite atmosférica. Por otro lado, aprovechando la simetría axial del flujo incidente sobre el aerogenerador, se ha resuelto el problema bidimensional en coordenadas cartesianas. La simplificación anterior permite la simulación de parques eólicos en dos dimensiones. En la figuras que acompañan a este texto se muestran, a modo de ejemplo, los campos de velocidad y presión que resultan de la simulación de un dominio de 900x600 m, en el que el flujo de aire ha sido perturbado por la presencia de un aerogenerador. (a) Campo de velocidades (b) Campo de presiones La aptitud del modelo ha sido estudiada por medio de simulaciones, que han permitido valorar la robustez del modelo ante variaciones del mallado y del parámetro (relacionado con la distribución de las fuerzas elementales) y la validez de la representación de la interacción de aerogeneradores ii Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos cercanos entre sí. Como parte de este estudio, se han llevado a cabo el análisis y contraste de resultados con las especificaciones técnicas y las valoraciones de otros trabajos. Entre las observaciones sobre el modelo, cabe destacar las consideraciones tomadas a la hora de valorar el salto de presiones que tiene lugar en el rotor del aerogenerador, dado que se trata de una variable que afecta a las predicciones de potencia generada. Se ha podido comprobar que, para evitar inestabilidades numéricas que conducen a estimaciones erróneas de la potencia generada, es preciso emplear los valores de presión de elementos del mallado que se encuentran a una distancia superior a la equivalente al 10 % de la diámetro del disco, ∆x. Por otro lado, la utilización de un valor = 2∆x conduce a la obtención de buenas estimaciones de potencia, por lo que éste ha sido el valor empleado en las simulaciones del presente trabajo. El resultado culminante de este proyecto es la correcta simulación de un conjunto de aerogeneradores existente, el parque eólico de Lillgrund. Dada la gran extensión del dominio simulado (3800x4160 m), los recursos de memoria disponibles han limitado el mallado a ∆x = 8m, definición que proporciona buenos resultados cualitativos pero no resulta suficientemente precisa a la hora de estimar la potencia. A continuación se muestran las representaciones gráficas de los resultados obtenidos. (c) Campo de velocidades (d) Campo de presiones Sobre el trabajo realizado en el desarrollo de este proyecto es posible llevar a cabo modificaciones que permitan prescindir de ciertas simplificaciones tomadas en la definición del modelo y, cabe esperar, aumenten la precisión de las resultados. Un claro ejemplo de modelización de mayor complejidad con un gran interés es la simulación en tres dimensiones del campo fluido. Se ha alcanzado el objetivo general de sentar los cimientos sobre los que será posible construir una herramienta de optimización de parques eólicos basada en el software IRIS. Inés Mateos Canals iii RESUMEN Palabras clave: Energía eólica, aerogeneradores, parque eólico, modelización de aerogeneradores, simulación, ecuaciones de Navier-Stokes, elementos finitos, IRIS. Códigos UNESCO: 220404 Mecánica de Fluidos 330101 Aerodinámica 330102 Cargas aerodinámicas 332202 Generación de energía 332203 Generadores de energía iv Escuela Técnica Superior de Ingenieros Industriales (UPM) Índice general Resumen I Introducción Energía eólica . . . . . . . . . . . . . . . . . . . . Simulación de aerogeneradores y parques eólicos Objetivos generales . . . . . . . . . . . . . . . . . Objetivos detallados . . . . . . . . . . . . . . . . VII . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii . viii . x . x 1. Aerodinámica de aerogeneradores 1.1. Teoría del disco actuador . . . . . . . . . . . . . . 1.2. Teoría del momento cinético . . . . . . . . . . . . . 1.3. Perfiles aerodinámicos . . . . . . . . . . . . . . . . 1.3.1. Sustentación . . . . . . . . . . . . . . . . . 1.3.2. Resistencia . . . . . . . . . . . . . . . . . . 1.3.3. Coeficientes adimensionales . . . . . . . . . 1.3.4. Torbellinos . . . . . . . . . . . . . . . . . . 1.4. Teoría del elemento de pala . . . . . . . . . . . . . 1.5. Teoría del momento sobre elemento de pala (BEM) 1.5.1. Pérdidas en la punta de la pala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 5 6 7 8 8 9 10 11 12 2. Modelización de aerogeneradores 2.1. Simulación de la capa límite atmosférica . . 2.2. Modelo del disco actuador (ADM) . . . . . 2.3. Modelo de la linea actuadora (ALM) . . . . 2.4. Modelo de la superficie actuadora (ASM) . 2.5. Comparación del ADM y el ALM . . . . . . 2.6. Implementación de modelos . . . . . . . . . 2.6.1. Parámetro de la función gaussiana 2.6.2. Implementación del ALM . . . . . . 2.6.3. Implementación del ADM . . . . . . 2.6.4. Método numérico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 14 14 15 16 16 17 17 18 19 19 3. Simulación con IRIS 3.1. Método de elementos finitos estabilizado 3.2. Modelo ADM . . . . . . . . . . . . . . . 3.3. Implementación del modelo en IRIS . . . 3.3.1. Simplificaciones . . . . . . . . . . 3.3.2. Elementos actuadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 21 22 23 23 24 Inés Mateos Canals . . . . . . . . . . v ÍNDICE GENERAL 3.3.3. Ampliación de IRIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Simulaciones realizadas 4.1. Aerogenerador NREL de 5 MW . . . . . . . . . 4.1.1. Propiedades estructurales de las palas . 4.1.2. Propiedades aerodinámicas de las palas 4.1.3. Curva de potencia . . . . . . . . . . . . 4.1.4. Parámetros: α, φ, CL y CD . . . . . . . 4.2. Simulación de un aerogenerador aislado . . . . 4.2.1. Condiciones de contorno y referencia . . 4.2.2. Velocidad . . . . . . . . . . . . . . . . . 4.2.3. Presión . . . . . . . . . . . . . . . . . . 4.2.4. Potencia . . . . . . . . . . . . . . . . . . 4.3. Simulaciones de aerogeneradores agrupados . . 4.3.1. Dos aerogeneradores en serie . . . . . . . . . . . . . . . . . . 5. Simulación de un parque eólico 5.1. Aerogenerador Siemens SWT-2.3-93 . . . . . . . 5.1.1. Propiedades aerodinámicas de las palas . 5.1.2. Curva de potencia . . . . . . . . . . . . . 5.1.3. Simulación del aerogenerador SWT-2.3-93 5.2. Resultados de la simulación de Lillgrund . . . . . 5.2.1. Velocidad . . . . . . . . . . . . . . . . . . 5.2.2. Presión . . . . . . . . . . . . . . . . . . . 5.2.3. Potencia . . . . . . . . . . . . . . . . . . . Conclusiones y líneas futuras 26 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 29 30 31 32 32 33 33 33 35 36 41 41 . . . . . . . . . . . . . . . aislado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 47 48 49 50 52 53 56 57 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Planificación temporal y presupuesto 65 Diagrama de Gantt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Presupuesto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Bibliografía 69 Índice de figuras 71 Índice de cuadros 73 Abreviaturas, unidades y acrónimos 75 Índice de símbolos 77 Glosario 79 vi Escuela Técnica Superior de Ingenieros Industriales (UPM) Introducción Energía eólica La Unión Europea apuesta de forma decidida por las energías renovables y en particular por la energía eólica. En el año de 2015 se alcanzaron los 141.6 GW de potencia instalada, 11 de los cuales corresponden a parques eólicos marinos. La generación de energía eólica ha superado a la hidráulica para situarse en el tercer puesto de las fuentes de energía más importantes en la UE, aportando un 15.6 % de la demanda energética total de la unión. Por otro lado, la instalación de parques eólicos supone un tercio de la totalidad de las instalaciones de generación de energía llevadas a cabo desde el año 2000 [1]. Figura 1: Evolución de la potencia eólica instalada en la Unión Europea. Fuente: EWEA [1] La energía eólica consiste en la obtención de electricidad por transformación de la energía cinética de las corrientes de aire. Al tratarse el viento de una fuente de energía inagotable, la energía eólica pertenece al grupo de las energías renovables. La potencia del viento ha sido empleada a lo largo de la historia en las velas de los barcos, en molinos y bombas de agua. Sin embargo, la transformación de la potencia mecánica en potencia eléctrica mediante un aerogenerador es característica de la energía eólica. El primer molino que se empleó como aerogenerador data de 1887 y fue construido por el catedrático James Blyth, del Anderson’s College, en Glasgow, Escocia [2]. El aerogenerador tenía una altura de 10 m y contaba con velas de tela. Por primera vez, se cargó un conjunto de baterías y se proporcionó electricidad de origen eólico a una vivienda. Con la evolución de la electricidad, la energía eólica contó con nuevas aplicaciones, como el suministro de electricidad en edificios aislados de la red eléctrica. A lo largo del siglo XX y de forma paralela, estaciones de generación de potencia eólica fueron construidas para generar energía Inés Mateos Canals vii INTRODUCCIÓN eléctrica para granjas o viviendas, así como aerogeneradores de mayor tamaño que podían ser conectados a la red. A día de hoy, los aerogeneradores operan en todos los rangos, desde instalaciones pequeñas para la recarga de baterías en lugares remotos hasta los grandes parques eólicos marinos que suministran electricidad a la red eléctrica nacional. Los aerogeneradores se pueden clasificar según las palas y la orientación del eje. El modelo más extendido es el aerogenerador de eje horizontal, conocido por su sigla en inglés HAWT (horizontalaxis wind turbine). El rotor de la turbina y el generador eléctrico se sitúan en lo alto de una torre de altura entre 60 y 90 m y la mayoría cuentan con un multiplicador que transforma la rotación lenta (entre 10 y 22 revoluciones por minuto) en una rotación más rápida y adecuada para el generador. Por lo general, las turbinas empleadas comercialmente para la producción de electricidad cuentan con tres palas y se orientan según la dirección del viento mediante motores controlados remotamente. Los parques eólicos son agrupaciones de aerogeneradores para la producción de energía eléctrica. Las turbinas están conectadas entre sí a media tensión (en torno a 34.5 kV), a un sistema de almacenamiento de energía y una red de comunicación. Una subestación se encarga de elevar la tensión mediante transformadores para su transporte. Un parque eólico de gran tamaño puede consistir en varios cientos de aerogeneradores distribuidos sobre grandes superficies, que pueden ser empleadas para agricultura. Este es el caso del Parque eólico de Gansu, en China, cuya potencia nominal alcanza los 6.000 MW [3]. Por otro lado, los parques eólicos pueden localizarse en alta mar, con el objetivo de aprovechar las corrientes de viento de gran intensidad y evitar el impacto ambiental que conlleva la construcción sobre tierra firme. Las dificultades constructivas y los elevados costes de mantenimiento suponen inconvenientes de estas instalaciones. El parque eólico marino (en inglés, offshore) de mayor tamaño es el London Array y suministra 630 MW [4]. Antecedentes: simulación de aerogeneradores y parques eólicos La investigación cuenta con un papel fundamental a la hora de diseñar y optimizar las instalaciones de generación de energía eólica. Empresas y universidades centran sus esfuerzos en llevar a cabo estudios que permitan predecir y analizar el comportamiento de las corrientes de aire y la producción de potencia posible por medio de aerogeneradores. Esto se consigue mediante el desarrollo de modelos de la interacción pala-fluido y la resolución de las ecuaciones de flujo de Navier-Stokes aplicadas de rotores de aerogeneradores. El análisis de mecánica de fluidos computacional cobra especial importancia como herramienta de simulación del flujo dinámico del aire en la capa límite atmosférica. Los métodos que tienen como objeto predecir el comportamiento del aire que interacciona con un rotor tienen su origen en el concepto de disco actuador, debido a Freude [5], que representa el rotor como una distribución de fuerzas equivalentes en la superficie de un disco permeable al aire de grosor nulo. Análisis posteriores, a principios del siglo XX, de Lanchester [6] y Betz [7], dieron a conocer la máxima extracción de energía posible mediante un aerogenerador, en torno al 59.3 % de la energía cinética del flujo incidente. La modelización del flujo de aire a través de rotores dio un gran paso hacia delante con las teorías del momento cinético generalizada (en inglés, General Momentum Theory) y del momento sobre el elemento de pala (BEM, del inglés Blade Element Momentum), de viii Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Glauert [8], en 1930. La aplicación de las ecuaciones de balance a tubos anulares de corriente proporciona una mejor aproximación a la distribución de esfuerzos a lo largo de las palas del aerogenerador. En años posteriores al desarrollo de estas teorías ha continuado el estudio de modelos, que han ido ganando complejidad según se acercaban a una mejor aproximación a la realidad. La capacidad y las prestaciones de los sistemas de computación mejoran año a año y, con ello, la calidad de las simulaciones. Aunque el modelo BEM es el más empleado en la industria, se han elaborado una gran variedad de modelos avanzados de rotores de aerogeneradores. Estos modelos pueden ser clasificados según las consideraciones relacionadas con la viscosidad del aire. El uso de un modelo u otro tiene consecuencias muy relevantes sobre el tiempo de simulación e incluso la dificultad de aplicación al campo fluido, debido a las singularidades y los fallos de convergencia. El método del disco actuador generalizado supone una extensión de la teoría BEM con consideraciones de no viscosidad. La principal diferencia es que la independencia entre los anillos de corriente se ve sustituida por la resolución de las ecuaciones de Euler o Navier-Stokes. Versiones con simetría axial de este método han sido desarrolladas y resueltas de forma analítica y semi-analítica o por métodos de elementos finitos (FEM). Los métodos de elementos finitos facilitan el desarrollo natural e inestable de las estelas, mientras que los métodos analíticos suelen ser resueltos para condiciones estables. Con el objetivo de evitar los problemas derivados de emplear parámetros corregidos o calibrados, algunos modelos que consideran la viscosidad han sido desarrollados para evaluar el flujo de aire en torno a un aerogenerador; estos modelos se deben a Sørensen [9], Sankar [10] y a Hansen [11] y sus colaboradores, entre otros. También se han llevado a cabo simulaciones en tres dimensiones mediante las ecuaciones de Navier-Stokes promediadas de Reynolds (RANS). A pesar de que los métodos basados en las ecuaciones RANS son capaces de reproducir las condiciones previas a la incidencia sobre el rotor, aún fallan en el tratamiento de las condiciones tras el plano del rotor, debido a una modelización de la turbulencia y una resolución de mallado inadecuadas. La fiabilidad de los modelos de rotor está condicionada por las hipótesis de partida y los datos de los perfiles aerodinámicos sobre los que trabajan. Discrepancias considerables son resultado de las limitaciones en los datos, de una modelización de simplicidad excesiva o de hipótesis no válidas del modelo. Algunas de las hipótesis a las que se hace referencia pueden ser consideraciones de simetría axial, independencia de los anillos de corriente entre sí o ajustes de tipo experimental, relacionados con las pérdidas en punta de pala u otras correcciones. Las predicciones del modelo de Sørensen [12], que sólo emplea datos de la geometría de la pala y de la corriente fluida incidente, resultan acertadas para pares moderados. La simulación de la turbulencia supone un tema clave en estos modelos, puesto que los requisitos computacionales las simulaciones numéricas directas (DNS, del inglés Direct Numerical Simulation) y de grandes torbellinos (LES, del inglés Large-Eddy Simulations) están por encima de la capacidad de los ordenadores de los que se dispone. Las limitaciones impuestas por los ordenadores han llevado al desarrollo de modelos híbridos, como el de Johansen [13], que aúnan características de de modelos convencionales de turbulencia y de LES. Estos métodos, conocidos como DES (del inglés Detached-Eddy Simulation), incluyen en mayor medida el comportamiento dinámico y tridimensional del flujo. Los métodos basados en las ecuaciones de Navier-Stokes en tres dimensiones tienen un futuro prometedor, aunque ligado al avance de las herramientas de simulación. Inés Mateos Canals ix INTRODUCCIÓN En la actualidad, a la espera de ordenadores con mayor poder de computación, gran parte de los esfuerzos se destina a la mejora de métodos más simples, campo en el que vuelve a destacar el trabajo debido a Sørensen y su equipo [14]. Los cambios realizados van desde una mayor aproximación a ciertas geometrías de rotor hasta correcciones de punta e influencia de la presión. Objetivos generales Este trabajo se enmarca dentro de un proyecto académico más amplio que busca construir una herramienta de simulación para la optimización de un parque eólico, sobre el programa del elementos finitos IRIS y los medios de computación disponibles en el Departamento de Mecánica Estructural y Construcciones Industriales de la Escuela Técnica Superior de Ingenieros Industriales de Madrid. Este objetivo general se aborda en tres etapas: Modelización del aerogenerador aislado. Simulación de aerogeneradores aislados y agrupados. Análisis de escenarios de distribución con vistas a la optimización de un parque eólico. En la modelización se pretende alcanzar un equilibrio entre simplicidad y precisión. Es importante señalar que el modelo debe facilitar la comparación de disposiciones geométricas de agrupaciones de generadores con vistas a que la interacción entre ellos se minimice y la potencia generada por el conjunto alcance un valor máximo. La validez de estas comparaciones es más importante que la precisión numérica de las estimaciones de potencia. Por otro lado, las restricciones de capacidad de computación (tiempo de ejecución y demanda de memoria) también imponen la simplificación. Objetivos detallados En el desarrollo del presente trabajo conviene establecer la siguiente secuencia de objetivos detallados: Comprender el mecanismo de extracción de energía del viento por medio de un aerogenerador por aplicación de las ecuaciones de conservación y los modelos teóricos propuestos hasta hoy. Realizar un estudio sobre los métodos de modelización de aerogeneradores disponibles para la simulación mediante herramientas de Mecánica de Fluidos Computacional (CFD). Simular la evolución de las propiedades de una corriente fluida que incide sobre un rotor haciendo uso de un método de modelización adecuado a los medios disponibles y empleando el programa de elementos finitos IRIS. Analizar y comparar los resultados experimentales de las simulaciones realizadas con las especificaciones técnicas correspondientes y con aquellos resultados obtenidos en otros estudios a través de modelizaciones similares. Estudiar la influencia de los parámetros considerados sobre las estimaciones de potencia que resultan de las simulaciones. Estimar la potencia generada por un parque eólico existente y contrastar resultados. x Escuela Técnica Superior de Ingenieros Industriales (UPM) Capítulo 1 Aerodinámica de aerogeneradores A continuación se expone el mecanismo básico de obtención de energía en un aerogenerador de eje horizontal según el enfoque analítico clásico. La potencia eléctrica generada en un aerogenerador es el resultado de la extracción de energía cinética del viento. Para conocer la magnitud de la energía extraíble, es necesario conocer la energía cinética total que posee una corriente de aire. El análisis que se desarrolla a continuación se basa en la Teoría del Momento Cinético, que requiere hacer una serie de hipótesis de situación ideal del flujo: Aire como fluido ideal, sin viscosidad. Esta hipótesis no se aleja mucho de la realidad, puesto que el movimiento del aire alrededor del aerogenerador tiene lugar a elevados números de Reynolds, es decir, las fuerzas viscosas resultan despreciables frente a las fuerzas inerciales. Incidencia del viento unidimensional, a presión, densidad y velocidad uniformes. Se desprecia la influencia de obstáculos, incluido el suelo, y la turbulencia que a la que éstos darían lugar. Aire como fluido incompresible. Los efectos disipativos debidos a la compresibilidad del flujo no empiezan a ser apreciables hasta que el número de Mach inicidente a la pala es muy próximo a la unidad (velocidades cercanas a la del sonido). Flujo estacionario. Indpendencia temporal. Estela sin rotación. El intercambio de par en la turbina resulta en un movimiento de rotación alrededor del eje por parte del fluido tras atravesarla, en la situación idealizada, no se considera este efecto, que en realidad afecta a la eficiencia. 1.1. Teoría del disco actuador La teoría del disco actuador permite calcular el empuje ejercido por el viento sobre un rotor ideal de aerogenerador, así como la producción de potencia del mismo. El aire que pasa a través del aerogenerador cederá, en la situación ideal, la misma cantidad de energía en cualquier punto de la corriente. Según esta teoría, se considera el rotor como el disco circunscrito a los extremos de las palas, a través del cual circula la corriente, cediendo energía a su paso. En adelante, el subíndice 0 hará referencia a las condiciones iniciales de la corriente de aíre, el 1, a las condiciones en las cercanías del disco actuador y el 2 a aquellas aguas abajo. Por aplicación Inés Mateos Canals 1 CAPÍTULO 1 1.1. TEORÍA DEL DISCO ACTUADOR de las leyes de conservación al flujo que atraviesa el disco de sección A = πΦ2 /4: Continuidad. El caudal másico debe mantenerse constante a lo largo del tubo de corriente. Sólo se considera la componente axial de la velocidad: ṁ = ρA0 v0 = ρA1 v1 = ρA2 v2 A0 v0 = Av1 = A2 v2 (1.1) Dado que, justo antes y justo después de atravesar el disco, la sección transversal del flujo es la misma, la velocidad delante y dentrás del rotor debe ser la misma. En adelante, en este texto se empleará la notación simplificada de v, para hacer referencia a v1 . Cantidad de movimiento. La fuerza ejercida por el disco sobre el fluido (D) da lugar a una variación de cantidad de movimiento del flujo. D = −ṁ(v − v0 ) = ρAv(v − v0 ) (1.2) Energía. El flujo está caracterizado por se estacionario, ideal e incompresible, por lo que resulta aplicable la ecuación de Bernouilli a una linea de corriente. Según esta ecuación, la presión de remanso, de parada o total, que se define como la suma de las presiones estática y dinámica, p + 21 ρv 2 , se conserva en una linea de corriente. A través del disco, la ecuación de conservación anterior no es aplicable, por que no se conserva la energía por unidad de masa, sin embargo, esta sí es aplicable entre las secciones transversales inicial y anterior al disco, así como entre la posterior al mismo y la final: 1 p1 + ρv12 = p+ + 2 1 − p + ρv12 = p2 + 2 1 2 v 2 1 2 v 2 2 (1.3) Restando las ecuaciones anteriores miembro a miembro, se obtiene la expresión que permite conocer el salto de presiones que tiene lugar a través del disco: 1 p+ − p− = ρ(v02 − v22 ) 2 (1.4) Equilibrio del disco. El sumatorio de las fuerzas sobre el disco debe anularse dada su situación estacionaria de velocidad nula. D = (p+ − p− )A (1.5) Igualando las expresiones de la resistencia aerodinámica (D) que dan las ecuaciones 1.2 y 1.5, se obtiene la expresión: ρAv(v − v0 ) = (p+ − p− )A (1.6) e, introduciendo la expresión del salto de presiones de 1.4: 1 ρAv(v − v0 ) = ρ(v02 − v22 )A 2 2 (1.7) Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Figura 1.1: Evolución de la velocidad y la presión del flujo [15]. La simplificación de esta igualdad lleva a la expresión: v= v02 + v22 2 (1.8) El disco constituye un obstáculo para la corriente incidente, causando una divergencia de las lineas de corriente, que sólo puede ser mantenida con gradientes de presión en los que las lineas se curvan hacia zonas de menor presión. La presión delante del disco, p+ , es superior a la de la corriente libre, p0 . El aumento de presión trae consigo una disminución de velocidad en el plano de disco, que resulta menor que la de la corriente libre, v < v0 . La extracción de energía cinética en el disco causa una caída de presión detrás de éste, por lo que p− < p+ . Las ecuaciones de conservación anteriores establecen que la presión p− es menor que la de la corriente incidente, p0 , por lo que la sección transversal del tubo de flujo aumenta progresivamente tras el disco, de forma que la velocidad disminuya y la presión se eleve. La igualdad de presiones se consigue aguas abajo, donde p2 = p0 , sin embargo, esta igualdad de presiones no implica igualdad de velocidades. La velocidad aguas abajo, v2 , es menor que la da incidente, v0 , puesto que el fluido ha perdido energía al atravesar el disco y la presión de parada es menor que la de la corriente incidente. Se cumple: v0 > v > v 2 (1.9) En efecto, la ecuación 1.8 indica que la velocidad inducida aguas abajo, v2 , es el doble que en el plano del disco, por lo que el tubo de corriente debe duplicar su sección transversal. En este punto, resulta útil introducir el parámetro adimensional a, definido como la fracción de reducción de la velocidad local en el entorno del disco, es decir, el defecto de velocidad, se conoce como coeficiente de inducción axial. v = v0 (1 − a) (1.10) v2 = v0 (1 − 2a) (1.11) De acuerdo con esta definición: Inés Mateos Canals 3 CAPÍTULO 1 1.1. TEORÍA DEL DISCO ACTUADOR Figura 1.2: Relación de velocidades y areas en la corriente [16]. Considerando la variación de energía cinética entre aguas arriba y aguas abajo como el trabajo elemental producido por la diferencia de presiones en la superficie del disco, la potencia extraída es: P = (p+ − p− )Av = Dv (1.12) Sustituyendo el valor de la diferencia de presiones por el visto en la ecuación 1.4: 1 1 P = ρ(v02 − v22 )Av = ṁ (v02 − v22 ) 2 2 (1.13) Esta igualdad indica que la potencia extraída es la diferencia de flujo de energía cinética entre la corriente incidente y la saliente, que es coherente en cuanto a la conservación de la energía. En función de la inducción axial: P = ρv0 A20 (1 − a)(v0 − v2 ) (1.14) y, de acuerdo con 1.11: 1 P = ρv03 A[4a(1 − a)2 ] (1.15) 2 El cociente entre la potencia extraída y la disponible a través del disco se denomina coeficiente de potencia Cp , en función de a: P 1 3 2 ρv0 A = 4a(1 − a)2 = Cp (1.16) La diferenciación de la ecuación anterior para obtener el máximo resulta: ∂Cp 1 16 = 4(1 − a)2 − 8a(1 − a) = 0 → a = ; Cpmax = w 0, 5926 (1.17) ∂a 3 27 El límite superior de Cp es conocido como límite de Betz e indica que aproximadamente un 60 % de la energía del viento es convertible en potencia por medio de una turbina como máximo. Por otro lado, el coeficiente de resistencia (en inglés, drag coefficient), se define como la razón de la fuerza resistente y la fuerza total de la corriente: CD = 4 D 1 2 2 ρv A = 4a(1 − a) (1.18) Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Figura 1.3: Velocidad angular impuesta sobre una sección anular [17] Que, en el caso de potencia máxima, toma el valor: CD = D 1 2 2 ρv A = 8 9 (1.19) Es necesario aplicar un valor próximo a la presión dinámica de la corriente sobre la superficie del rotor para conocer la carga de empuje aerodinámico sobre la turbina en el punto de máxima extracción de potencia. Este valor es inferior al coeficiente de resistencia aerodinámica de una placa plana circular impermeable perpendicular a la corriente. 1.2. Teoría del momento cinético La teoría de la cantidad de movimiento no permite determinar el par Q que ejerce el fluido sobre las palas y el par reacción que ésas ejercen sobre el fluido, causando un cambio en la cantidad de movimiento con respecto al eje de la corriente, que coincide con el del rotor. La componente de la velocidad capaz de dar momento es la tangencial, vθ en un sistema de coordenadas cilíndrico (r, θ, z). Se supone una situación de uniformidad tangencial, es decir, velocidad tangencial independiente del ángulo θ, aunque dependiente del radio, r, de la forma vp = Ωr. Queda definido entonces un coeficiente de inducción tangencial, a0 , que mide la velocidad angular, ω, impuesta a la corriente como una fracción de la velocidad angular de giro del rotor, Ω: a0 = ω vθ /r = Ω vp /r (1.20) Para conocer la potencia extraída, se aplica la conservación del momento cinético a un volumen de control constituido por una sección transversal dA antes del rotor con forma de corona circular de radio r y espesor dr, dA = 2πrdr, el tubo anular de corriente que pasa por su perímetro exterior e interior y la sección transversal dA resultante inmediatamente aguas abajo. Justo delante del rotor, la velocidad de la corriente carece de componente tangencial y justo detrás, la componente tangencial tiene un valor medio vθ . Según la teoría del momento cinético, el diferencial de par, dQ Inés Mateos Canals 5 CAPÍTULO 1 1.3. PERFILES AERODINÁMICOS aplicado por los elementos de la pala entre r y r + dr sobre el fluido corresponde con la variación de momento angular del flujo. Teniendo en cuenta 1.10 y 1.20: dQ = vθ rdṁ dṁ = ρvdA ⇒ dQ = 2πr3 ρv0 (1 − a)a0 Ωdr (1.21) Esta relación permite comprobar que la contribución de las partes de la pala más cercanas al extremo contribuyen más al par que las próximas al buje, dada la dependencia con r3 . La potencia es el producto del par por la velocidad angular, por lo que el diferencial de potencia producida por el elemento de pala anterior es: dP = ΩdQ (1.22) La integración a lo largo de la pala y la suma de las contribuciones de cada pala resultarían en el par y la potencia totales extraídos. Al no conocer los valores de a y a0 es necesario recurrir a modelos más detallados. Estos modelos consideran tanto a como a0 funciones del radio, lo que resulta posible dados los elevados números de Reynolds del fujo, que implican la existencia de gradientes de velocidad sin gran oposición por parte de la viscosidad. La dependencia de a con el radio no se considera en la teoría de la cantidad de movimiento, en la que la velocidad axial se consideraba constante en toda la superficie del disco. La rotación de la estela supone una reducción en la energía extraíble, puesto que la energía cinética de giro no es aprovechable por los aerogeneradores de eje horizontal. La rotación de la estela y el frenado de la corriente se manifiestan por un sistema de torbellinos ligados a las palas, que se desprenden de ellas por el buje y por los extremos de las mismas. Tras desprenderse, los torbellinos quedan libres y se ven arastrados por la corriente, que gira tras el plano de movimiento de las palas, por lo que describen trayectorias helicoidales. 1.3. Perfiles aerodinámicos Las teorías descritas previamente no tienen en cuenta la geometría de las palas o su número. Resulta interesante aproximar el flujo alrededor de palas esbeltas (en las que el radio del rotor es muy superior al ancho de la pala) a aquel que existiría si el ancho de la pala fuera infinitesimal con respecto al radio, resultando un flujo bidemensional en una sección a distancia r = cte. del eje del rotor. El perfil de la pala es el corte de la misma con una superficie r = cte., debido a que las palas de las turbinas rápidas no son muy anchas y son esbeltas, esta sección es prácticamente plana. En la zona cercana al buje esto no se cumple, pero en esa zona la actividad aerodinámica es muy baja. Los perfiles empleados en la obtención de energía eólica siguen la tecnología aeronáutica de perfiles de alas y hélices de baja velocidad, aunque se han desarrollado perfiles específicos para turbinas de viento. A velocidades subsónicas, la forma adecuada de las palas sigue las siguientes lineas generales: Borde de ataque o de entrada: borde enfrentado a la corriente, redondeado y de forma lisa y suave. Permite al perfil actuar con elevado rendimiento a distintos ángulos de orientación a la corriente. 6 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Borde de fuga o de salida: borde en el extremo de sotaviento, afilado. Su forma tiene por objeto evitar que la corriente rodee el perfil, salvo con un desprendimiento intenso. Dirige la corriente y permite reducir la resistencia al avance. Cuerda del perfil: línea recta de longitud c que une el borde de ataque con el borde de fuga. Línea de curvatura: línea equidistante entre ambos lados del perfil, el extradós, el más convexo, y el intradós, el menos, pudiendo llegar a ser cóncavo. La distancia máxima a la cuerda define la curvatura máxima del perfil cmax , que suele estar entre el 25 % y el 50 % de la cuerda, va de cero (perfil simétrico) a un 15 % de la cuerda. Distribución de espesor: distancia entre extradós e intradós. Normalmente es una curva suave que alcanza su máximo tmax entre el 20 % al 40 % de la cuerda. Cuanto más grueso es un perfil, mayor resistencia aerodinámica tiene, pero permite una estructura más rígida para soportar cargas. Los perfiles con curvatura tienen un óptimo de espesor por consideraciones aerodinámicas y de resistencia mecánica. Figura 1.4: Geometría de un perfil aerodinámico subsónico [15] Las fuerzas aerodinámicas que pueden aparecer sobre un perfil son: Sustentación, L (de lift en inglés): componente perpendicular a la corriente incidente aguas arriba, considerada positiva si es hacia extradós. Resistencia, D (de drag en inglés): componente en la misma dirección y sentido que la corriente incidente sin perturbar. 1.3.1. Sustentación La componente de sustentación aparece si el perfil forma un ángulo de ataque α con la corriente sin perturbar, ángulo que se mide respecto a la línea de sustentación nula, aquella dirección de la corriente relativa al perfil que no provoca sustentación y forma un ángulo α0 con la cuerda. Un perfil con curvatura presenta una línea de sustentación nula que pasa por su borde de fuga y por un punto de la línea de curvatura próximo a su máximo. Al inicir la corriente sobre la pala, aparece también un momento de encabritamiento, M , que tiende a aumentar el ángulo de ataque. Si el perfil aerodinámico se enfrenta a corrientes con ángulos de ataque pequeños, la corriente le rodea suavemente describiendo un flujo laminar aproximadamente bidimensional, salvo en una Inés Mateos Canals 7 CAPÍTULO 1 1.3. PERFILES AERODINÁMICOS capa muy delgada, denominada capa límite, donde los efectos viscosos son dominantes. Al margen de la capa límite, predominan las fuerzas inerciales y es aplicable la ecuación de Berouilli. La sustentación es consecuencia de la mayor aceleración experimentada por el fluido en contacto con el extradós, frente al intradós, al suponer mayor estechamiento a su paso. El flujo genera una distribución de presión sobre el perfil, cuya integración es la sustentación neta. Del orden de dos tercios de la sustentación se generan por la depresión en el extradós (succión) y un tercio por la sobrepresión en el intradós. La deflexión de la corriente es local, se circunscribe a las inmediaciones del perfil. El flujo tiende al de la corriente sin perturbar al aumentar su distancia al perfil. Si el ángulo de ataque α es excesivo, superior a 15◦ o 20◦ , se produce el desprendimiento de la corriente del extradós, dejando de ejercer succión. El movimiento adquiere carácter turbulento tridimensional, altera la distribución de presiones y el funcionamiento del perfil. Se pierde sustentación y aumenta la resistencia, el perfil está en lo que se conoce como pérdida (stall en inglés). 1.3.2. Resistencia La resistencia del perfil se debe a la distribucón de presiones y a los esfuerzos de cortadura en la capa límite, originados por el elevado gradiente de velocidades en ella, que hacen que las velocidades de la corriente se anulen en contacto con la pala. 1.3.3. Coeficientes adimensionales Como queda patente por el análisis anterior, las presiones sobre la pala son proporcionales a la presión dinámica de la corriente sin perturbar, ρv 2 /2, por lo que las fuerzas de sustentación, L, y resistencia, D, así como el par de encabritamiento, M , también lo serán. Además, serán proporcionales a la superficie de la pala, S, resultado del producto de la cuerda, c, por la envergadura, l: S = cl. Se definen los coeficientes adimensionales de sustentación, resistencia y de momento, de la forma siguiente: L tmax cmax , ) c c D tmax cmax CD = 1 2 = CD (α, Re, M, , ) c c 2 ρv S M tmax cmax CM = 1 2 = CL (α, Re, M, , ) c c 2 ρv S CL = 1 2 2 ρv S = CL (α, Re, M, (1.23) Adicionalmente, se define la eficiencia aerodinámica como el cociente CL /CD , que mide la capacidad de proporcionar sustentación frente a la resistencia asociada. Toma valores superiores a 150 en casos ideales e inferiores a 100 en caso prácticos. Se emplea para evaluar las actuaciones de los diferentes perfiles. A continuación se analiza la dependencia de los coeficientes definidos con los diferentes parámetros: Ángulo de ataque: el coeficiente de sustentación, CL , crece de forma aproximadamente lineal con el ángulo de ataque hasta que se produce el desprendimiento, en torno a los 16◦ , punto en el que alcanza un máximo superior a la unidad y cae bruscamente debido a la entrada 8 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos en pérdida. El coeficiente de resistencia, CD , tiene un comportamiento aproximadamente parabólico con un mínimo para un ángulo de ataque del orden de 0,08, disparándose su valor en las proximidades de la situación de pérdida. Por último, el coeficiente de momento, CM , es prácticamente constante. Número de Reynolds: valores elevados de este parámetro adimensional suponen una influencia baja de las fuerzas viscosas en el comportamiento del perfil. A Re elevados aumentan CL /CD , CLmax y CDmax . Sin embargo, el Re en aeroturbinas es relativamente bajo, por lo que ha de tenerse en cuenta su influencia. Número de Mach: salvo en palas muy rápidas y vientos fuertes no suele ser necesario tener en cuenta su efecto sobre el comportamiento del perfil, pudiéndose introducir correcciones. Parámetros geométricos: se consideran generalmente el espesor y la curvatura máximos, como cmax porcentaje de la cuerda, tmax c y c , respectivamente. Los perfiles con espesor bajo tienen un buen comportamiento en un rango pequeño de ángulos de ataque y dan lugar a una pala débil, mientras que un espesor grande resulta en una pala robusta y proporciona una entrada en pérdida gradual, aunque con un elevado coeficiente de resistencia. Normalmente se emplean perfiles con el espesor mínimo compatible con los requerimientos estructurales, la palas suelen tener pequeño espesor en la punta y éste crece hacia la raíz. En cuanto a la curvatura, si ésta es mayor, aumenta el rendimiento, con coeficientes de sustentación y correspondientes ángulos de ataque mayores. 1.3.4. Torbellinos Dado que la pala no se extiende hasta el infinito, el fluido intenta bordearla cerca de la punta para igualar la sobrepresión en el intradós a la succión en el extradós. Este movimiento transversal causa un torbellino libre, que nace en el borde de la pala y es arrastrado por la corriente. Este torbellino hace que el coeficiente de sustentación del perfil del extremo de la pala sea nulo, independientemente de su geometría. La circulación alrededor de un contorno cerrado arbitrario es nula si el flujo es ideal (sin viscosidad ni conductibilidad) e incompresible, denominándose entonces flujo irrotacional. El flujo alrededor de un perfil aerodinámico a elevados números de Reynolds se aproxima a un flujo irrotacional salvo en la capa límite, donde los esfuerzos viscosos dominan. Debido a esto, la sustentación de un perfil es resultado de la circulación embebida en la capa límite o en el interior del propio perfil, situación asociada a un torbellino ligado al perfil. Este torbellino no puede acabar en el fluido, sino que se desprende por los dos extremos de la pala, inyectando rotacionalidad en el fluido. Los torbellinos desprendidos no dan lugar a sustentación, pues son arrastrados por la corriente, uniéndose ambos en el infinito y formando un circuito cerrado. Pueden desprenderse a lo largo de la longitud de la pala si la sustentación a lo largo de ella varía, tienen tendencia a arrollarse sobre sí mismos formando uno más intenso y localizable en la punta o en la raíz de la pala. Los torbellinos justifican la anulación de la sustentación en el final de la pala y resultan coherentes con el bordeo que realiza el fluido patra tratar de igualar las presiones en intradós y extradós. El flujo alrededor de las palas se concibe, según lo anterior, como la superposición de una corriente uniforme y un sistema de torbellinos. Este modelo no considera la viscosidad pero resulta coherente con ella. La viscosidad tiene un efecto de disipación, haciendo que los torbellinos pierdan intensidad y desaparezcan cierta distancia aguas abajo de la turbina. Por otra parte, es causante Inés Mateos Canals 9 CAPÍTULO 1 1.4. TEORÍA DEL ELEMENTO DE PALA de inestabilidades en el sistema de torbellinos, haciendo que pierda la simetría axial. La teoría de sustentación basada en torbellinos provee resultados satisfactorios y conclusiones teóricas útiles, por lo que es frecuentemente empleada en el diseño de aeroturbinas. 1.4. Teoría del elemento de pala Según la teoría del elemento de pala, el área barrida por el rotor se concibe como un conjunto de áreas anulares concéntricas, la acción de la pala es la suma de las acciones independientes en cada uno de estos anillos. La pala se divide longitudinalmente mediante elementos que se consideran independientes entre sí, de forma que el balance entre la tasa de variación del momento del fluido y las fuerzas elementales de la pala puede ser establecido para cada área anular. La teoría básica se debe a Glauert y tiene el objetivo de analizar la variación radial de la carga a lo largo de las palas. Figura 1.5: Teoría del elemento de pala [16] Un elemento de pala supone un corte de la pala a una distancia r y dr del eje. La velocidad relativa al perfil de la corriente incidente, w, resulta de restar a la velocidada del viento local, v, considerada axial, la velocidad de arrastre debida al giro de la pala, vp = Ωr. Entonces, introduciendo las inducciones axiales y tangenciales obtenidas de las teorías de cantidad de movimiento y de momento cinético: w2 = v 2 + vp2 v = v0 (1 − a) vp = Ωr 1 + a0 2 s ⇒w= v02 (1 − a)2 + Ω2 r 2 a0 1+ 2 2 = (1 − a)v0 sin ϕ (1.24) Donde ϕ se denomina ángulo de la corriente, con referencia al plano de giro y queda definido por: (1 − a)v0 1 tan ϕ = = (1.25) a0 a0 (1 + 2 )Ωr (1 + 2 )λ(r) 10 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos y λ(r) es conocida como la rapidez local de la pala: λ(r) = vp Ωr = v (1 − a)v0 (1.26) El motivo por el que se emplea a0 /2 en vez de a0 es la discontinuidad de vθ a través del disco, siendo nula delante cobrando valor no nulo detrás por acción de las palas. Se supone la velocidad tangencial inducida en el perfil como la media de la que existe justo delante y justo detrás. Según lo anterior, las expresiones del impulso, el par y la potencia del elemento de pala son: 1 1 dT = ρw2 (CL cos ϕ + CD sin ϕ)cN dr = ρw2 CN cN dr 2 2 1 2 1 (1.27) dQ = ρw (CL sin ϕ − CD cos ϕ)cN dr = ρw2 CT cN rdr 2 2 1 1 dP = ρw2 (CL sin ϕ + CD cos ϕ)cN dr = ρw2 CN cN Ωrdr 2 2 Donde N representa el número de palas y tanto c como w son funciones del radio, r. Por otro lado, CN y CT son los coeficientes de fuerza normal y en el plano del disco, respectivamente, que también son función del radio, r. Figura 1.6: Velocidades en un elemento de pala [15] 1.5. Teoría del momento sobre elemento de pala (BEM) En inglés, Blade Element Momentum Theory, también conocida como BEM Theory. Es la teoría más ampliamente empleada en métodos de diseño práctico y en computación para predecir cargas y actuación de aerogeneradores. En relación con las técnicas de modelización, está muy ligado a las líneas de investigación actuales, especialmente a la mecánica de fluidos computacional (CFD, del ingés Computational Fluid Dynamics). Esta teoría persigue expresar las acciones aerodinámicas en una pala en función del comportamiento de los perfiles que la forman y hacerlo coincidir con el comportamiento global obtenido mediante las teorías anteriores. La teoría BEM combina las ecuaciones de la Teoría del momento cinético y la Teoría del elemento de pala. Inés Mateos Canals 11 CAPÍTULO 1 1.5. TEORÍA DEL MOMENTO SOBRE ELEMENTO DE PALA (BEM) Las ecuaciones obtenidas a partir de la teoría del momento cinético, considerando el impulso como la tasa de variación del momento lineal del flujo que atraviesa un anillo de radio r y anchura dr son: Impulso dT = 4πρrv 2 a(1 − a)F dr (1.28) y, análogamente en el caso del par como tasa de variación del momento angular: Par dQ = 4πρr3 va0 ω(1 − a)F dr (1.29) En las ecuaciones anteriores se ha hecho uso del factor de efecto punta, F , que se discutirá más adelante en este texto. Por otro lado, se define la solidez como: σ= Nc 2πr (1.30) Igualando la ecuación 1.28 con la correspondiente al impulso obtenida mediante la teoría del elemento de pala, en la ecuación 1.27, se obtiene: a σ(CL + CD tan ϕ) = 1−a 4F tan ϕ sin ϕ (1.31) Análogamente, igualando la ecuación 1.29 con la correspondiente al par obtenida mediante a la teoría del elemento de pala, en la ecuación 1.27: a0 σ(CL tan ϕ − CD ) = 0 1+a 4F sin ϕ (1.32) Normalmente, un proceso iterativo se emplea para resolver las ecuaciones anteriores para cada elemento de pala de anchura ∆r. A partir de estos resultados, la integración de las ecuaciones diferenciales de impulso y el par de 1.27 permite conocer el impulso y el par totales del rotor. 1.5.1. Pérdidas en la punta de la pala En la punta de la pala se producen pérdidas. Estas pérdidas se tienen en cuenta en la teoría BEM por medio de un factor de corrección, que toma valores comprendidos entre 0 y 1 y caracteriza la reducción de las fuerzas a lo largo de la pala. Un método aproximado para estimar el efecto de las pérdidas de punta de pala fue desarrollado por L. Prandtl, que obtuvo la siguiente expresión del factor de pérdidas en punta: ( " −(B/2)[1 − 2 F = cos−1 exp π ( Rr ) sin ϕ r R] #) (1.33) De acuerdo con esta corrección, para un determinado perfil aerodinámico y un ratio punta-velocidad específico según la longitud de la pala, la geometría de la pala puede diseñarse para la optimización del rotor. Mediante estos parámetros geométricos es posible analizar el rendimiento del rotor. 12 Escuela Técnica Superior de Ingenieros Industriales (UPM) Capítulo 2 Modelización de aerogeneradores Durante los últimos años, en la mecánica de fluidos computacional (en inglés, Computational Fluid Dinamics, CFD) ha adquirido gran importancia la resolución de las ecuaciones de NavierStokes para fluidos incompresibles. El flujo de en aerogeneradores es incompresible a velocidades entre 5 y 25 m/s. Los efectos de compresibilidad están presentes solo en la punta de la pala, debido a las elevadas velocidades causadas por la rotación. En cualquier caso, el número de Mach en la punta sigue siendo relativamente bajo (inferior a 0.25), por lo que se puede considerar flujo incompresible en la resolución de las palas. Las ecuaciones de Navier-Stokes para flujo incompresible, que son aplicables a la modelización de aeroturbinas, son: ∇·u=0 (2.1) ∂u 1 + (u · ∇)u = − ∇p + v∇2 u ∂t ρ (2.2) En el análisis de flujos de la capa límite atmosférica (ABL, del ingles Atmospheric Boundary Layer), se suele emplear la aproximación de Boussinesq y es necasaria la resolución de una ecuación de temperatura. El efecto de Coriolis puede ser relevante en el estudio de grandes aerogeneradores y parques eólicos. Para flujos en la ABL, las escalas varían entre 1 km y 1 mm, disparidad que hace imposible la utilización de simulaciones numéricas directas (DNS, del inglés Direct Numerical Simulation) dados los elevados requisitos computacionales que éstas presentan. En la investigación, para resolver los flujos medios, se han adoptado métodos basados en las ecuaciones de Navier-Stokes promediadas por Reynolds (Reynolds-Averaged Navier-Stokes, RANS): ∂ ū 1 + (ū · ∇)ū = − ∇p + ν∇2 ū − ∇ · (u0 u0 ) ∂t ρ (2.3) Donde u = ū + u0 es la descomposición de Reynolds y u0 u0 es el tensor de Reynolds, que es modelizado en términos de las magnitudes medias del flujo. Algunos investigadores de la comunidad de la energía eólica han adoptado la formulación LES (Large Eddy Simulations, que en español se traduciría como simulaciones de grandes torbellinos). El flujo en los aerogeneradores es altamente cambiante y, con las ecuaciones RANS, estas fluctuaciones se pierden dado el uso de magnitudes medias. El uso de LES proporciona una aproximación más adecuada, que resuelve la dependencia temporal del flujos. Como su propio nombre indica, LES calcula los torbellinos grandes, mientras que los más pequeños son modelados. Las ecuaciones de Inés Mateos Canals 13 CAPÍTULO 2 2.1. SIMULACIÓN DE LA CAPA LÍMITE ATMOSFÉRICA Navier-Stokes restringidas empleadas son: ∂ ū 1 f − ũũ) + (ū · ∇)ū = − ∇p̃ + ν∇2 ū − ∇ · (uu ∂t ρ (2.4) f − ũũ) suele modelizarse mediante: Donde el término (uu f − ũũ) = −νSGS (∇ũ + (∇ũ)T ) τSGS = (uu (2.5) νSGS es la viscosidad de escala subred (subgrid-scale) y la operación de restricción se aplica en u = ũ + u00 . Durante los últimos años, RANS y LES se han empleado para simular el flujo en aerogeneradores. Sin embargo, en investigación se está popularizando el uso de LES porque permite una mejor predicción de la capa límite atmosférica y de la estela tras las turbinas. 2.1. Simulación de la capa límite atmosférica Los factores que afectan en mayor medida el comportamiento de la capa límite atmosférica (ABL) son el viento geostrófico, la aspereza del terreno, el efecto Coriolis y la estratificación térmica, que puede ser categorizada como estable, inestable o neutra. En particular, la estratificación térmica inestable ocurre cuando existe un calentamiento significativo de la superficie terrestre, que causa la elevación del aire caliente cercano a la misma. El ascenso del volumen de aire supone una caída de la presión sobre él, su expansión y enfriamiento casi adiabático, de forma que el resultado es una mayor mezcla vertical y una ABL de mayor grosor. Por el contrario, un terreno frío elimina el desplazamiento vertical del aire, dando lugar a una estratificación térmica estable; por último, la estratificación neutra se asocia al equilibrio térmico del viento con el aire circundante en su ascenso. La simulación de la estela tras un aerogenerador conlleva el inclusión en los modelos de unas condiciones de contorno que permitan imitar las características de la ABL. En los modelos se consideran el perfil de velocidad transversal, la anisotropía de la turbulencia y la naturaleza no estacionaria del flujo. La superficie inferior se suele describir como una condición de contorno sin desplazamiento o como una función de pared (en inglés, wall function) para considerar la aspereza del terreno. En el límite superior, se suelen emplear condiciones de flujo libre, mientras que en la dirección de la corriente, se recurre habitualmente a la periodicidad para simular parques eólicos infinitos o como precursor en la generación de la turbulencia de la ABL. Varios estudios han tratado la capa límite atmosférica en presencia de aerogeneradores mediante simulaciones de grandes torbellinos LES (Porté-Agel [18], Calaf [19]). La resolución espacial empleada no permite resolver por completo la geometría de las palas del aerogenerador y el flujo en la ABL en torno a ellas dadas las limitaciones de los medios de computación disponibles. Esta restricción es una de las razones impulsoras para el desarrollo de modelos que buscan simular aerogeneradores sin que ello implique la solución completa de su geometría. 2.2. Modelo del disco actuador (ADM) ADM (Actuator Disk Model) simula la turbina como un disco en el seno del fluido que impone una serie de fuerzas sobre él y se basa en la resolución de las ecuaciones de Navier-Stokes. Existen 14 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos numerosas versiones de este modelo, algunas aplican fuerzas axial (empuje, T ) y tangencial, mientras que otras aplican sólo empuje. Es un modelo frecuentemente utilizado en investigación para simular turbinas y estelas, a pesar de no crear los vórtices de punta de pala que llegan a la estela. El modelo sin fuerza tangencial se basa en la teoría de la cantidad de movimiento para vuelo axial. La fuerza del disco actuador se distribuye uniformemente sobre el disco circular, de forma que la fuerza ejercida sobre el disco es: 1 T = ρũ20 ACT (2.6) 2 Donde T es la fuerza de empuje, ρ es la densidad, ũ0 es la velocidad no perturbada del flujo axial incidente sobre el centro del disco, A es el área barrida por las palas del rotor y CT es el coeficiente de empuje. En algún proyecto de investigación se ha empleado la velocidad promediada sobre la región del disco en vez de la velocidad en el centro del rotor. Esta versión de ADM constituye una aproximación unidimensional a la fuerza de empuje y no considera la rotación del flujo que se produce como consecuencia del movimiento de las palas, tampoco tiene en cuenta la dependencia radial del empuje. Por otro lado, el ADM que sí tiene en cuenta la fuerza tangencial se implementa por integración de las fuerzas de resistencia y sustentación de las palas según la resolución espacial y temporal de las LES. La fuerza que ejerce el disco sería: f2D = dF 1 2 Nc = ρvrel (CL eL + CD eD ) dA 2 2πr (2.7) Donde N es el número de palas, vrel es la velocidad local, CL y CD son los coeficientes de sustentación y resistencia, respectivamente, c es la longitud de la cuerda y eL y eD son vectores unitarios en la dirección de sustentación y resistencia, respectivamente. Como se ha visto, existen varias versiones de ADM, a las que pueden incorporarse correcciones para tener en cuenta las pérdidas en el extremo de la pala o la conicidad del rotor. El modelo proporciona resultados generales razonables, sin embargo, no tiene en cuenta todas las características del flujo. Una de sus mayores desventajas es la utilización de fuerzas promediadas sobre la superficie barrida por las palas, mientras que las fuerzas reales actúan de forma diferente en cada localización instantánea de la pala. Por otro lado, ADM no considera los vórtices en la punta de la pala, que son relevantes en el estudio de la estela próxima de un aerogenerador. Dadas estas limitaciones, este modelo está siendo sustituido en investigación por modelos que reflejen con mayor precisión los fenómenos que tienen lugar en el flujo de una corriente de aire con un aerogenerador en su trayectoria. 2.3. Modelo de la linea actuadora (ALM) En este caso, las palas del aerogenerador se modelizan como una serie de elementos de pala a lo largo de ellas. El modelo emplea datos de los perfiles aerodinámicos para calcular la sustentación y la resistencia en cada elemento actuador y ejerce estas fuerzas sobre el flujo. El ALM no resuelve la geometría de las palas por completo. Inés Mateos Canals 15 CAPÍTULO 2 2.4. MODELO DE LA SUPERFICIE ACTUADORA (ASM) Las fuerzas son calculadas en base a las características de los perfiles aerodinámicos, de forma que la fuerza total sobre un elemento de pala es: f2D = dF 1 2 = vrel c(CL eL + CD eD ) dA 2 (2.8) donde CL y CD son los coeficientes de sustentación y resistencia, respectivamente, vrel es la velocidad local del viento y eL y eD son vectores unitarios con la dirección de la sustentación y la resistencia, respctivamente. Las fuerzas sobre los elementos de pala son proyectadas mediante la siguiente convolución: f = f ⊗ η (2.9) donde η es una función de regularización. Se suele emplear la función gaussiana: η = 1 3 π 3/2 exp[−(r/)2 ] (2.10) donde establece la amplitud de la función. El modelo fue validado por comparación con los datos de la turbina Nordtank. El ALM proporciona resultados precisos pero tiene limitaciones. Cada fuerza elemental es la fuerza total ejercida sobre un elemento de pala, mientras que las fuerzas reales están distribuidas uniformemente sobre la cuerda del elemento de pala. Para distribuir el punto de aplicación de la fuerza, se hace uso de la función gaussiana, método que elimina las inestabilidades numéricas pero que no reproduce exactamente la distribución real de la fuerza aplicada. El uso de otro tipo de funciones para proyectar la fuerza sobre la superficie del elemento de pala es posible y objeto de estudio en investigación. 2.4. Modelo de la superficie actuadora (ASM) Se trata de un modelo reciente que calcula las fuerzas sobre un perfil aerodinámico en 2D como función de la cuerda, evitando de esta forma una de las mayores simplificaciones del ALM. Esto se consigue mediante fórmulas empíricas en las que las fuerzas son función de la posición de la cuerda. El ASM supone un enfoque prometedor en la modelización de aerogeneradores para la simulación mediante CFD. La mayor ventaja del modelo es el carácter más realista de la distribución de fuerzas, mientras que una de sus limitaciones es la dependencia de cálculos externos de la distribución de fuerzas en función de la posición a lo largo de la cuerda del perfil aerodinámico. 2.5. Comparación del ADM y el ALM Como se ha mencionado en secciones anteriores, existen dos versiones del modelo de disco actuador, una que ejerce sobre el flujo empuje y fuerzas tangenciales (ADM-R) y otra que impone sólo empuje (AMD-NR). Por otr lado, el modelo de línea actuadora (ALM) calcula las fuerzas de sustentación y resistencia para secciones discretas de las palas y las distribuye a lo largo de las mismas. El ALM, a diferencia del ADM, tiene en cuenta los vórtices en la punta de la pala y las 16 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos estructuras helicoidales en la zona de la estela próxima al aerogenerador. Las fuerzas de estos modelos son proyectadas suavemente, por medio de una convolución con una función gaussiana, sobre la corriente de aire para evitar puntos singulares e inestabilidad numérica. Además, la góndola y la torre se modelizan como fuerzas resistivas e incluyen sus respectivos coeficientes de resistencia. Los modelos numéricos han sido comparados con resultados experimentales de mediciones realizadas en túneles de viento en el Laboratorio de Saint Anthony Falls, donde se caracterizaron estelas de aerogeneradores por la velocidad e intensidad de turbulencia del flujo. Las predicciones del ALM y el ADM-R fueron corroboradas por los datos del experimento, mientras que la predicción de la distribución de velocidades del ADM-NR resultó inexacta.En cuanto a la intensidad de la turbulencia y la zona en la que ésta se maximiza, tanto el ALM como el ADM-R demostraron proporcionar resultados correctos. 2.6. Implementación de modelos Tanto el modelo del disco actuador (ADM) como el de línea actuadora (ALM) predicen las fuerzas de las palas según la velocidad local del fluido en cada elemento actuador. Estas fuerzas son proyectadas sobre el campo fluido e impuestas en las ecuaciones de Navier-Stokes. Las fuerzas calculadas por los modelos tienen carácter puntual y son ejercidas individualmente por cada elemento de pala. Para distribuir las fuerzas puntuales por el volumen circundante, se proyectan mediante una ecuación gaussiana: F fturb = 3 3/2 exp[−(r/)2 ] (2.11) π donde fturb es la fuerza proyectada por unidad de volumen, F es la fuerza de un elemento actuador, establece la magnitud de la proyección y r es la distancia al punto central del elemento actuador. Ambos modelos permiten determinar datos de la pala como el ángulo de ataque, la sustentación y la resistencia a través de la velocidad local en cada elemento actuador. El ángulo de ataque local se define como el ángulo entre la cuerda y la línea de incidencia del flujo local. Este supuesto se justifica dado que el punto actuador se encuentra en el centro del torbellino ligado a él, donde las perturbaciones del flujo a su paso por el perfil son despreciables. Los coeficientes de sustentación y resistencia están tabulados para los diferentes perfiles aerodinámicos. La potencia extraída se calcula por integración del par creado respecto al centro del rotor, multiplicado por la velocidad angular del rotor. 2.6.1. Parámetro de la función gaussiana La función gaussiana se introduce en el modelo para conseguir una mejor aproximación a la realidad al ejercer la fuerza sobre la corriente incidente. Como consecuencia, las fuerzas aplicadas mediante los elementos actuadores del disco pierden el carácter puntual se parecen más a una fuerza continua impuesta sobre el volumen fluido. La función evita las inestabilidades numéricas que causarían las fuerzas puntuales sobre las secciones discretas de la pala. Inés Mateos Canals 17 CAPÍTULO 2 2.6. IMPLEMENTACIÓN DE MODELOS El parámetro determina la amplitud de la función gaussiana, su papel es muy relevante en los modelos de disco y linea actuadores. Simulaciones llevadas a cabo en un estudio del Laboratorio Nacional de Energías Renovables de EE.UU. (NREL) [20] mostraron que la potencia que predicen ambos modelos es menor cuanto menor sea el valor de . El estudio permitió observar que un valor de pequeño en relación al tamaño de la red conduce a inestabilidades numéricas en la resolución de la ecuaciones del flujo. Análogamente, se constató que un valor de igual al tamaño de la red no proporciona resultados convergentes. Para evitar inestabilidades y divergencias en la resolución de las ecuaciones, se emplean valores del parámetro por encima de cierto umbral, que se establece como el doble del tamaño del mallado, como fue sugerido por Troldborg [21] Sin embargo, las simulaciones con valores elevados de frente al tamaño de la red, a pesar de dar soluciones convergentes, dan lugar a valores de potencia que se alejan de las predicciones de la Teoría BEM. El incremento de supone un ensanchamiento de la capa de cizalladura en torno a la estela, disminuyendo las turbulencias del flujo. El estudio concluyó, según lo anterior, que existe un valor óptimo de este parámetro relacionado con la geometría de la pala. Este valor proporcionaría las mejores aproximaciones de potencia generada en las simulaciones que emplean tanto ADM como ALM, es decir, valores muy cercanos a los que predice la teoría BEM. La implementación de diferentes valores de en cada elemento de pala, teniendo en cuenta la variación de la cuerda a lo largo de la misma, proporcionaría una mejor aproximación de la distribución de la fuerza sobre su superficie. 2.6.2. Implementación del ALM La geometría de la turbina se define en un sistema de coordenadas cartesianas, basada en tres localizaciones puntuales: góndola, buje y suelo. Las palas están constituidas por un conjunto de puntos a lo largo de su eje y cada de estos puntos es el centro de un elemento actuador. La primera pala se define paralela al plano del suelo y el resto son definidas mediante una matriz de rotación: u2x + (1 − ux2 )c ux uy (1 − c) − uz s ux uz (1 − c) + uy s x0 h i u2y + (1 − uy 2 )c uy uz (1 − c) − ux s y0 = x1 y1 z1 ux uy (1 − c) + uz s z0 ux uz (1 − c) − uy s uy uz (1 − c) + ux s u2z + (1 − uz 2 )c (2.12) donde c = cos θ y s = sin θ, siendo θ el ángulo de rotación, u = ux ex + uy ey + uz ez es el vector unitario normal al plano de rotación, x0 , y0 , z0 las coordenadas del punto a rotar y x1 , y1 , z1 las coordenadas del punto rotado. Las palas giran cada intervalo de tiempo según la velocidad angular del rotor. La localización de los puntos a lo largo de cada pala se rota mediante la matriz de rotación anterior. Los puntos actuadores alineados definen la localización de las secciones de las palas que tienen unas componentes de sustentación y resistencia determinadas, según el perfil aerodinámico, el ángulo de giro, la longitud de la cuerda y las características del flujo incidente. Un vector de fuerzas de sustentación y resistencia es asignado a cada punto de la línea de actuación. Una fuerza de igual magnitud y dirección y sentido impuesto se impone en la ecuación de la cantidad de movimiento. La sustentación y la resistencia se calculan en base a la velocidad local del viento y el ángulo de ataque: 1 L = CL (α)ρv 2 cw 2 18 (2.13) Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos 1 D = CD (α)ρv 2 cw (2.14) 2 Ecuaciones en las que α es el ángulo de ataque respecto a la cuerda, CL y CD son los coeficientes de sustentación y resistencia, respectivamente, ρ es la densidad, v es la velocidad local del viento, c es la cuerda y w es el ancho de la sección. Como se mencionó con anterioridad, las fuerzas se proyectan sobre el flujo mediante la función Gaussiana, 2.11. El plano del vector de fuerza es perpendicular a la las líneas de actuación. El producto vectorial de un vector con la dirección de la linea de elementos actuadores, ~z, y otro con la dirección del eje en torno al que giran las palas, proporciona un vector ~y , cuya dirección es tangencial a la rotación de las palas. El producto vectorial de este vector tangencial ~y con ~z resulta en el vector de dirección ~x. Los vectores ~x y ~y definen el plano en el que se situan los vectores de fuerza. Un nuevo conjunto de vectores fuerza es definido en cada instante para cada pala, según la rotación de las palas. El vector velocidad es proyectado sobre el nuevo plano mediante el producto escalar. 2.6.3. Implementación del ADM Como se ha mencionado en secciones anteriores, el ADM simula la turbina como un disco en el seno del fluido que ocupa el área barrida por las palas del aerogenerador. Esto se consigue mediante la división del área en un gran número de elementos. La diferencia en la geometría del ADM respecto al ALM reside en que las lineas de actuación son secciones circunferenciales. Cada punto representa un elemento del disco. La fuerza en cada elemento se escala mediante el factor de solidez, definido como la razón del área ocupada por las palas y el área total del disco: σ= N AP AR (2.15) Donde N es el número de palas, AP es el área de cada sección y AR es el área barrida por las palas. Las fuerzas de sustentación y resistencia se definen a continuación: 2.6.4. 1 L = σ CL (α)ρv 2 cw 2 (2.16) 1 D = σ CD (α)ρv 2 cw 2 (2.17) Método numérico En la resolución de los modelos mencionados, tanto ALM como ADM, es necesario emplear una herramienta de mecánica de fluidos computacional (CFD). El código LES resuelve las ecuaciones de Navier-Stokes para fluidos incompresibles. El flujo alrededor de las turbinas es incompresible, con velocidades entre 5 y 25 m/s, por lo que las ecuaciones aplicables a la modelización de turbinas son: Inés Mateos Canals ∇ · ũ = 0 (2.18) ∂ ũ 1 + (ũ · ∇)ũ = − ∇p̃ + ν∇2 ũ − ∇ · τSGS + fturb ∂t ρ (2.19) 19 CAPÍTULO 2 2.6. IMPLEMENTACIÓN DE MODELOS donde ũ y p̃ son la velocidad y la presión restringidas, respectivamente, ν es la viscosidad dinámica, fturb es la fuerza de la turbina, ∇ · τSGS = −νSGS (∇ũ + (∇ũ)T ) (2.20) y νSGS es la viscosidad a escala subred (subgrid-scale) y es modelada mediante el modelo a escala subred de Smagorinsky. 20 Escuela Técnica Superior de Ingenieros Industriales (UPM) Capítulo 3 Simulación con IRIS IRIS es un software genérico de elementos finitos desarrollado por el Departamento de Mecánica Estructural y Construcciones Industriales de la Escuela Técnica Superior de Ingenieros Industriales de Madrid. IRIS se emplea para resolver el problema incompresibilidad por medio del método de elementos finitos estabilizado por SUPG. 3.1. Método de elementos finitos estabilizado La resolución del problema de flujo incompresible por métodos de elementos finitos cuenta con dos potenciales fuentes de inestabilidades asociadas a la formulación de Galerkin del problema. Una de las fuentes se debe a la presencia de términos de advección en las ecuaciones, que puede resultar en falsas oscilaciones nodo a nodo, principalmente en el campo de velocidades. Dichas oscilaciones se vuelven más aparentes en flujos de desplazamiento horizontal o advección dominante (como es el caso de flujos con elevados números de Reynolds). La otra fuente de inestabilidades se encuentra en el uso de combinaciones inadecuadas de funciones de interpolación para representar los campos de presión y velocidad. Estas inestabilidades aparecen generalmente como oscilaciones en el campo de presiones. De hecho, las mencionadas inestabilidades numéricas no son realmente inherentes a la formulación de elementos finitos y aparecen también las versiones estándar de otras técnicas de discretización como los métodos de diferencias finitas y volúmenes finitos. Las formulaciones de elementos finitos estabilizadas han sido diseñadas para evitar posibles inestabilidades numéricas como las que acaban de ser descritas. La técnica de estabilización empleada por IRIS se conoce como SUPG, sigla en inglés de streamline-upwind/Petrov-Galerkin. El método SUPG consigue la estabilización de problemas de flujo incompresible mediante la adición a la formulación de Galerkin de una serie de términos. Estos términos implican el producto del residuo de la ecuación de la cantidad de movimiento por el operador de advección que actúa en la función de forma. La formulación anterior fue presentada por Hughes y Brooks [22]. La estabilización SUPG de la formulación de la función de flujo de vorticidad para problemas de flujo incompresible, incluyendo aquellos con dominios múltiplemente conectados, se debe a Tezduyar et al. [23]. Inés Mateos Canals 21 CAPÍTULO 3 3.2. 3.2. MODELO ADM Modelo ADM El modelo escogido para simular el flujo de corriente en torno a un aerogenerador en este proyecto es el modelo del disco actuador (ADM) axisimétrico. Se considera el rotor como un disco cuya superficie coincide con el área barrida por las palas del aerogenerador. Según la Teoría del momento sobre elemento de pala (BEM), la carga sobre las palas del rotor varía radialmente a lo largo de las mismas y se consideran, por tanto, tubos de corriente anulares con independencia radial entre sí. Además, se asume la inexistencia de fuerzas de presión sobre el volumen de control y que el flujo presenta simetría axial. Las fuerzas aerodinámicas que actúan sobre el rotor están gobernadas por las velocidades locales y se determinan por medio de las características de los perfiles aerodinámicos en 2D. La velocidad relativa a al perfil aerodinámico elemental se determina según: 2 vrel = (v0 − wy )2 + (Ωr + wθ )2 (3.1) donde Ω es la velocidad angular. Se define φ como el ángulo formado por la vrel y la dirección de incidencia axial, de forma que se puede calcular según la siguiente expresión: φ = tan−1 v o − wz Ωr + wθ (3.2) Localmente, el ángulo de ataque es función del ángulo φ anterior y el ángulo de paso, γ: α = φ−γ. Según esto, las fuerzas de sustentación y resistencia por elemento de pala en función de los datos de los perfiles aerodinámicos: 1 2 (L, D) = ρvrel cB(CL eL , CD eD ) 2 (3.3) donde CL y CD , los coeficientes de sustentación y resistencia, respectivamente, son función del ángulo de ataque, α y del número de Reynolds, Re. B es el número de palas y c es la cuerda. La fuerza sobre un elemento de pala es la resultante de las dos fuerzas aerodinámicas anteriores: F=L+D (3.4) La proyección de la fuerza resultante F sobre las direcciones axial y tangencial del rotor proporciona las componentes: Fz = L cos φ + D sin φ Fθ = L sin φ − D cos φ ) ⇒ F = Fθ + Fz = Fθ eθ + Fz ez (3.5) Estas son las fuerzas sobre las palas que equilibran los cambios de cantidad de movimiento en las direcciones axial y tangencial, por lo que: ∆T = Fz ∆r; 22 ∆Q = Fθ r∆r (3.6) Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Figura 3.1: Fuerzas [24] 3.3. Implementación del modelo en IRIS La implementación del modelo ADM se ha llevado a cabo mediante el uso conjunto de MATLAB e IRIS. A continuación se detallan el proceso seguido y las consideraciones tomadas, así como las simplificaciones que han sido necesarias. 3.3.1. Simplificaciones La simulación de un parque eólico supone la generación de mallados de dominios relativamente extensos que comprendan el conjunto completo de aerogeneradores que lo integran. Con vistas a posibilitar trabajos futuros de optimización de un parque eólico, la modelización empleada en el presente trabajo se centra en desarrollar un modelo que permita, simultáneamente, aproximar adecuadamente el mecanismo de extracción de potencia que tiene lugar en un aerogenerador y limitar las exigencias en cuanto a capacidad y tiempo de computación. Con este objetivo en mente, se han llevado a cabo una serie de simplificaciones, que se describen a continuación. La modelización empleada en las simulaciones está restringida al rotor de un aerogenerador o a un conjunto de rotores agrupados. Esto implica que el modelo no considera la torre del aerogenerador y se centra en la dinámica de las palas del rotor, en la zona del espacio barrida por ellas al girar. La torre tiene cierto efecto sobre la carga aerodinámica de las palas, efecto que se manifiesta en forma de cambios bruscos de la misma, así como de la deflexión de la pala. Estas variaciones son producidas por el paso de las palas por delante de la torre, que altera la circulación del aire. Se considera en este estudio que las variaciones en la carga aerodinámica no son suficientemente significativas como para alterar la generación de potencia en la zona del disco afectada. Tampoco se considera en la modelización realizada la influencia de la topografía del terreno sobre el campo del fluido. Esta hipótesis alejaría al modelo de la realidad de forma drástica en el caso de parques eólicos situados en zonas con importantes diferencias de relieve. En estos parques, los ejes de los rotores no se encontrarían en el mismo plano horizontal, por lo que la interacción entre las estelas sería diferente a la del modelo y éste no resultaría valido. En la simulación de zonas Inés Mateos Canals 23 CAPÍTULO 3 3.3. IMPLEMENTACIÓN DEL MODELO EN IRIS con poco relieve, las consecuencias de la simplificación serían poco importantes. Otra simplificación importante que se está tomando en el modelo escogido es el obviar los efectos del suelo sobre la distribución de velocidades en una sección vertical de la corriente de aire. El suelo frena el viento en la zona próxima a él. Esta simplificación es importante y sólo se justifica en cierta medida dada la altura de las torres de los aerogeneradores similares, en torno a los 90 m. Se supone, en este caso, que la distancia mínima posible de la punta de la pala a la zona de la corriente afectada por la presencia del suelo es suficientemente grande como para despreciar el efecto. No se tienen en cuenta los efectos de la capa límite en las condiciones de contorno. Tampoco se consideran las variaciones en las condiciones del aire que tienen lugar en la capa límite atmosférica. Las dos consideraciones anteriores sobrepasan el alcance de este proyecto y exceden la capacidad de los recursos computacionales disponibles. Las turbulencias del flujo producen variaciones en la generación de potencia de los aerogeneradores. La influencia sobre la potencia generada de los flujos no uniformes es complicada, dado que implica que cada sección de la pala está sujeta a cargas aerodinámicas espacial y temporalmente variables. Debido a esto, la simulación de la turbulencia supone un reto de la modelización del flujo en torno a un aerogenerador, dada la elevada capacidad de computación que ésta exige. Los medios de los que se dispone y el objetivo a largo plazo de este proyecto obligan a descartar modelos que reproduzcan con fidelidad los efectos de la turbulencia. En lo que respecta concretamente a la modelización de un aerogenerador, se ha dado un paso más en la simplificación del modelo ADM axisimétrico. La simetría axial de la corriente incidente respecto al eje giro del rotor permite no considerar la componente acimutal, es decir, hace que el problema sea bidimensional (componentes axial y radial en coordenadas cilíndricas). En la implementación del modelo en IRIS se hace uso de un sistema de coordenadas cartesianas en vez de cilíndricas, lo que permite que el programa simule dominios que incluyan parques eólicos y no sólo un aerogenerador aislado. En el cambio a cartesianas de la ecuación de Navier-Stokes de la cantidad de movimiento, se desprecia el término fuente debido a la geometría. 3.3.2. Elementos actuadores Según el modelo del disco actuador, las fuerzas que se ejercen sobre el fluido se aplican sobre elementos distribuidos uniformemente sobre la superficie del rotor. El número de elementos actuadores influye sobre la magnitud de la potencia calculada y sobre la continuidad de la fuerza integrada que ejerce el disco en la totalidad de su superficie. Se ha empleado MATLAB para definir la disposición de los elementos actuadores a lo largo del diámetro del disco. La fuerza sobre el fluido en sentido axial que se implementa por medio de elementos actuadores es la siguiente: F = Fθ 3 π 3/2 exp[−(r/)2 ] (3.7) donde, 1 2 Fθ = ρvrel cB(CL sin φ − CD cos φ) 2 24 (3.8) Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Mediante la convolución con una función gaussiana se tiene en cuenta la distribución de la fuerza aplicada sobre toda la superficie del elemento de pala, de forma que ésta no se aplique de forma puntual sobre el volumen fluido. Se puede definir una constante K de la siguiente forma: 2 Fθ = Kρvrel ; 1 K = cB(CL sin φ − CD cos φ) 2 (3.9) La cuerda es característica del perfil aerodinámico y varía radialmente, es decir, a lo largo de la pala. El valor correspondiente a cada elemento actuador se obtiene por interpolación a partir de la tabla de propiedades aerodinámicas de las palas del aerogenerador. La constante anterior agrupa una serie de términos de la fuerza ejercida por cada elemento actuador, entre los que se encuentra la cuerda. El valor de K se asigna mediante MATLAB, mientras que la densidad y la velocidad relativa son implementadas por medio de IRIS, de ahí la agrupación. Se ha trabajado con una densidad del aire de valor ρ = 1,205kg/m3 . Por otro lado, la velocidad relativa, vrel se define como: 2 vrel = vr2 + u2 = (Ωr)2 + u2 (3.10) es decir, la velocidad relativa es la combinación de la velocidad del flujo incidente, u, y la velocidad acimutal impuesta por el giro del las palas del rotor, vr . La velocidad del fluido es diferente en cada elemento de pala, lo mismo ocurre en el caso de la velocidad acimutal, que varía según la distancia del centro del elemento de pala al centro del rotor. Para considerar la variación radial de la componente acimutal de la velocidad, IRIS hace uso del parámetro relación de velocidad de punta, que se define como la razón entre la velocidad acimutal en la punta de la pala, vR y la velocidad local del flujo incidente, u, es decir: λ= vR ΩR = u u (3.11) El coeficiente de potencia de una turbina alcanza el valor máximo con cierto valor de la relación de velocidad de punta λ. Este es el valor que da el mejor rendimiento y que ha sido empleado en las simulaciones. Según lo anterior, la expresión de la velocidad relativa empleada por IRIS es la que sigue a continuación: " 2 # r 2 r 2 λr 2 2 2 2 vrel = ΩR · + u = λu +u =u 1+ (3.12) R R R Las carácterísticas del disco que modeliza el rotor y la distribución de elementos actuadores sobre el mismo, así como los valores de la constante K asociada a cada uno de los elementos anteriores, se pasan a IRIS a través de un archivo de texto generado por MATLAB. El archivo tiene un aspecto similar al ejemplo de la figura 3.2, en el caso de un aerogenerador aislado: Inés Mateos Canals 25 CAPÍTULO 3 3.3. IMPLEMENTACIÓN DEL MODELO EN IRIS Figura 3.2: Archivo de texto generado con MATLAB como input para IRIS La primera fila del archivo contiene información sobre el número de aerogeneradores que se simulan en el dominio, el tamaño del mallado, la mitad la distancia entre elementos actuadores y el valor de . En la segunda fila se especifican los detalles del disco actuador: la localización de su centro en el dominio simulado, que en este caso es (x0 , y0 ) = (600, 300), el diámetro (D = 126m), el valor de la relación de velocidad de punta (λ = 7,55), y el número de elementos actuadores en los que se subdivide el disco (419). Las filas siguientes (en este caso, 419) proporcionan información sobre la posición del centro de cada elemento actuador (x, y) y a continuación, el valor de la constante K asociada, con valor negativo por ser la fuerza aplicada de oposición al sentido de la corriente. En simulaciones posteriores, en dominios que incluyen agrupaciones de aerogeneradores, se han empleado archivos de texto similares, con la especificación del número de rotores que integran el conjunto y la información particular de cada uno de ellos. Las características de cada aerogenerador y los elementos actuadores asociados a su disco se suceden de forma ordenada en el archivo de texto. 3.3.3. Ampliación de IRIS Antes de este proyecto, IRIS ya resolvía las ecuaciones de Navier-Stokes para flujo incompresible por el método de elementos finitos estabilizado. Debido a la presencia del aerogenerador, se introduce un término nuevo, de forma que las nuevas ecuaciones de Navier-Stokes que debe resolver IRIS son las siguientes: ∂u + ρ(u · ∇)u − ρν∇2 u = −∇p + ρfturb ∂t ∇·u=0 ρ (3.13) El término asociado a la fuerza de la turbina en la ecuación de la cantidad de movimiento supone la aparición del siguiente término en la resolución de las ecuaciones: ρfturb N A 26 (3.14) Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos así como los términos debidos a la estabilización: τ ρfturb (∇N A ) · u (3.15) τ ρ(fturb · ∇N A ) El algoritmo iterativo de Newton-Raphson, empleado por IRIS para resolver el sistema linealizado, precisa de las derivadas de los términos anteriores, cuyo desarrollo matemático se presenta a continuación: fturb fD < u, u > 0 = ; 0 Ke 3 π 3/2 fD = K = exp[−(r/)2 ] ∂(−ρfturb N A ) ∂fturb N A ∂N A ∂fturb = −ρ = −ρ fturb ⊗ + NA ∂u ∂u ∂u ∂u ∂ < u, u > = ∂u ∂u ∂u T u+ ∂u ∂u T u=2 ∂u ∂u T ! (3.17) <u,u> u ∂u 0 = Σ2KN B 0 (3.18) < u, u > < u, u > ∂ ∂ ∂fturb 0 0 = K = K =K ∂u ∂u ∂u 0 0 (3.16) 0 0 T u ∂ u = 2 ΣN B 0 u = 2ΣN B IT u ∂u 0 = Σ2N B IT u = Σ2N B u (3.19) " ∂ ∂ (∇N A · u) ∂fturb [−τ ρfturb (∇N A · u)] = −τ ρ [fturb (∇N A · u)] = −τ ρ fturb ⊗ + (∇N A · u) ∂u ∂u ∂u ∂u # u B A B = −τ ρ fturb ⊗ Σ∇N u + (∇N · u)Σ2KN 0 0 u = −Στ ρN B fturb ⊗ ∇N A + (∇N A · u)2K 0 0 (∇N A · u) = ∂u ∂u ∂u T A ∇N + ∇N A ∂ fturb [−τ ρfturb ∇N A ] = −τ ρ ∂u ∂u h −τ ρ Σ2KN B [u Inés Mateos Canals 0 ∂∇N A ∂u = −τ ρ i (3.20) !T u = ΣN B IT ∇N B ∇N A ∂∇N A ∂u !T fturb + 0]∇N A = −Σ2τ ρKN B [u 0 ∂fturb ∂u 0]∇N A T (3.21) ∇N A (3.22) 27 Capítulo 4 Simulaciones realizadas En este capítulo se detalla el proceso seguido en la modelización y simulación del comportamiento de una corriente de aire que incide sobre un aerogenerador, dadas una determinadas condiciones de operación. A continuación se define el aerogenerador que se ha empleado en la simulación. 4.1. Aerogenerador NREL de 5 MW En las simulaciones de este proyecto se han empleado las dimensiones y especificaciones de una turbina real y estandarizada, de forma que los resultados resultaran aplicables a un caso real en la mayor medida posible. La turbina que ha sido escogida se trata de aerogenerador de gran tamaño, con una potencia generada de varios megavatios y resulta representativo del grupo de turbinas marítimas. La elección del tamaño y la potencia de la turbina, que pueden considerarse elevados, se justifica debido a los altos costes de instalación y soporte de un parque eólico en alta mar. La potencia generada por las aeroturbinas debe ser suficientemente alta como para justificar la inversión y resultar rentable. Teniendo en cuenta lo anterior, el rango de potencias de las turbinas que entran dentro de esta categoría varía de 5 a 20 MW. En estudios y proyectos de investigación similares, las turbinas empleadas generan en torno a 5 o 6 MW, por lo que la turbina simulada en este proyecto también es de 5 MW. El aerogenerador escogido es el modelo estándar de aerogenerador marítimo desarrollado por el Laboratorio Nacional de Energía Renovable de EE.UU. (NREL), que comparte las especificaciones fundamentales con el modelo REpower 5M. Sin embargo, al tratarse este último de un prototipo, la información disponible públicamente sobre la turbina resulta incompleta, por lo que se ha completado con las propiedades, disponibles en mayor detalle, de modelos conceptuales empleados en otros proyectos. En particular, las similitudes entre la turbina base escogida y el modelo empleado en el proyecto DOWEC hacen que este resulte la mejor opción a la hora de completar las especificaciones de la turbina NREL de 5 MW. El aerogenerador NREL de 5 MW ha sido empleado para establecer las especificaciones de referencia en varios proyectos de investigación aprobados por el Programa de Tecnologías Eólica e Hidráulica del Departamento de Energía de EE.UU.. Además, tanto el programa de investigación UpWind de la U.E. como la Agencia Internacional de Energía (IEA), han adoptado este modelo como referencia. En definitiva, es frecuente la elección como referencia de este modelo por parte de Inés Mateos Canals 29 CAPÍTULO 4 4.1. AEROGENERADOR NREL DE 5 MW equipos de investigación a escala mundial para estandarizar las especificaciones de turbinas marítimas y cuantificar los beneficios de la tecnología eólica, ya sea en tierra o alta mar. A continuación se detallan las características del aerogenerador de alta mar de NREL de 5 MW [25]: Potencia Rotor: Orientación, configuración Control Tren motriz Diámetro: rotor, buje Altura del buje Velocidad del viento: arranque, nominal, parada Velocidad del rotor: arranque, nominal Velocidad nominal en punta de pala Voladizo, inclinación del eje, precono Masa del rotor Masa de la góndola Masa de la torre Coordenadas del centro de masas 5 MW A barlovento, 3 palas velocidad variable, paso colectivo alta velocidad, multiplicador de varias etapas 126 m, 3 m 90 m 3 m/s, 11.4 m/s, 25 m/s 6.9 rpm, 12.1 rpm 80 m/s 5 m, 5◦ , 2.5◦ 110,000 kg 240,000 kg 347,460 kg (-0.2 m, 0.0 m, 64 m) Cuadro 4.1: Cacterísticas generales del aerogenerador NREL de 5 MW El aerogenerador REpower 5M cuenta con palas con curvatura previa incorporada, con el objetivo de aumentar el espacio libre sin necesidad de que el rotor se encuentre más saliente. Gran parte de las herramientas de simulación no tienen en cuenta la curvatura previa, por lo que el modelo presenta un precono de 2.5o , de acuerdo con un diseño de precono limitado y mayor curvatura previa en las palas. El diámetro del rotor indicado en la tabla no tiene en cuenta el efecto de precono de las palas, que reduce el diámetro total de la zona barrida por las palas. El diámetro exacto del rotor, asumiendo que las palas no se encuentran desviadas, es realmente (126m) cos(2, 5o ) = 125,88m y el área barrida por las palas es: π (125,88m)2 = 12, 445,3m2 (4.1) 4 4.1.1. Propiedades estructurales de las palas El aerogenerador de 5 MW de NREL cuenta con tres palas. La distribución de propiedades estructurales de las palas se basa en aquella de la pala LM Glasfiber (fibra de vidrio) de 62.6 m de longitud, que se emplea en el estudio DOWEC. Dado que las palas del estudio DOWEC eran 1.1 m más largas que las del modelo REpower de 5 MW, en el modelo NREL se acortan a 61.5 m. Las propiedades estructurales en la punta de la pala se obtienen por interpolación entre los valores conocidos para las palas de 61.2 m y 61.7 m. En la siguiente se detallan las principales características estructurales de las palas: 30 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Longitud Factor de escala de masa Masa (integrada) total Segundo momento de inercia (resp. raíz) Primer momento de inercia (resp. raíz) Localidación del centro de masas Ratio de amortiguamiento estructural 61.5 m 4.536 % 17,740 kg 11,776,047 km · m2 363,231 kg · m 20.475 m 0.477465 % Cuadro 4.2: Propiedades estructurales de las palas del aerogenerador NREL de 5 MW 4.1.2. Propiedades aerodinámicas de las palas De forma análoga al caso de las propiedades estructurales de las palas, las propiedades aerodinámicas del aerogenerador NREL de 5 MW se basan en las características de las palas del modelo empleado por el proyecto DOWEC. La tabla resume las propiedades aerodinámicas distribuidas de las palas, en los nodos, que se encuentran en el centro de cada elemento de la pala: Nodo (-) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 RNodos (m) 2.8667 5.6000 8.333 11.7500 15.8500 19.9500 24.0500 28.1500 32.2500 36.3500 40.4500 44.5500 48.6500 52.7500 56.1667 58.9000 61.6333 Torsión aerod. (◦ ) 12.308 13.308 13.308 13.308 11.480 10.162 9.011 7.795 6.544 5.361 4.188 3.125 2.319 1.526 0.863 0.370 0.126 DRNodos 2.7333 2.7333 2.7333 4.1000 4.1000 4.1000 4.1000 4.1000 4.1000 4.1000 4.1000 4.1000 4.1000 4.1000 2.7333 2.7333 2.7333 Cuerda (m) 3.542 3.854 4.167 4.557 4.652 4.458 4.249 4.007 3.748 3.502 3.256 3.010 2.764 2.518 2.313 2.086 1.419 Cuadro 4.3: Propiedades aerodinámicas de las palas del aerogenerador NREL de 5 MW Los nodos están numerados del 1 al 17, en la columna encabezada por RNodos se detalla su localización a lo largo del eje de la pala, desde el centro del rotor a las secciones transversales de las palas. Las longitudes de los elementos está recogida en la columna DRNodos y la suma de todas ellas resulta en la longitud total de la pala, 61.5 m, como indica la tabla 4.2. Por integración de la distribución de las cuerdas a lo largo de la pala, se obtiene una solidez del rotor de aproximadamente 5,16 %. Inés Mateos Canals 31 CAPÍTULO 4 4.1.3. 4.1. AEROGENERADOR NREL DE 5 MW Curva de potencia 6 Potencia (MW) 5 4 3 2 1 0 0 5 10 15 20 25 Velocidad (m/s) Figura 4.1: Potencia generada en función de la velocidad del aire [25] En el gráfico 4.1 se ha representado la potencia que es posible extraer en el rotor del aerogenerador según la velocidad de la corriente de aire incidente. 4.1.4. Parámetros: α, φ, CL y CD El ángulo de ataque tiene una influencia importante sobre el valor de los coeficientes de sustentación y resistencia, como se vio en capítulos anteriores. Según los gráficos que relacionan los coeficientes corregidos en función del ángulo de ataque el flujo incidente para perfiles aerodinámicos NACA64, la sustentación se maximiza para ángulos de ataque en torno a 10◦ . Para las simulaciones realizadas, se ha tomado este valor y se ha generalizado a la longitud íntegra de la pala. Los valores asociados de los coeficientes de sustentación y resistencia son CL = 1,5 y CD ≈ 0, respectivamente. Por otro lado, se ha considerado un ángulo de paso, γ, con valor 1◦ , dadas las características aerodinámicas del aerogenerador NREL. Las consideraciones anteriores conducen a un valor del ángulo φ de alrededor de 11◦ . 32 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Figura 4.2: Coeficientes corregidos del perfil NACA64 según el ángulo de ataque [25] 4.2. Simulación de un aerogenerador aislado Se ha llevado a cabo la simulación de un campo fluido de densidad ρ = 1,205kg/m3 , viscosidad µ = 1,808 · 10−5 P a · s y dimensiones 900x600 m con un mallado de ∆x = 3m y = 6m. En el punto de coordenadas (300,300) se sitúa el buje del aerogenerador NREL de 5 MW, cuyas palas barren una sección circular de 126 m de diámetro. Sobre el plano de giro del rotor incide una corriente de aire según la dirección del eje del rotor, a una velocidad de 8 m/s. 4.2.1. Condiciones de contorno y referencia La corriente de aire simulada se desplaza según la dirección del eje x, dirección horizontal en la vista escogida, de izquierda a derecha. Se establecen como nulas las componentes del vector velocidad según los ejes y (vertical) a la entrada del dominio. La componente vertical de la velocidad en los nodos varía, aunque en muy pequeña mediada, en el campo del fluido. Sin embargo, se impone una componente z nula en todos los nodos del dominio estudiado. Por otro lado, se fija la referencia de presión, p = 0,0 P a, en el extremo final (derecho) del campo fluido estudiado. Todos los valores de presión representados en las subsecuentes secciones comparten esta referencia. Estas condiciones de contorno y referencia se han empleado en todas las simulaciones llevadas a cabo en el desarrollo de este proyecto. 4.2.2. Velocidad En la figura 4.3 se muestra el campo de velocidades del fluido en una sección horizontal a la altura del eje del rotor. Inés Mateos Canals 33 CAPÍTULO 4 4.2. SIMULACIÓN DE UN AEROGENERADOR AISLADO Figura 4.3: Campo de velocidades, dominio 900x600 m, mallado ∆x = 3m La figura 4.4 recoge representaciones gráficas de la velocidad del aire en dos puntos de la estela tras el plano de giro de las palas del aerogenerador. Se ha superpuesto la distribución de velocidades obtenida por otro estudio [20] mediante una simulación de las mismas características: velocidad de incidencia 8 m/s y las mismas condiciones de contorno. En dicho proyecto se emplea otro método de resolución (volúmnes finitos) de una formulación alternativa de las ecuaciones de Navier-Stokes, se lleva a cabo la simulación en 3D y se tiene en cuenta la turbulencia y la variación radial del ángulo de ataque en las palas. 500 500 450 450 400 400 Posición y (m) Posición y (m) Cualitativamente, las gráficas son similares. Las disparidades entre ambas distribuciones se deben, a las diferencias entre los modelos empleados y, en menor medida, a la ligera diferencia entre mallados y valores del parámetro épsilon. El anterior estudio emplea un mallado de ∆x = 4,2m y = 8,4, frente a los valores ∆x = 3m y = 6 de la simulación realizada. 350 300 250 350 300 250 200 200 150 150 100 100 3 4 5 6 7 8 9 10 3 4 5 Velocidad (m/s) 6 7 8 9 10 Velocidad (m/s) (a) A un D del rotor (b) A 4D del rotor Figura 4.4: Velocidad en la estela 34 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos 4.2.3. Presión La figura 4.5 se trata de una representación análoga del campo de presiones en el mismo plano que en el caso de la velocidad. 600 600 500 500 400 400 Posición y (m) Posición y (m) Figura 4.5: Campo de presiones, dominio 900x600 m, mallado ∆x = 3m 300 300 200 200 100 100 0 -2 -1 0 1 2 3 0 -2 -1 Presión (Pa) 0 1 2 3 Presión (Pa) (a) A un D del rotor (b) A 4D del rotor Figura 4.6: Presión en la estela Las gráficas recogidas en la figura 4.6 muestran las variaciones frente a la presión de referencia establecida en cortes de la estela del aerogenerador. Se observa que las pendientes la distribución de presiones tienen un carácter mucho más acusado en el caso de la zona de la estela más próxima al aerogenerador. Según aumenta la distancia, la distribución de presiones se aplana, dado que en esta zona la influencia del aerogenerador es muy inferior. Inés Mateos Canals 35 CAPÍTULO 4 4.2.4. 4.2. SIMULACIÓN DE UN AEROGENERADOR AISLADO Potencia IRIS proporciona la solución de las ecuaciones de Navier-Stokes para flujo incompresible de forma que, tras la simulación, se conocen la presión y la velocidad vectorial en cada uno de los nodos del mallado. A partir de estos datos, es posible calcular la potencia extraída del viento por el aerogenerador. Para un elemento de pala, que puede denotarse como elemento i, la potencia sigue la siguiente expresión: Pi = ∆pvi Ai (4.2) Es decir, la potencia generada por la incidencia del viento sobre un elemento de pala es el producto de la variación de presión que tiene lugar entre la parte delantera y posterior del disco (respecto al sentido del flujo), por la velocidad del aire en ese punto y el área del elemento de pala en cuestión. El término de superficie de los elementos actuadores permite ponderar la contribución a la potencia total de cada uno. Dadas las consideraciones de simetría que se han tomado para la simulación del aerogenerador en IRIS, se ha realizado una media ponderada de los elementos actuadores superiores e inferiores al centro del rotor. Se cumple que: A= πD2 = ΣAi = Σπri ∆r 4 (4.3) donde ri es la distancia radial desde el centro del rotor hasta el centro del elemento actuador y ∆r es la distancia entre elementos consecutivos, que coincide con la magnitud de su dimensión radial. La potencia generada en la totalidad del disco es el resultado de la integración de los incrementos de potencia de cada elemento de pala: P = ΣPi (4.4) La magnitud de la potencia estimada mediante la simulación varía de acuerdo a las consideraciones que se tomen en cuanto al mallado y el término de la función gaussiana. A continuación se analiza la casuística. Influencia de la posición de los nodos Como se vio en el análisis aerodinámico del aerogenerador, la obstaculización del flujo de una corriente de aire da lugar una sobrepresión delante del obstáculo y una depresión detrás, en la parte próxima de la estela. La diferencia entre los valores de presión delante y detrás del rotor es una de las variables que intervienen en el cálculo de la potencia, de ahí que la selección de los nodos que se emplean en el cálculo de la misma tenga especial relevancia. En la representación del campo de presiones se aprecian ambos casos, la sobrepresión en tonos rojos y la depresión en tonos azules. Se observa asímismo una franja vertical de presiones cercanas a la presión de referencia (en color verde) entre las zonas de sobrepresión y depresión. La posición de esta franja coincide con la del plano de giro de las palas del rotor y es consecuencia del modelo empleado. La variación de presión se calcula entre los elementos del mallado que se encuentran delante y detrás del aerogenerador. Sin embargo, se comprueba gráficamente que las presiones de los elementos directamente adyacentes al plano de giro no toman valores esperados, pues se encuentran en la franja verde a la presión de referencia, de forma que no pueden ser empleados en el cálculo de la 36 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Figura 4.7: Campo de presiones en la zona próxima al rotor con mallado ∆x = 3m potencia. Para evitar el elevado margen de error al que conduciría el considerar las primeras columnas de nodos en el cálculo de la potencia estimada, se prescinde de los elementos que se encuentran en las proximidades del plano de giro de las palas del rotor, puesto que en ellos las presiones no alcanzan los valores esperados. Es fundamental determinar la distancia óptima entre el plano del rotor y las primeras columnas de elementos del mallado aguas arriba y aguas abajo que proporcionan una diferencia de presiones adecuada, es decir, una variación de presión que permita obtener resultados adecuados de potencia generada. Para conseguir este objetivo, se han lanzado simulaciones análogas con cuatro mallados diferentes en las que toma siempre un valor doble al de ∆x. La velocidad de la corriente de aire incidente es 8 m/s en todos los casos, por lo que la potencia esperada está en torno a 1.9 MW. En la figura 4.8 se ha representado gráficamente la evolución de la potencia que predice la simulación según el conjunto de nodos empleados para calcular la variación de presión para cuatro posibles mallados. El parámetro m designa la posición de una columna de nodos relativa al plano de giro del rotor. Como se ha mencionado, por razones de estabilidad, se prescinde en todos los casos de la columna de elementos adyacentes al plano del rotor, por lo que el valor mínimo que toma m es 2. Se observa que, en todos los casos, la potencia alcanza valores máximos para cierto valor de m. La tabla 4.4 que se presenta a continuación refleja los resultados: Inés Mateos Canals 37 CAPÍTULO 4 4.2. SIMULACIÓN DE UN AEROGENERADOR AISLADO 2.5 ∆x = 4 m ∆x = 3 m ∆x = 5 m ∆x = 6 m Potencia (MW) 2 1.5 1 0.5 2 3 4 5 6 7 m Figura 4.8: Potencia según la columna de nodos considerada (m) para mallas de diferente ∆x m ∆x = 3 ∆x = 4 ∆x = 5 ∆x = 6 2 0,9176 1,0062 0,983 0,9169 3 1,5032 1,4408 1,2938 1,1409 Potencia 4 1,8774 1,632 1,3851 1,1783 5 1,9893 1,6382 1,3519 1,1304 6 1,9655 1,5797 1,2876 1,0681 Cuadro 4.4: Potencia calculada según los nodos considerados Se han resaltado en negrita los valores de potencia máxima calculada, dadas las condiciones de la simulación. Para cada mallado, la posición en el eje x de los nodos empleados es: m ∆x = 3 ∆x = 4 ∆x = 5 ∆x = 6 2 6 8 10 12 Distancia 3 4 9 12 12 16 15 20 18 24 (m) 5 15 20 25 30 6 18 24 30 36 Cuadro 4.5: Distancia de los nodos considerados al plano de giro de las palas De nuevo, se ha resaltado en negrita la distancia en metros al aerogenerador de los nodos que proporcionan los valores de la diferencia de presiones más adecuados y que mejor permiten estimar la potencia. Los cuadros 4.4 y 4.5 permiten concluir que los elementos que proporcionan estimaciones de potencia más aproximadas al valor esperado se encuentran a una distancia superior a 12 metros del plano de giro de las palas del rotor. Esta distancia supone en torno a un 10 % de la longitud característica del problema, que es el diámetro del aerogenerador (126 m), lo que permite concluir 38 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos que se trata de una consideración válida del modelo. En adelante las potencias calculadas se obtienen considerando variaciones de presión promediadas de dos columnas de nodos consecutivas a distancias superiores al 10 % del diámetros en la dirección de incidencia del viento. Influencia del mallado La simulación con un mallado más fino proporciona valores de potencia más próximos a la realidad. 2.5 Potencia (MW) 2 1.5 1 0.5 0 3 3.5 4 4.5 5 5.5 6 ∆x (m) Figura 4.9: Potencia según el mallado (∆x, con = 2∆x) En la figura 4.9 se han representado los resultados de cuatro simulaciones de un aerogenerador con mallados en los que ∆x toma valores 3, 4, 5 y 6 m. La corriente de aire incide sobre el rotor a 8 m/s, velocidad a la que, según las especificaciones técnicas de la turbina, le corresponde una potencia generada en torno a 1.9 MW. En la gráfica se observa que la mejor aproximación a este valor es proporcionada por el mallado de ∆x = 3m y que mallados más gruesos dan lugar a valores que se alejan cada vez más de lo esperado según el tamaño del elemento aumenta. Influencia de los elementos de pala Los puntos sobre los cuales actúan las fuerzas elementales en la superficie actuadora están dispuestos de manera equidistante a lo largo del diámetro del disco. Entre ellos existe una distancia, que se ha denominado ∆b, de 0.3 m. Dado que el rotor tiene un diámetro de 126 m, la elección de ∆b resulta en un conjunto de 419 elementos actuadores. En las simulaciones realizadas, se ha mantenido el número de elementos actuadores constante con independencia del mallado empleado por IRIS para resolver el problema. Una alternativa a esto consiste en considerar la distancia entre elementos actuadores, ∆b, como función de la distancia entre nodos del mallado elegido, ∆x. Sin embargo, según los resultados de análisis [20] de la variación de la potencia como función de la relación entre ∆b y ∆x, las predicciones logradas convergen para valores del ratio ∆b/∆x inferiores a 0.5. En el caso de las simulaciones realizadas en este trabajo, las limitaciones de memoria no han permitido el uso de un mallado de mayor finura que ∆x = 3m, Inés Mateos Canals 39 CAPÍTULO 4 4.2. SIMULACIÓN DE UN AEROGENERADOR AISLADO por lo que no ha existido la posiblidad de superar el umbral del ratio ∆b/∆x a partir del cual la influencia de ∆b es relevante. Influencia de epsilon, Se han llevado a cabo una serie de simulaciones de un dominio con un único aerogenerador con la idea de comprobar las consecuencias que tiene para el método empleado el uso de valores de poco apropiados. La simulación de un campo fluido de 900x600 m con un mallado de ∆x = 3m y = 2m no da resultados convergentes. Queda patente que el uso de valores de inferiores al mallado dan lugar a inestabilidades numéricas. La misma simulación, con un valor de = 3m sí que permite alcanzar resultados convergentes, aunque son necesarias más iteraciones y un mayor esfuerzo de computación. Las figuras 4.10 y 4.11 muestran los campos de velocidad y presión. Se observa que la simulación no es correcta y existen inestabilidades. Figura 4.10: Campo de velocidades 900x600 m, mallado ∆x = 3m y = 3m 40 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Figura 4.11: Campo de presiones 900x600 m, mallado ∆x = 3m y = 3m Los resultados anteriores, así como análisis detallados de la influencia de sobre la resolución de las ecuaciones que se han llevado a cabo en otros estudios [20] llevan a seleccionar un valor óptimo de . En las simulaciones sucesivas, se emplea un valor de = 2∆x. Dicho valor proporciona resultados apropiados y no compromete la estabilidad de la solución 4.3. 4.3.1. Simulaciones de aerogeneradores agrupados Dos aerogeneradores en serie La simulación de dos aerogeneradores alineados permite llevar a cabo un análisis de la interacción de las estelas y la influencia que la modificación de las distribuciones de presiones y velocidades producida tiene sobre la potencia generada. En particular, la distribución simulada se trata de un caso analizado por otro estudio, cuyos resultados permiten analizar la validez del modelo como representación de la influencia de la posición relativa de aerogeneradores agrupados. Se pretende cuantificar el efecto que tiene situar un aerogenerador alineado con otro, de forma que la velocidad del flujo incidente sobre el aerogenerador que se encuentra en segundo plano no sea aquella de la corriente no perturbada en el infinito, por encontrarse el segundo rotor en la estela causada por el primero. Para lograr este objetivo, se ha llevado a cabo la simulación de un campo fluido de dimensiones 1500x600 m en cuyo seno se sitúan dos aerogeneradores NREL de 5 MW y rotor de 126 m. Los planos de giro de las palas están separados 630 m entre sí, distancia equivalente a cinco diámetros (5D). La velocidad del flujo incidente sobre el primer aerogenerador es 11.4 m/s. El mallado considerado es de ∆x = 3m y = 6m, como en el caso de un aerogenerador. Inés Mateos Canals 41 CAPÍTULO 4 4.3. SIMULACIONES DE AEROGENERADORES AGRUPADOS Velocidad Figura 4.12: Campo de velocidades, dominio 1500x600 m, mallado ∆x = 3m La figura 4.12 muestra el campo de velocidades que resulta de la simulación. Se observa que la estela provocada por el primer aerogenerador se prolonga una longitud superior a la distancia entre aerogeneradores, de forma que la velocidad de incidencia del aire sobre el segundo aerogenerador es muy inferior a aquella sobre el primero. 500 500 450 450 400 400 Posición y (m) Posición y (m) En las figuras 4.13 y 4.14 se presentan las distribuciones de velocidad en el campo del fluido correspondientes a las estelas del primer y segundo aerogenerador, respectivamente. Se trata de cortes según el eje y a medio diámetro (0,5D = 63 m) y cuatro diámetros (4D = 504 m) aguas abajo del aerogenerador correspondiente, respecto al plano de giro de las palas. 350 300 250 350 300 250 200 200 150 150 100 100 4 6 8 10 12 Velocidad (m/s) 14 4 6 8 10 12 14 Velocidad (m/s) (a) 0.5D (b) 4D Figura 4.13: Primer aerogenerador (A1). 42 Escuela Técnica Superior de Ingenieros Industriales (UPM) 500 500 450 450 400 400 Posición y (m) Posición y (m) Modelización por elementos finitos de parques eólicos 350 300 250 350 300 250 200 200 150 150 100 100 4 6 8 10 Velocidad (m/s) 12 14 4 6 8 10 12 14 Velocidad (m/s) (a) 0.5D (b) 4D Figura 4.14: Segundo aerogenerador (A2). En el caso del primer aerogenerador, la distribución de velocidades a lo largo del diámetro es cualitativamente muy similar a aquella de un simulador aislado. Sin embargo, la distribución de velocidades en la estela del segundo aerogenerador tiene un aspecto significativamente diferente a la representación equivalente de un aerogenerador aislado. Tanto la mayor amplitud de las variaciones de la velocidad en la sección estudiada como la forma de la distribución reflejan la interacción de las estelas de los aerogeneradores que tiene lugar aguas abajo del segundo aerogenerador. Presión Figura 4.15: Campo de presiones, dominio 1500x600 m, mallado ∆x = 3m El campo de presiones del fluido se ha representado en la figura 4.15. Es posible observar que, mientras que en el entorno del primer aerogenerador se producen una fuerte sobrepresión y depresión, delante y detrás de él, respectivamente, esto no ocurre en el caso del segundo aerogenerador. La presión en la estela del primer aerogenerador no alcanza valores bajos respecto a la referencia en las cercanías del segundo aerogenerador. Sin embargo, el aire incide sobre el segundo aerogenerador con una velocidad, que, aunque varía radialmente sobre la superficie del disco, es netamente inferior a la que incide sobre el primero, por lo que la sobrepresión delante del segundo disco es muy inferior. Inés Mateos Canals 43 4.3. SIMULACIONES DE AEROGENERADORES AGRUPADOS 600 600 500 500 400 400 Posición y (m) Posición y (m) CAPÍTULO 4 300 300 200 200 100 100 0 0 0 5 10 15 20 0 5 Presión (Pa) 10 15 20 15 20 Presión (Pa) (a) 0.5D (b) 4D 600 600 500 500 400 400 Posición y (m) Posición y (m) Figura 4.16: Primer aerogenerador (A1). 300 300 200 200 100 100 0 0 0 5 10 15 20 0 5 Presión (Pa) 10 Presión (Pa) (a) 0.5D (b) 4D Figura 4.17: Segundo aerogenerador (A2). En las figuras 4.16 y 4.17 se presentan las distribuciones de presión en el campo del fluido correspondientes a las estelas del primer y segundo aerogenerador, respectivamente. Se trata de cortes verticales a distancias de la mitad y cuatro diámetros del rotor respecto al plano de giro del mismo (0.5D y 4D). Es posible observar diferencias claras entre las gráficas anteriores, que cuentan con los mismos ejes para facilitar su interpretación. La primera (figura 4.16, (a)) correspondiente al corte próximo del aerogenerador situado en primer plano, representa una distribución de velocidades muy similar a aquella de un aerogenerador aislado. En el caso del corte equivalente de la estela del segundo aerogenerador (figura 4.17, (a)) el resultado es cualitativamente similar, aunque la variación de presiones en la estela es menos pronunciada. Esto es de esperar, teniendo en cuenta la situación del segundo aerogenerador, en la zona de influencia de la estela del primero. 44 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos En el caso del corte a 4D del rotor del primer aerogenerador se observan resultados que reflejan la influencia del segundo aerogenerador sobre el desarrollo normal (similar al de un aerogenerador aislado) de la estela. La sobrepresión provocada por la presencia un obstáculo en el flujo compensa totalmente la depresión que cabría esperar en la gráfica de una estela lejana, es más, existe cierta sobrepresión ya en esta zona, que se encuentra tan solo a la distancia de un diámetro del segundo aerogenerador. Por otro lado, la distribución de presiones en el corte lejano de la estela del segundo aerogenerador (figura 4.17, (b)) presenta desviaciones muy pequeñas frente a la presión de referencia. Potencia La magnitud de la velocidad en el plano del rotor y la variación de presiones que experimenta el fluido al atravesarlo presentan diferencias significativas que influyen de manera importante sobre la potencia generada por cada uno de los aerogeneradores. Según la curva de potencia del aerogenerador NREL de 5 MW, ante una corriente incidente de 11.4 m/s, el rotor del aerogenerador permite extraer una potencia que se encuentra en torno a los 5 MW. En el caso del primer aerogenerador, se espera un resultado de la simulación cercano a este valor. Sin embargo, dada la localización del segundo aerogenerador, la corriente de aire que incide sobre él presenta una velocidad y una presión diferentes a las del viento inalterado, al tratarse de la estela del primera aerogenerador, por lo que la potencia esperada es inferior a la primera. La estimaciones de potencia que resultan de la simulación son las siguientes: Aerogenerador A1 A2 Potencia 4.92 MW 1.62 MW Cuadro 4.6: Potencias de aerogeneradores alineados, distanciados 5D El resultado de la simulación permite concluir que existe una influencia considerable de la estela del primer aerogenerador sobre el funcionamiento del segundo. Las zonas de turbulencia y la magnitud de las mismas se ven incrementadas aguas abajo de la segunda turbina. La aparición de turbulencia tiene efectos significativos sobre el rendimiento y se une a la reducción de la velocidad incidente y de la presión aguas arriba del generador para resultar en una potencia generada muy inferior a aquella del primer aerogenerador. Inés Mateos Canals 45 Capítulo 5 Simulación de un parque eólico Lillgrund es un parque eólico offshore que se encuentra en una zona de poca profundidad (entre 4 y 8 m) de Öresund, a 7 km de la costa de Suecia y a otros 7 km del puente de Öresund, que conecta Suecia y Dinamarca. La empresa Vattenfall Vindkraft AB es propietaria y operadora del parque, que es un proyecto piloto apoyado por la Agencia de Energía de Suecia (STEM). El proceso de licitación finalizó en 2005 y el parque se construyó entre 2006 y 2007, desde Diciembre de 2007 se encuentra operacional. La velocidad media del viento en la zona toma valores medios de 8.5 m/s a la altura del eje del rotor. El parque integra 48 aerogeneradores iguales, con una potencia nominal de 2.3 MW, de forma que la potencia total de la planta alcanza los 110 MW [26]. Figura 5.1: Disposición de los aerogeneradores en Lillgrund [27] 5.1. Aerogenerador Siemens SWT-2.3-93 En el cuadro 5.1 quedan recogidas las características del aerogenerador de Siemens de potencia nominal 2.3 MW que opera en el parque eólico de Lillgrund. Inés Mateos Canals 47 CAPÍTULO 5 5.1. AEROGENERADOR SIEMENS SWT-2.3-93 Rotor Diámetro Área barrida Velocidad del rotor Regulación potencia 93 m 6,800m2 6-16 rpm Regulación de paso con velocidad variable Palas Tipo B45 Longitud 45 m Freno aerodinámico Tipo Paso de extensión completa Activación Activo, hidráulico Sistema de transmisión Tipo de multiplicador Planetario/helicoidal de 3 etapas Relación del multiplicador 1:91 Freno mecánico Tipo Freno de disco hidráulico Generador Tipo Asíncrono Potencia nominal 2.300 kW Tensión 690 V Sistema de refrigeración Intercambiador de calor integrado Sistema de orientación Tipo Activo Sistema de control Sistema SCADA WebWPS Control remoto Control pleno de la turbina Torre Tipo Tubular cilíndrico y/o cónico Altura del cubo 80 m o específico del emplazamiento Datos operativos Velocidad del viento de conexión 4 m/s Velocidad del viento potencia nominal 13-14 m/s Velocidad del viento de desconexión 55 m/s Pesos Rotor 60 toneladas Góndola 82 toneladas Torre según el emplazamiento Cuadro 5.1: Especificaciones técnicas del aerogenerador Siemens SWT-2.3-93 [28]. 5.1.1. Propiedades aerodinámicas de las palas El perfil aerodinámico del aerogenerador Siemens de Lillgrund pertenece al conjunto NACA 63. Aunque no se conoce con exactitud el perfil concreto dentro del conjunto anterior, pues no está especificado en la literatura técnica sobre el parque eólico, se considera que la opción más probable es el perfil NACA 63-415. Se trata de un perfil frecuentemente empleado en aerogeneradores modernos, que tiene las características de sustentación deseadas y ha sido empleado en otros estudios de 48 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos simulación de Lillgrund [29]. De este mismo estudio se han obtenido las longitudes de las cuerdas a lo largo de las palas del aerogenerador, que quedan recogidas en el diagrama 5.2. c (m) 6 4 2 0 0 5 10 15 20 25 30 35 40 45 50 r (m) Figura 5.2: Cuerda (c) en función de la distancia al centro del rotor En cuanto a la relación de velocidad de punta (λ), el coeficiente de potencia de la turbina NREL de 5 MW alcanza el valor máximo cuando este parámetro toma un valor de λ = 6,2329. Este es el valor que ha sido empleado en las simulaciones. 5.1.2. Curva de potencia 2.5 Potencia (MW) 2 1.5 1 0.5 0 0 5 10 15 20 25 Velocidad (m/s) Figura 5.3: Potencia generada en función de la velocidad del aire [26] En el gráfico 5.3 se ha representado la potencia que es posible obtener en el rotor del aerogenerador según la velocidad de la corriente de aire incidente. Inés Mateos Canals 49 CAPÍTULO 5 5.1.3. 5.1. AEROGENERADOR SIEMENS SWT-2.3-93 Simulación del aerogenerador SWT-2.3-93 aislado Como paso previo a la simulación del parque eólico al completo, se ha llevado a cabo una simulación del aerogenerador SWT-2.3-93 de Siemens aislado con mallado fino, de ∆x = 2m, con = 4m. El objetivo de la simulación es comprobar la validez del modelo desarrollado aplicado a un aerogenerador de especificaciones diferentes al empleado como base (NREL 5 MW). El dominio escogido coincide con el empleado en la simulación previa de un aerogenerador aislado, 900x600 m, con el centro del rotor situado en el punto de coordenadas (600,300). La velocidad de incidencia del viento es de 8 m/s, la misma que se empleó en la simulación del aerogenerador de 5 MW aislado. Velocidad En la figura 5.4 se muestra el campo de velocidades del fluido en una sección horizontal a la altura del eje del rotor. Figura 5.4: Campo de velocidades, dominio 900x600 m, mallado ∆x = 2m La figura 5.5 recoge representaciones gráficas de la velocidad del aire en dos puntos de la estela tras el plano de giro de las palas del aerogenerador. Existe una clara diferencia entre las distribuciones de velocidad en la estela producida por el aerogenerador NREL de 5 MW (figura 4.4) y aquella correspondiente al SWT-2.3-93. Dicha variación es consecuencia de los perfiles de las palas. Un El aerogenerador de diámetro de 126 m cuenta con palas NACA 64 y el de 93 m, con un perfil NACA 63-415, esto implica que las cuerdas correspondientes a los elementos de pala son diferentes. Como se ha visto, las cuerdas influyen directamente sobre la magnitud de la fuerza implementada por un elemento de pala y, como consecuencia, sobre la distribución de velocidades en la estela del aerogenerador. 50 Escuela Técnica Superior de Ingenieros Industriales (UPM) 500 500 400 400 Eje y (m) Eje y (m) Modelización por elementos finitos de parques eólicos 300 200 300 200 100 100 4 5 6 7 8 9 10 4 5 6 Velocidad (m/s) 7 8 9 10 Velocidad (m/s) (a) A un D del rotor (b) A 4D del rotor Figura 5.5: Velocidad en la estela Presión La figura 5.6 se trata de una representación del campo de presiones en el mismo plano que en el caso del campo de velocidades. Figura 5.6: Campo de presiones, dominio 900x600 m, mallado ∆x = 2m Potencia Según las especificaciones técnicas del aerogenerador, la potencia esperada para una velocidad de incidencia de 8 m/s es 0.9 MW. La potencia calculada mediante los resultados de la simulación toma un valor de 0.829 MW, según las consideraciones explicadas en capítulos anteriores. El resultado anterior permite concluir que se trata de una buena aproximación. Inés Mateos Canals 51 CAPÍTULO 5 5.2. 5.2. RESULTADOS DE LA SIMULACIÓN DE LILLGRUND Resultados de la simulación de Lillgrund La figura 5.1 mostraba la disposición Norte-Sur de los aerogeneradores, según esta orientación, los aerogeneradores están dispuestos diagonalmente, con los rotores orientados perpendicularmente al viento incidente. Para la simulación se emplea un campo fluido de dimensiones rectangulares y suficiente expansión como para contener en su dominio todos los aerogeneradores del parque eólico. Los límites laterales son paralelos a los planos de giro de las palas todos los rotores y perpendiculares a la dirección de la velocidad de la corriente que incide sobre ellos. En la literatura técnica del parque, cada uno de los aerogeneradores queda identificado mediante letras y números. Las filas se identifican por una letra, de la A (y = 200m) a la H (y = 2348,3m), empezando la asignación por la fila inferior. En lo que respecta a las columnas, la numeración avanza de derecha a izquierda. Se observan dos huecos en la distribución, que corresponderían a los aerogeneradores D05 y E05, que no pudieron ser instalados dada la imposibilidad de acceso de los buques de instalación por la insuficiente profundidad del agua en la zona. Para una lograr un mejor análisis de la variación de la velocidad que tiene lugar en la estela de los aerogeneradores, se ha representado la velocidad en secciones distanciadas un diámetro (1D = 93m) de cada aerogenerador. El siguiente código de colores ha sido empleado para representar a los aerogeneradores: 2500 Eje y (m) 2000 1500 1000 500 0 0 500 1000 1500 2000 2500 3000 Eje x (m) Figura 5.7: Disposición de los aerogeneradores en Lillgrund, MATLAB Los centros de los rotores de los 48 aerogeneradores que integran el conjunto están representados en la figura 5.7. Cada color identifica a un conjunto de aerogeneradores cuyos rotores se sitúan en una misma columna, es decir, en el mismo punto sobre el eje de abscisas. El mismo código de colores se ha empleado más adelante en el estudio de la distribución de velocidades en la zona próxima a los aerogeneradores de las estelas. Se ha simulado un dominio de 3800x3280 m para asegurar que el conjunto de aerogeneradores se encuentra alejado de los extremos, situación que podría conducir a inestabilidades. Debido a los restrictivos recursos de memoria de los equipos disponibles, la definición del mallado se ha visto 52 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos limitada. Los resultados más extensos, y los que quedan reflejados a continuación, corresponden a un mallado de ∆x = 8m, con un valor de = 16m. Estas limitaciones de definición tienen una gran influencia sobre los resultados de potencia estimada, que se alejan de los valores esperados. En lo que respecta a las condiciones de contorno y referencia, se han empleado los mismos valores que en el resto de simulaciones realizadas. 5.2.1. Velocidad En la figura 5.8 se ha representado el campo de velocidades que resulta de la simulación de una corriente incidente a 9 m/s. Es posible observar los efectos que tiene sobre la distribución de velocidades la interacción entre las estelas de los aerogeneradores. El efecto de la presencia de los aerogeneradores sobre la corriente es acumulativo, alcanzándose las mayores desviaciones en la zona final del dominio simulado. La simulación del mismo dominio con un mallado más fino ofrecería una mejor representación de la deceleración que tiene lugar en el entorno de los rotores. Figura 5.8: Campo de velocidades 3800x4160 m, mallado ∆x = 8m Las figuras 5.9 y 5.10, son representaciones gráficas de la distribución de velocidades en la estela de los aerogeneradores, a una distancia de un diámetro (93m) del plano de giro de las palas. Se han superpuesto las gráficas asociadas a las estelas de cada uno de los aerogeneradores que integran una fila a lo largo del eje x en el parque eólico. Los colores permiten diferenciar la posición en el eje x de los rotores (x0 ), según se indicó en la codificación reflejada en la figura 5.7. La misma codificación se emplea en todas las figuras. Inés Mateos Canals 53 CAPÍTULO 5 5.2. RESULTADOS DE LA SIMULACIÓN DE LILLGRUND 400 700 x0 =200 m 350 650 x0 =599.9 m x0 =999.8 m 300 600 250 x0 =1799.6 m Eje y (m) Eje y (m) x0 =1399.7 m x0 =2199.5 m 200 x0 =2599.4 m x0 =2999.3 m 150 550 500 450 100 400 50 350 300 0 3 4 5 6 7 8 9 10 3 4 5 6 7 Velocidad (m/s) Velocidad (m/s) (a) Fila A (b) Fila B 8 9 10 Figura 5.9: Distribución de velocidades (I). Azul asociado a la estela del aerogenerador en x0 = 200m, verde a x0 = 599,9m, negro a x0 = 999,8m, rojo a x0 = 1399,7m, amarillo a x0 = 1799,6m, magenta a x0 = 2199,5m, cian a x0 = 2599,4m, violeta a x0 = 2999,3m. La figura 5.9 corresponde a las filas inferiores de la distribución, en las que las estelas de los aerogeneradores presentan cierta desviación respecto a los ejes de los rotores. A pesar de encontrarse los aerogeneradores alineados, los puntos de velocidad mínima, generalmente asociados al buje, se hayan desplazados. Las distribuciones de velocidad representadas en las subfiguras (a), (b), (c) y (d) de la figura 5.10 muestran menor evidencia de la interacción entre las estelas de aerogeneradores situados en paralelo frente al flujo de aire incidente. Las zonas de menor velocidad se encuentran alineadas con el centro del rotor. Cabe destacar la presencia en las gráficas de líneas sin apenas desviación frente a la velocidad en el fluido circundante. Dichas curvas se asocian a huecos vacíos en la distribución de aerogeneradores y aparecen dado que se ha representado la velocidad en ciertas zonas del dominio, aunque no hubiese un aerogenerador situado un diámetro por delante de la sección estudiada en concreto. En el caso de los aerogeneradores situados en primer y segundo plano, la ausencia de perturbación en el fluido proporciona distribuciones de velocidades muy homogéneas. Sin embargo, en el caso particular de huecos en la distribución, las curvas de velocidad reflejan el desarrollo normal de la estela, sin que éste se vea modificado por la presencia de otro aerogenerador. En cuanto a la velocidad del aire en las estelas de las últimas filas de aerogneradores (filas G y H), que corresponden con las subfiguras (e) y (f) de la figura 5.10, se pueden observar de nuevo signos de desviación e interacción de las estelas en paralelo. Se trata de un efecto análogo al que tiene lugar en las filas inferiores del parque eólico. 54 Escuela Técnica Superior de Ingenieros Industriales (UPM) 1100 1400 1000 1300 900 1200 Eje y (m) Eje y (m) Modelización por elementos finitos de parques eólicos 800 700 1000 600 900 3 4 5 6 7 8 9 10 3 4 5 6 7 Velocidad (m/s) Velocidad (m/s) (a) Fila C (b) Fila D 1700 2000 1600 1900 1500 1800 Eje y (m) Eje y (m) 1100 1400 1300 8 9 10 8 9 10 8 9 10 1700 1600 1200 1500 3 4 5 6 7 8 9 10 3 4 5 6 7 Velocidad (m/s) Velocidad (m/s) (c) Fila E (d) Fila F 2550 2300 2500 2200 2100 Eje y (m) Eje y (m) 2450 2000 2400 2350 2300 2250 1900 2200 1800 2150 3 4 5 6 7 8 9 10 3 4 5 6 7 Velocidad (m/s) Velocidad (m/s) (e) Fila G (f) Fila H Figura 5.10: Distribución de velocidades (II). Azul asociado a la estela del aerogenerador en x0 = 200m, verde a x0 = 599,9m, negro a x0 = 999,8m, rojo a x0 = 1399,7m, amarillo a x0 = 1799,6m, magenta a x0 = 2199,5m, cian a x0 = 2599,4m, violeta a x0 = 2999,3m. Inés Mateos Canals 55 CAPÍTULO 5 5.2.2. 5.2. RESULTADOS DE LA SIMULACIÓN DE LILLGRUND Presión En la figura 5.11 se ha representado el campo de presiones en el dominio simulado. Es considerable la diferencia que existe entre la distribución de presiones en los entornos de aerogeneradores situados en primer plano y aquellos situados hacia el final del dominio. Sobre los tres primeros aerogeneradores incide el aire con una velocidad de 9 m/s, sin embargo, la presencia de aerogeneradores alineados causa la deceleración de la corriente, que alcanza los últimos aerogeneradores con velocidades muy inferiores a la del fluido no perturbado. La velocidad de incidencia tiene un fuerte efecto sobre la magnitud de la sobrepresión que se produce frente al rotor, causando velocidades bajas sobrepresiones muy inferiores. Por otro lado, como se mencionó en el caso de la simulación de dos aerogeneradores en serie, la potencia del aerogenerador en primer plano también se ve afectada por la presencia de los aerogeneradores que se encuentran aguas abajo. La caída de presión que tiene lugar en el rotor es inferior, dada la sobrepresión causada por el aerogenerador que se encuentra más adelante en el sentido del flujo. El salto de presiones (y con ello la potencia) en los aerogeneradores situados en primer plano es consecuentemente inferior al que proporcionaría el mismo aerogenerador en el seno de una corriente de similares características si se encontrara aislado. Figura 5.11: Campo de presiones 3800x4160 m, mallado ∆x = 8m 56 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos 5.2.3. Potencia La potencia especificada para una corriente de aire incidente a 9 m/s sobre un aerogenerador SWT-2.3-93 de Siemens está en torno a 1.2 MW según su curva de potencia. En ninguno de los aerogeneradores del parque eólico se alcanzaría dicha potencia ante un viento de 9 m/s, dadas las características del campo de presiones comentadas en el anterior apartado. La potencia extraída en los primeros aerogeneradores alcanzaría los valores más cercanos, aunque inferiores, a 1.2 MW. Debido al grosor del mallado empleado en la simulación, ∆x = 8m, y al elevado valor del parámetro épsilon asociado a éste, = 16m, la potencia estimada está muy lejos de ofrecer una buena aproximación. Sin embargo, los resultados de la simulación aportan información cualitativa sobre la interacción entre los aerogeneradores que integran el parque eólico. La potencia extraída por los rotores que no se encuentran en un primer plano respecto a la incidencia del viento no es fácilmente estimable mediante la curva de potencia, puesto que la velocidad con la que el aire choca con las palas del aerogenerador presenta una distribución no uniforme a lo largo del diámetro del rotor, al encontrarse éste en la estela de aerogeradores aguas arriba de él. En las figuras 5.12,5.13, 5.14 y 5.15 se ha representado la evolución de la potencia generada por los aerogeneradores que integran una fila del parque eólico. Como cabía esperar, la potencia decrece de forma importante en los aerogeneradores situados en la zona de influencia de aerogeneradores previos. 0.28 0.28 0.26 Potencia (MW) Potencia (MW) 0.26 0.24 0.22 0.24 0.22 0.2 0.18 0.2 0.16 0.18 500 0.14 1000 1500 2000 2500 3000 0 500 1000 1500 Eje x (m) Eje x (m) (a) Fila A (b) Fila B 2000 2500 3000 Figura 5.12: Distribución de velocidades (I) Inés Mateos Canals 57 5.2. RESULTADOS DE LA SIMULACIÓN DE LILLGRUND 0.3 0.3 0.25 0.25 Potencia (MW) Potencia (MW) CAPÍTULO 5 0.2 0.15 0.1 0.2 0.15 0.1 0.05 0.05 0 500 1000 1500 2000 2500 3000 0 500 1000 1500 Eje x (m) Eje x (m) (a) Fila C (b) Fila D 2000 2500 3000 Figura 5.13: Distribución de velocidades (II) Los resultados obtenidos describen cualitativamente bien la influencia que tiene sobre la potencia generada la alineación de los aerogeneradores en la dirección de la corriente de aire incidente. Esto se ha podido constatar por contraste con los resultados de otro estudio [20], que llevó a cabo una simulación de idénticas características, aunque haciendo uso de otro modelo. Sin embargo, los valores estimados de potencia generada distan mucho de ser buenas aproximaciones. Cabe esperar que, con un mallado fino (∆x = 2m), las predicciones en cuanto a potencia generada ofrecerían información realista sobre la capacidad del parque eólico ante unas determinadas condiciones del viento. 0.28 0.3 0.26 0.25 Potencia (MW) Potencia (MW) 0.24 0.2 0.15 0.22 0.2 0.18 0.16 0.1 0.14 0.05 500 1000 1500 2000 2500 3000 0.12 1000 1500 2000 Eje x (m) Eje x (m) (a) Fila E (b) Fila F 2500 3000 Figura 5.14: Distribución de velocidades (III) 58 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos 0.28 0.27 0.26 0.26 Potencia (MW) Potencia (MW) 0.25 0.24 0.22 0.2 0.24 0.23 0.22 0.21 0.18 0.16 1200 0.2 1400 1600 1800 2000 2200 2400 2600 0.19 1600 1800 2000 2200 Eje x (m) Eje x (m) (a) Fila G (b) Fila H 2400 2600 Figura 5.15: Distribución de velocidades (IV) Inés Mateos Canals 59 Valoración de impactos y de aspectos de responsabilidad legal, ética y profesional La energía eólica es una fuente de energía renovable, si bien es cierto que no carece de impacto ambiental, ya sea en términos de ruido, impacto visual o climático. La instalación de un parque eólico tiene consecuencias relevantes sobre el entorno [30]. Sin embargo, el presente trabajo tiene un carácter fundamentalmente teórico y la obtención de resultados se consigue mediante la realización de simulaciones, por lo que su impacto ambiental es mínimo. El consumo de recursos asociado al desarrollo de este proyecto se limita a la disposición de electricidad para el funcionamiento de los ordenadores. Dicho consumo eléctrico está dentro del margen razonable dada la magnitud del proyecto. En la ejecución de las tareas que integra este trabajo y la posterior redacción del presente documento se ha hecho uso de artículos y documentación técnica disponible al público, que se encuentra debidamente referenciada y recogida en el capítulo de Bibliografía. Por otro lado, el código empleado en las simulaciones en IRIS es resultado del trabajo del alumno y del personal del Departamento de Mecánica Estructural y Construcciones Industriales de la Escuela Técnica Superior de Ingenieros Industriales de Madrid. Lo mismo aplica al caso del código de MATLAB. En el desarrollo del proyecto y, en particular, de la modelización de un aerogenerador, se han llevado a cabo simplificaciones, que han quedado en todo caso reflejadas en el texto, justificadas y su influencia, analizada en lo posible. Las consideraciones anteriores se han abordado desde un punto de vista crítico y con la intención de alcanzar los objetivos marcados en el planteamiento del Trabajo de Fin de Grado (TFG). Inés Mateos Canals 61 Conclusiones y líneas futuras A continuación se expone un análisis del grado de alcance de los objetivos propuestos al inicio de este proyecto. En cuanto a la modelización de un aerogenerador, se ha realizado un estudio de modelos básicos de la aerodinámica de un aerogenerador con el objetivo de identificar un modelo adecuado que permita la futura optimización de un parque eólico. El modelo debe presentar un balance conveniente entre simplicidad, en términos de exigencia de poder computacional, y precisión, pues debe ser capaz de predecir con suficiente aproximación la potencia generada por un aerogenerador y, por extensión, por un parque eólico. El modelo del disco actuador (ADM) axisimétrico se fundamenta en la aplicación de fuerzas promediadas sobre la superficie barrida por las palas para facilitar la resolución de las ecuaciones de Navier-Stokes. El ADM logra una adecuada modelización del mecanismo de extracción de energía que tiene lugar en el rotor de un aerogenerador, si bien es cierto que no tiene en cuenta todas las características del flujo. Se trata de una opción viable que proporciona resultados de los campos de velocidad y presión en el dominio, así como de potencia extraída cercanos a los valores especificados, siempre en cuando la definición del mallado lo permita. La simulación de la interacción de una corriente de aire y el rotor de un aerogenerador se ha llevado a cabo mediante la implementación del modelo de ADM descrito, con ciertas modificaciones, en el programa de elementos finitos IRIS. Cabe destacar las siguientes consideraciones sobre el modelo empleado: Se trata de un modelo del buje y las palas del rotor, no considera el efecto de la torre sobre la carga aerodinámica en las palas. No se considera la topografía del terreno, el modelo sólo es válido en el caso de parques eólicos en el que los aerogeneradores están a la misma altura. El modelo no tiene en cuenta el efecto del suelo sobre la distribución de velocidades en el campo del fluido, ni las características asociadas a la capa límite atmosférica. La simetría axial del flujo incidente sobre el rotor permite resolver un problema bidimensional y, despreciando el término fuente debido a la geometría, emplear coordenadas cartesianas en vez de cilíndricas. Esta simplificación posibilita la simulación de un conjunto de aerogeneradores con IRIS. Debido a la consideración anterior, en la implementación del modelo no han sido contempladas las turbulencias y los efectos de punta de pala. La modelización de dichos fenómenos se lleva a cabo en tres dimensiones y exigen una gran capacidad de computación. Inés Mateos Canals 63 Las simulaciones han permitido realizar las siguientes observaciones sobre el modelo: La variación de presiones a ambos lados del disco (rotor) debe ser calculada haciendo uso de elementos del mallado alejados al menos un 10 % del diámetro del disco respecto al plano de giro de las palas, para evitar inestabilidades numéricas que lleven a la obtención de resultados incorrectos. Valores elevados del parámetro épsilon, , (relacionado con la distribución de las fuerzas sobre los elementos de pala) frente a las dimensiones del mallado, ∆x, proporcionan peores resultados en cuanto a potencia estimada. Por otro lado, la elección de valores de próximos a ∆x conduce a inestabilidades numéricas en las simulaciones. Se ha comprobado que la utilización de un valor de = 2∆x proporciona buenos resultados, por lo que se ha empleado en las simulaciones realizadas. El modelo refleja adecuadamente el comportamiento esperable del fluido ante una disposición en serie de aerogeneradores, reproduciendo la interacción entre las estelas y su efecto sobre la potencia generada. Finalmente, se ha simulado un parque eólico existente, el parque eólico marino de Lillgrund. Debido a las grandes dimensiones del dominio de simulación (3800x4160 m)y las limitaciones de memoria, el mallado de mejor definición empleado ha sido de ∆x = 8m. Dicho mallado aporta una aproximación cualitativa de las características del campo fluido, si bien es cierto que no proporciona estimaciones adecuadas de la potencia generada. Se ha llevado a cabo un análisis de la distribución de la velocidad en las estelas de los aerogeneradores y una comparación de las potencias estimadas para los aerogeneradores. En el desarrollo de este trabajo se ha pretendido allanar el camino para la futura utilización de IRIS en la optimización de un parque eólico. Este ha sido desde un principio el objetivo a largo plazo, que justifica la elección del modelo empleado y la simplificaciones que han sido efectuadas sobre él. El trabajo que se ha llevado a cabo abre la puerta a la creación de un código que permita optimizar la disposición de los aerogeneradores de un parque eólico con vistas a maximizar la potencia generada. La continuación de los esfuerzos iniciados por este proyecto tendría como requisito la ampliación de los recursos computacionales de los que dispone el departamento, en concreto, de la memoria. Existe la posibilidad de emplear un solver iterativo en vez de directo para la resolución del sistema lineal, que no exigiría tanta memoria. Esto permitiría la obtención de mejores aproximaciones con el modelo estudiado y la simulación de grandes dominios con mallados de mayor definición. Es más, la disposición de mayor capacidad de computación permitiría prescindir de ciertas simplificaciones del modelo y aumentar la precisión del mismo. En particular, no supone un gran salto conceptual la modelización en tres dimensiones con IRIS, en vez de en dos, como se ha hecho en las simulaciones de este trabajo. La simulación en tres dimensiones haría posible prescindir de simplificaciones importantes del modelo, como por ejemplo, las relacionadas con la variación del ángulo de ataque a lo largo de la pala. Es de esperar que la simulación en tres dimensiones mejoraría sustancialmente las predicciones cuantitativas del modelo. 64 Escuela Técnica Superior de Ingenieros Industriales (UPM) Planificación temporal y presupuesto Este proyecto se ha llevado a cabo durante el segundo cuatrimestre del curso académico 20152016 de forma simultánea con el octavo cuatrimestre del Grado en Ingeniería de las Tecnologías Industriales. Debido a esto, el desarrollo del proyecto no ha sido del todo lineal, sino que han tenido lugar pausas entre tareas, correspondientes a las convocatorias a exámenes, tanto ordinaria (enero y mayo/junio) como extraordinaria. A continuación se presenta la organización jerárquica y temporal de las diferentes tareas que comprende este proyecto. Estructura de descomposición del proyecto (EDP) Figura 5.16: EDP del proyecto Inés Mateos Canals 65 Diagrama de Gantt El siguiente diagrama de Gantt expone gráficamente el tiempo dedicado a cada una de las actividades del Trabajo de Fin de Grado. Figura 5.17: Diagrama de Gantt del proyecto 66 Escuela Técnica Superior de Ingenieros Industriales (UPM) Modelización por elementos finitos de parques eólicos Presupuesto En el desarrollo del presente TFG se han visto involucradas tres personas con distintos grados de implicación: el alumno y dos tutores, que han co-dirigido el proyecto. La siguiente tabla detalla el desglose de los costes que ha conllevado la realización del trabajo. Alumno Tutor 1 Tutor 2 Total Coste unitario (e/h) 12,5 45 45 Unidades (h) 320 80 48 Coste total (e) 4000 3600 2160 9760 Cuadro 5.2: Presupuesto del proyecto Notas: en la tabla anterior no se han incluido los costes del software empleado por las razones que se detallan a continuación: MATLAB: los alumnos de la UPM disponen de licencia académica gratuita. IRIS: se trata de un programa desarrollado por el propio departamento en el que se ha llevado a cabo el trabajo, su uso no conlleva desembolso. Paraview: software libre, gratuito. LaTex y Texstudio: software libre. Inés Mateos Canals 67 Bibliografía [1] EWEA: The European Wind Energy Association. Wind in power. 2015 european statistics. 2016. [2] Trevor J. Price. James blyth – britain’s first modern wind power engineer. Wind Engineering, 29:191–200, 2005. [3] Huang C. Watts, Jonathan. Winds of change blow through china as spending on renewable energy soars. The Guardian, 1992. [4] London array: the world’s largest offshore wind farm. Brochure. [5] Freude R.E. On the part played in propulsion by difference in pressure. Transactions of the Royal Institution of Naval Architects, 30:390–423, 1889. [6] Lanchester FW. A contribution to the theory of propulsion and the screw propeller. Transaction of the Institution of Naval Architects, 56:98–116. [7] Betz A. Das maximum der theoretisch möglichen ausnützung des windes durch windmotoren. Zeitschrift für das gesamte Turbinewessen, 56:307–309, 1920. [8] Glauert H. Airplane propellers. Aerodynamic Theroy, 4(Division L):169–269, 1963. Dover, New York. [9] Sørensen JN. Three-level, viscous-inviscid interaction technique for the prediction of separated flow past rotating wing. PhD Dissertation, (AFM 86-03,), 1986. [10] Berezin CR Torok MS. Berkman ME, Sankar LN. A navier-stokes/full potential/free wake method for rotor flows. AIAA Paper 97-0401, 1997. [11] Michelsen JA Sørensen NN. Hansen MOL, Sørensen JN. A global navier-stokes rotor prediction model. AIAA Paper 97-0970, 1997. [12] Schreck S. Sørensen NN, Michelsen JA. Navier-stokes predictions of the nrel phase-vi rotor in the nasa ames 80ft x 120 ft wind tunnel. Wind Energy, 5:151–169, 2002. [13] Michelsen JA Schreck S. Johansen J, Sørensen NN. Detached-eddy simulation of flow around the nrel phase vi blade. Wind Energy, 5:185–197, 2002. [14] Madsen HA. Sørensen JN, Gervang B. Vindmølle aerodynamik: Status og perspektiver. (Danish). Fluid Mechanics, (Report ET-AFM97-01). [15] A. Lecuona Neumann. La energía eólica: Princpios básicos y tecnología. Universidad Carlos III, Madrid. Inés Mateos Canals 69 [16] P. Jamieson. Innovation in Wind Turbine Design. Wiley, Chichester, United Kingdom, 2011. [17] E. Kulunk. Aerodynamics of wind turbines, fundamental and advanced topics in wind power. Technical report, 978-953-307-508-2, 2011. [18] Lu H. y Wu Y.T. Porté-Agel, F. A large-eddy simulation framework for wind energy applications. The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA, 2010. [19] Meneveau C. y Meyers J. Calaf, M. Large eddy simulation study of fully developed wind-turbine array boundary layers. Physics of Fluids, 22(1):015110, 2010. [20] S. Leonardi L.A. Martínez Tossas. Wind turbine modeling for computational fluid dynamics. Technical Report Subcontract Report NREL/SR-5000-55054, 2013. [21] Sørensen J.N y Mikkelsen R. Trodlborg, N. Numerical simulations of wake characteristics of a wind turbine uniform flow. Wind Energy, 13:86–99, 2010. [22] Hughes T.Jr. Brooks, A.N. Streamline-upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. Computer methods in applied mechanics and engineering, 1982. [23] T.E. Tezduyar. Finite element formulation for the voritcity-stream function form of the incompressible euler equations on multiply-connected domains. Computer methods in applied mechanics and engineering, 73:331–339, 1982. [24] Sørensen J. Nørkær Mikkelsen, R. Flemming. Actuator disk methods applied to wind turbines. 2004. [25] W. Musial y G. Scott J. Jonkman, S. Butterfield. Definition of a 5-mw reference wind turbine for offshore system development. Technical report, 2009. [26] Larsen P.E. Larsson Å. Jeppsson, J. Technical description of lillgrund wind power plant. Technical report, Vattenfall Vindkraft AB, 2008. [27] J.-Å. Dahlberg. Assessment of the lillgrund windfarm: Power performance wake effects. Technical report, Vattenfall Vindkraft AB, 2009. [28] Siemens Wind Power A/S. Technical specifications 2.3 mw mkii. Technical Report PG-R-0310-0000-0002-04, Siemens, 2005. [29] Wolf-Gerrit F Maguire A.E. Creech, A. Simulations of an offshore wind farm using large eddy simulation and torque-controlled actuator disc model (postprint). Institute of Energy Systems, University of Edimburgh, Scotland, 2014. [30] Yang-Y. Leung, D.Y.C. Wind energy development and its environmental impact: A review. Renewable and Sustainable Energy Reviews, 16:1031–1039, 2012. 70 Escuela Técnica Superior de Ingenieros Industriales (UPM) Índice de figuras 1. Evolución de la potencia eólica instalada en la Unión Europea. Fuente: EWEA [1] . vii 1.1. 1.2. 1.3. 1.4. 1.5. 1.6. Evolución de la velocidad y la presión del flujo [15]. . Relación de velocidades y areas en la corriente [16]. . Velocidad angular impuesta sobre una sección anular Geometría de un perfil aerodinámico subsónico [15] . Teoría del elemento de pala [16] . . . . . . . . . . . . Velocidades en un elemento de pala [15] . . . . . . . . . . . . . [17] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 4 5 7 10 11 3.1. Fuerzas [24] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Archivo de texto generado con MATLAB como input para IRIS . . . . . . . . . . . . 23 26 4.1. Potencia generada en función de la velocidad del aire [25] . . . . . . . . . . . 4.2. Coeficientes corregidos del perfil NACA64 según el ángulo de ataque [25] . . . 4.3. Campo de velocidades, dominio 900x600 m, mallado ∆x = 3m . . . . . . . . 4.4. Velocidad en la estela . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5. Campo de presiones, dominio 900x600 m, mallado ∆x = 3m . . . . . . . . . . 4.6. Presión en la estela . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7. Campo de presiones en la zona próxima al rotor con mallado ∆x = 3m . . . . 4.8. Potencia según la columna de nodos considerada (m) para mallas de diferente 4.9. Potencia según el mallado (∆x, con = 2∆x) . . . . . . . . . . . . . . . . . . 4.10. Campo de velocidades 900x600 m, mallado ∆x = 3m y = 3m . . . . . . . . 4.11. Campo de presiones 900x600 m, mallado ∆x = 3m y = 3m . . . . . . . . . 4.12. Campo de velocidades, dominio 1500x600 m, mallado ∆x = 3m . . . . . . . . 4.13. Primer aerogenerador (A1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14. Segundo aerogenerador (A2). . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15. Campo de presiones, dominio 1500x600 m, mallado ∆x = 3m . . . . . . . . . 4.16. Primer aerogenerador (A1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17. Segundo aerogenerador (A2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ∆x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 33 34 34 35 35 37 38 39 40 41 42 42 43 43 44 44 5.1. 5.2. 5.3. 5.4. 5.5. 5.6. 5.7. 5.8. . . . . . . . . . . . . . . . . 47 49 49 50 51 51 52 53 Disposición de los aerogeneradores en Lillgrund [27] . . . . . . Cuerda (c) en función de la distancia al centro del rotor . . . . Potencia generada en función de la velocidad del aire [26] . . . Campo de velocidades, dominio 900x600 m, mallado ∆x = 2m Velocidad en la estela . . . . . . . . . . . . . . . . . . . . . . . Campo de presiones, dominio 900x600 m, mallado ∆x = 2m . . Disposición de los aerogeneradores en Lillgrund, MATLAB . . Campo de velocidades 3800x4160 m, mallado ∆x = 8m . . . . Inés Mateos Canals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.9. Distribución de velocidades (I). Azul asociado a la estela del aerogenerador en x0 = 200m, verde a x0 = 599,9m, negro a x0 = 999,8m, rojo a x0 = 1399,7m, amarillo a x0 = 1799,6m, magenta a x0 = 2199,5m, cian a x0 = 2599,4m, violeta a x0 = 2999,3m. 5.10. Distribución de velocidades (II). Azul asociado a la estela del aerogenerador en x0 = 200m, verde a x0 = 599,9m, negro a x0 = 999,8m, rojo a x0 = 1399,7m, amarillo a x0 = 1799,6m, magenta a x0 = 2199,5m, cian a x0 = 2599,4m, violeta a x0 = 2999,3m. 5.11. Campo de presiones 3800x4160 m, mallado ∆x = 8m . . . . . . . . . . . . . . . . . . 5.12. Distribución de velocidades (I) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13. Distribución de velocidades (II) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14. Distribución de velocidades (III) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15. Distribución de velocidades (IV) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16. EDP del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17. Diagrama de Gantt del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 54 55 56 57 58 58 59 65 66 Escuela Técnica Superior de Ingenieros Industriales (UPM) Índice de cuadros 4.1. 4.2. 4.3. 4.4. 4.5. 4.6. Cacterísticas generales del aerogenerador NREL de 5 MW . . . . . . . . . . Propiedades estructurales de las palas del aerogenerador NREL de 5 MW . Propiedades aerodinámicas de las palas del aerogenerador NREL de 5 MW Potencia calculada según los nodos considerados . . . . . . . . . . . . . . . Distancia de los nodos considerados al plano de giro de las palas . . . . . . Potencias de aerogeneradores alineados, distanciados 5D . . . . . . . . . . . . . . . . . 30 31 31 38 38 45 5.1. Especificaciones técnicas del aerogenerador Siemens SWT-2.3-93 [28]. . . . . . . . . . 5.2. Presupuesto del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 67 Inés Mateos Canals . . . . . . . . . . . . . . . . . . . . . . . . 73 Abreviaturas, unidades y acrónimos HAWT: del inglés, “horizontal- axis wind turbine”, en español, aerogenerador de eje horizontal. Teoría BEM: del inglés, “blade element momentum”, teoría del momento sobre el elemento de pala. ADM: del inglés, “actuator disk model”, modelo del disco actuador. ALM: del inglés, “actuator line model”, modelo de la linea actuadora. ASM: del inglés, “actuator surface model”, modelo de la superficie actuadora. ABL: del ingles, “atmospheric boundary layer”, capa límite atmosférica. CFD: del inglés, “computational fluid dynamics”, mecánica de fluidos computacional. DNS: del inglés, “direct numerical simulation”, simulación directa. LES: del inglés, “large eddy simulation”, simulación de grandes torbellinos. N-S: Navier-Stokes. RANS: del inglés, “Reynolds-averaged Navier-Stokes”, ecuaciones de Navier-Stokes promediadas de Reynolds. SGS: del ingles, “subrgrid-scale”, escala subred. MEF: método de elementos finitos. NREL: del inglés, “National Renewable Energy Laboratory”, Laboratorio Nacional de Energía Renovable de EE.UU.. Inés Mateos Canals 75 Índice de símbolos Aerodinámica de aerogeneradores a a0 A AP AR c CL CD CN CT Cp Cpmax D D e f F K L ṁ m N p P Q r t t T u v vtheta vp vrel w — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — Inés Mateos Canals Coeficiente de inducción axial. Coeficiente de inducción tangencial. Área o sección (m2 ). Área del perfil (m2 ). Área barrida por las palas (m2 ). Cuerda (m). Coeficiente de sustentación (-). Coeficiente de resistencia (-). Coeficiente de fuerza normal (-). Coeficiente en el plano del disco (-). Coeficiente de potencia (-). Coeficiente de potencia máximo o Límite de Betz (-). Fuerza de resistencia aerodinámica (del inglés, Drag) (N). Diámetro del rotor (m). Vectores de dirección unitarios: ex , ey y ez (-). Fuerza diferencial. Fuerza (N). Constante característica de la fuerza elemental. Fuerza de sustentación aerodinámica (del inglés Lift) (N). Flujo, gasto másico (m3 /s). Identificador de columna de elementos (-). Número de palas (-). Presión (Pa). Potencia (W). Par (N·m). Variable radio (m). Espesor (m). Tiempo (s). Impulso (N). [ux uy uz ] Velocidad (m/s). Velocidad (m/s). Componente tangencial de la velocidad (m/s). Componente radial de la velocidad (m/s). Velocidad relativa (m/s). Anchura, del inglés width (m). 77 Letras griegas α λ λ π ρ Φ ϕ ∂ ω Ω Σ ∇ ν η θ γ φ ∆x — — — — — — — — — — — — — — — — — — — Ángulo de ataque (◦ ). Rapidez local de la pala (-). Relación de velocidad de punta (-). 3.1416 (-). Densidad (kg/m3 ). Diámetro (m). Ángulo de la corriente (◦ ). Operador derivada parcial (-). Velocidad angular del fluido (rad/s). Velocidad angular del rotor (rad/s). Solidez (-). Operador nabla. Viscosidad (P a · s). Parámetro épsilon (m). Función gaussiana. Ángulo de rotación (rad). Ángulo de paso (rad). Ángulo (rad). Dimensión del mallado (m). Subíndices 0 — 1 2 — — Condiciones aguas abajo del aerogenerador, es decir, de la corriente de aire previa a la perturbación (-). Condiciones en las cercanías del disco actuador (-). Condiciones de la corriente de aire aguas abajo del aerogenerador (-). Números adimensionales Re M 78 — — Número de Reynolds (-). Número de Mach (-). Escuela Técnica Superior de Ingenieros Industriales (UPM) Glosario Aerogenerador: generador eléctrico que funciona convirtiendo la energía cinética del viento en energía mecánica a través de una hélice y en energía eléctrica gracias a un alternador. Góndola o nacelle: alojamiento de los elementos mecánicos y eléctricos del aerogenerador. Rotor: parte de un aerogenerador (que incluye el buje y las palas) que gira y permite transformar la energía cinética del viento en par en el eje del generador eléctrico. Estela: rastro que deja tras de sí un cuerpo en movimiento en el seno de un fluido. Advección: desplazamiento horizontal, meridiano o zonal de una masa de aire, lo que provoca cambios de tiempo y transferencias de calor de unas zonas a otras de la superficie terrestre. Turbulencia: movimiento desordenado de un fluido en el cual las moléculas, en vez de seguir trayectorias paralelas, describen trayectorias sinuosas y forman torbellinos. Línea de corriente: aquella curva cuya tangente en cualquier punto coincide con la dirección de la velocidad del fluido en dicho punto. Capa límite: en mecánica de fluidos, zona donde el movimiento de un fluido es perturbado por la presencia de un sólido con el que está en contacto. Método numérico: conjunto de algoritmos que, a través de números y reglas matemáticas simples, simula procesos matemáticos más complejos. Modelización: uso de modelos físicos, matemáticos o representaciones lógicas de un sistema, entidad, fenómeno o proceso para llevar a cabo simulaciones y obtener datos con el objetivo de tomar decisiones técnicas. Axisimétrico: que presenta simetría axial, es decir, respecto a un eje. Inés Mateos Canals 79
© Copyright 2024