Source: http://docplayer.es/1427208-.html
Timestamp: 2017-08-23 16:12:30+00:00

Document:
Óscar Suárez Ávila
7 Capítulo 1 Introducción y planteamiento del proyecto Motivación del proyecto Objetivos del proyecto Metodología de trabajo Recursos a utilizar... 4 Capítulo 2 Introducción al efecto Magnus Enunciación del efecto Magnus Aplicación del flujo potencial al estudio del Efecto Magnus Aplicación de la teoría de la capa límite de Prandtl al estudio del Efecto Magnus Estudios experimentales sobre el comportamiento del Efecto Magnus Efecto Magnus a bajos números de Reynolds (Re<300) Efecto Magnus a altos números de Reynolds ( 10 ) Capítulo 3 Introducción a la dinámica de fluidos computacional Historia de la dinámica de fluidos computacional Métodos de discretización del problema Métodos de discretización de las diferencias finitas Métodos de discretización de los volúmenes finitos Métodos de discretización de los elementos finitos Mallado Tipos de malla: Elección de malla adecuada: Comprobación de la malla Modelos de turbulencia Reynolds-averaged Navier-Stokes (RANS): Reynolds Stress Model (RSM): Direct numerical simulation (DNS) Modelo LES (Large Eddy Simulation) Detached Eddy Simulation (DES)
8 Capítulo 4 Simulaciones propias del Efecto Magnus Simulación bidimensional del efecto Magnus Geometría: Mallado: Condiciones de contorno: Análisis de los resultados obtenidos: Simulación tridimensional del cilindro infinito Geometría: Mallado: Condiciones de contorno: Resultados: Simulación de cilindro con efecto de borde en 3D Geometría: Mallado: Condiciones de contorno: Resultados obtenidos: Capítulo 5 Optimización de la geometría Estudio de las variables que influyen en el rendimiento del aerogenerador Número óptimo de palas: Diámetro óptimo de las palas: Obtención de la expresión de la potencia producida por el aerogenerador Análisis de los resultados obtenidos para el diseño óptimo Valores resultantes de la optimización de ω y R Par en el eje Potencia Coeficiente de potencia Geometría de las palas Capítulo 6 Conclusiones Estudio del Efecto Magnus
9 6.2 Dinámica de fluidos computacional Simulaciones propias del Efecto Magnus Simulaciones bidimensionales Simulaciones tridimensionales Modelo numérico. Optimización Recomendaciones para futuros estudios Bibliografía Apéndice A Código Matlab
10 Capítulo 1 INTRODUCCIÓN Y PLANTEAMIENTO DEL PROYECTO.
11 Introducción y planteamiento del proyecto Motivación del proyecto. A lo largo de este proyecto se realizará un estudio profundo del efecto Magnus y en particular de su aplicación en la obtención de energía eólica. Para ello primero se hará acopio de la bibliografía disponible y partiendo ella se realizarán las simulaciones necesarias para poder estudiar el funcionamiento de este tipo de aerogeneradores y optimizar así las variables que influyen en su rendimiento. El efecto Magnus fue enunciado en 1852 por Heinrich Gustav Magnus y describe la fuerza de sustentación que experimenta un cilindro en rotación inmerso en una corriente de aire. Desde entonces ha despertado el interés de multitud de científicos. Este fenómeno influye fuertemente en la trayectoria de cualquier objeto en rotación desplazándose a través de un fluido, de ahí que los estudios más exhaustivos de este fenómeno estén relacionados con la balística, para evaluar las posibles desviaciones en trayectorias de proyectiles. En menor medida se ha investigado su aplicación al deporte de élite estudiando su influencia en el comportamiento de los cuerpos esféricos en movimiento utilizados. Debido a la expectación que despertó este efecto tras su descubrimiento, se realizaron investigaciones para aplicarlo a multitud de áreas. Aunque sin excesivo éxito, en los años veinte el alemán Anton Flettner intentó su aplicación en la propulsión naval inventando el rotor Flettner e implantándolo en un prototipo llamado Buckau o Baden-Baden. También se intentó utilizar para construir aviones con cilindros en vez de alas, como el prototipo de hidroavión modelo 921-V de 1930, fabricado en Estados Unidos, que fracasó nuevamente. De entre todas sus aplicaciones la más interesante para este estudio es la creación de aerogeneradores basados en dicho efecto. Han sido varios
12 Introducción y planteamiento del proyecto. 3 los proyectos que lo han intentado pero todos se han encontrado con problemas de ineficiencia inesperados. La empresa japonesa MECARO construyó en 2007 un prototipo de aerogenerador basado en el efecto Magnus con la particularidad de tener dispuesta una espiral enroscada a lo largo de los cilindros. Según sus estudios la superficie lisa de los cilindros exigía que giraran a muy altas revoluciones para tener una sustentación apreciable, lo que hacía consumir demasiada energía. Con las nuevas espirales parecen haber solucionado el problema y construido el primer aerogenerador basado en este efecto lo suficientemente eficiente como para comercializarse. Estos nuevos descubrimientos en la materia han llevado a la universidad a realizar proyectos de investigación sobre esta posible área de desarrollo para comprender mejor la tecnología y analizar entre otras cuestiones su viabilidad frente a los aerogeneradores actuales. 1.2 Objetivos del proyecto. El objetivo de este proyecto es continuar con el estudio de los aerogeneradores eólicos basados en el efecto Magnus iniciado en [GARC09]. Para ello se estudiará con el fin de optimizar la geometría y parámetros de funcionamiento del aerogenerador. Para ello se 1. Estudiar el efecto Magnus y en particular su aplicación a los aerogeneradores. 2. Analizar las curvas disponibles de y y su validez para emplearlas en el estudio de aerogeneradores. 3. Estudiar exhaustivamente técnicas dinámica de fluidos computacional (CFD), los diferentes modelos para la simulación de fluidos y los métodos de mallado.
13 Introducción y planteamiento del proyecto Realizar simulaciones de los perfiles de los cilindros para compararlas con las curvas disponibles y optimizar los parámetros del cilindro. 5. Comparación de los resultados obtenidos para el cilindro con los disponibles para otros perfiles alternativos. 6. Optimización de los parámetros del aerogenerador. 7. Con los resultados obtenidos se analizarán las posibles ventajas frente a los aerogeneradores actuales. 1.3 Metodología de trabajo Se estudiará el comportamiento de un aerogenerador basado en el efecto Magnus. Primeramente se realizará un estudio sobre el estado del arte del Efecto Magnus. Posteriormente se realizarán simulaciones mediante CFD (ANSYS) de un cilindro en rotación y con el fin de obtener las gráficas que describen el comportamiento de los coeficientes y. Se compararán las gráficas obtenidas con la bibliografía encontrada para contrastar su validez. Con los resultados obtenidos se estudiará se elaborará un modelo matemático del aerogenerador para estudiar la optimización de su geometría. Con las conclusiones que arroje el modelo se analizarán los parámetros que influyen en el rendimiento del aerogenerador con el fin de optimizarlos. 1.4 Recursos a utilizar Se recogerá toda la información posible sobre el efecto Magnus, curvas de C y C de cilindros en rotación y los estudios sobre la aplicación del Efecto Magnus en aerogeneradores. Para ello se consultará la bibliografía de Mecánica de Fluidos disponible así como de las publicaciones y patentes relacionadas con la materia. Posteriormente se
14 Introducción y planteamiento del proyecto. 5 utilizará la bibliografía disponible acerca de la mecánica de fluidos computacional para estudiar y evaluar el funcionamiento de esta herramienta. A lo largo del proyecto se utilizará el siguiente software: Microsoft Word Microsoft Excel ANSYS version (CFX Y FLUENT). Matlab versión (R2007b). GIMP
15 Capítulo 2 INTRODUCCIÓN AL EFECTO MAGNUS
16 Introducción al efecto Magnus Enunciación del efecto Magnus. El Efecto Magnus es el fenómeno físico por el cual un cilindro, esfera o cualquier cuerpo de revolución, en rotación sobre su eje de revolución, que se mueve a través de un fluido siguiendo una dirección no paralela a dicho eje, experimenta una fuerza aerodinámica de sustentación. Esta fuerza es perpendicular al plano formado por la dirección de traslación y el eje de rotación y su magnitud depende de varios factores, como por ejemplo la velocidad de desplazamiento del objeto o la velocidad de rotación, que serán estudiadas a lo largo de este capítulo. Esta fuerza y su momento respecto al centro de gravedad juegan un papel importante en la trayectoria de proyectiles o cualquier objeto en rotación que siga una trayectoria a través del aire, por lo que este efecto ha suscitado un gran interés en áreas como la aerodinámica o la balística. Fue observado por primera vez en 1672 por Sir Isaac Newton, quien en su escrito [NEWT72] hace referencia a las curiosas trayectorias que a menudo siguen las pelotas de tenis al imprimirlas ciertos golpes. En 1742 B. Robins, ingeniero de artillería británico, explicó ciertas desviaciones inesperadas en la trayectoria de proyectiles como una consecuencia del este efecto. En 1852 G. Magnus [MAGN52] realizó ensayos y mediciones experimentales con un cilindro girando sobre su eje de revolución inmerso en un flujo de aire perpendicular a dicho eje. Constató que, debido a la combinación del giro del cilindro y la velocidad del aire, se generaban diferentes presiones en la parte superior e inferior del cilindro. Para estudiar este fenómeno en esferas, utilizó balas de mosquete esféricas con su centro de masa desplazado de su centro geométrico. Al situar el centro de masa a la izquierda o derecha del eje del cañón del arma conseguía imprimir a la bala un spin hacia la derecha o la izquierda respectivamente. Observó que según la dirección de spin que imprimía a la bala ésta se
17 Introducción al efecto Magnus 8 desviaba hacia la izquierda o a la derecha, demostrando la influencia del spin en la trayectoria de proyectiles y enunciado el teorema que llevaría su nombre. 2.2 Aplicación del flujo potencial al estudio del Efecto Magnus. Tras los avances de G. Magnus, en 1877 Lord Rayleigh publicó un estudio sobre las trayectorias irregulares de las pelotas de tenis [RAYL77], en el que confirmaba la explicación del fenómeno propuesta por Magnus y enunciaba el primer modelo matemático del mismo. Utilizó la teoría del flujo potencial, válida exclusivamente para fluidos no viscosos, para obtener una solución del flujo ideal que interviene en este efecto. Demostró que al superponer una corriente uniforme (Figura 1) y un torbellino irrotacional (Figura 2), aparecía una fuerza perpendicular a dicha circulación equivalente a la producida por el Efecto Magnus. Esta aportación sustentación. sería la base del teorema de Kutta-Joukowsky de la
18 Introducción al efecto Magnus 9 Figura 1. Corriente uniforme en torno a un cilindro [PRAN25]. Figura 2. Torbellino irrotacional alrededor de un cilindro [PRAN25]. Lord Rayleigh confirmó la explicación del fenómeno dada por Magnus y advirtió que la fricción entre la superficie del cuerpo en rotación y el fluido jugaba un papel clave en la curvatura del flujo que causa la fuerza asociada al efecto Magnus, sin embargo su teoría sólo podía considerarse como una aproximación ideal y nunca como una explicación matemática de la física real que interviene en este proceso. Su modelo, aplicable exclusivamente a fluidos no viscosos, no tenía en cuenta la fricción y por tanto no podía establecer ningún mecanismo para establecer la circulación del fluido alrededor del cilindro ni predecir la fuerza externa que el cuerpo en rotación comunica al fluido. 2.3 Aplicación de la teoría de la capa límite de Prandtl al estudio del Efecto Magnus. El problema que suponía explicar la causa de la circulación del fluido alrededor del cuerpo en rotación fue resuelto en 1918 por L. Prandtl [PRAN25], con la aplicación de su teoría de la capa límite al estudio del Efecto Magnus. Prandtl explica que las fuerzas debidas a la fricción que se producen en el interior de los fluidos no viscosos son generalmente
19 Introducción al efecto Magnus 10 depreciables en comparación con las fuerzas de inercia, pero existe una fina región de fluido, situada junto a las paredes de cualquier cuerpo inmerso en dicho fluido, en que dichas fuerzas son del mismo orden que las fuerzas de inercia y por tanto deben ser tenidas en cuenta. Para afirmar esto, se basó en observaciones que se habían llevado a cabo en experimentos con fluidos viscosos. En ellos se apreciaba cómo la capa de fluido adyacente al cuerpo estaba relativamente en reposo. En las siguientes capas de fluido, a medida que éstas se alejaban de la pared, la velocidad aumentaba respecto a la capa anterior (Figura 3). Figura 3. Distribución de velocidades a lo largo de la capa límite Esto implica, por tanto, la existencia de una capa alrededor del cuerpo en la que se daba la transición en la velocidad desde el reposo hasta la velocidad del flujo que no se ve afectado por la viscosidad producida por las fuerzas de fricción. Esta teoría justifica la circulación en torno al cilindro y explica también el desprendimiento de la capa límite de la pared y la creación de vórtices. Prandtl explica que en un fluido no viscoso incidiendo perpendicularmente en un cilindro inmóvil (Figura 1), cada partícula seguiría una línea de corriente conservando su energía. En los fluidos viscosos, debido a la fricción, las partículas cercanas a la pared perderían parte de su energía no pudiendo así seguir su trayectoria ideal. En la
20 Introducción al efecto Magnus 11 Figura 1, una partícula de un fluido no viscoso que llega al punto A, iría a B para finalmente alcanzar C y seguir con su trayectoria horizontal. En los flujos viscosos reales la partícula, en su trayecto alrededor del cilindro, va perdiendo energía cinética debido a la fricción y no es capaz de alcanzar la presión necesaria para llegar a C, por lo que se desprende del cuerpo inmerso del fluido desviándose de su trayectoria ideal. Esta aportación supuso un paso importantísimo en el estudio del Efecto Magnus al permitir comprender la física que intervenía en él. No obstante todavía no se dispone de un modelo matemático exacto que tenga en cuenta todas las variables que influyen en este fenómeno y prediga su comportamiento, por lo que se hace imprescindible la adquisición de datos experimentalmente o medio de CFD. mediante simulación por 2.4 Estudios experimentales sobre el comportamiento del Efecto Magnus. Antes de analizar el efecto Magnus conviene observar las características de un flujo incidente perpendicular a un cilindro en reposo. Según el trabajo publicado por Kang y Choi [KANG99], el comportamiento del fluido presenta una gran dependencia del número de Reynolds del problema. Para bajos números de Reynolds ( 47), la estela que se forma tras el cilindro como consecuencia del desprendimiento de la capa límite está formada por una región de fluido en recirculación compuesta por dos torbellinos simétricos unidos al cilindro (Figura 4). Con el aumento del número de Reynolds los torbellinos aumentan su tamaño.
21 Introducción al efecto Magnus 12 Figura 4. Vórtices para bajo número de Reynolds [PRAN25]. Para números de Reynolds comprendidos entre el flujo se vuelve inestable y los vórtices se desprenden del cilindro en sus puntos superior e inferior, oscilando periódicamente; son los llamados Vórtices de Von Karman (Figura 5). Esto produce una fuerza de sustentación oscilante de la misma frecuencia y con valor medio nulo. Figura 5. Vórtices de Von Karman para (a) Re=54, (b) Re=65 y (c) Re=102[BIFA03]. Para números de Reynolds 200 el flujo se vuelve turbulento y tridimensional y los vórtices se desprenden de manera turbulenta e incontrolada.
22 Introducción al efecto Magnus 13 Figura 6. Desprendimiento turbulento de la estela para Re= 8 10 [BIFA03]. El desprendimiento de estos vórtices puede dar lugar a vibraciones indeseadas, acarreando problemas inesperados, por lo que es importante controlar su desprendimiento. La rotación necesaria para darse el Efecto Magnus, que modifica el desprendimiento de la estela tras el cilindro, es capaz de anular estas vibraciones, al hacer que los vórtices se desprendan siempre del mismo lado del cilindro. Kang y Choi también aportan interesante información sobre este aspecto. El control de la estela turbulenta es función del número de Reynolds y de la relación entre la velocidad de rotación del cilindro y la del fluido incidente. Para 47, a pesar de que, como ya se ha indicado, los vórtices no se desprenden, la rotación del cilindro retrasa el desprendimiento de la capa límite pudiendo llegar a evitarlo. Para números de Reynolds 47 se observa la existencia de un, dependiente del número de Reynolds, para el que las oscilaciones de la estela desaparecen. Como se puede apreciar en la Figura 7 el parámetro depende logarítmicamente del número de Reynolds, pudiendo considerarse constante e igual a 2 para valores de 200.
23 Introducción al efecto Magnus 14 Figura 7. Curva que determina el paso de flujo oscilatorio en la estela del cilindro a régimen estacionario. Medidas tomadas con desprendimiento de vórtices. ( Me ó [KANG99]. Se puede apreciar que un aumento en el número de Reynolds contribuye a la inestabilidad del flujo tras el cilindro, mientras que un incremento de la velocidad de rotación del cilindro tiende a estabilizarlo. Estas oscilaciones en el flujo implican oscilaciones en las fuerzas que intervienen en el Efecto Magnus, por lo que es importante tenerlas en cuenta y, en la medida de lo posible intentar evitarlas. A lo largo de las diversas investigaciones que se han llevado a cabo se ha podido constatar que las fuerzas producidas como consecuencia del efecto Magnus dependen de una gran cantidad de parámetros. A lo largo de este apartado se partirá de la bibliografía disponible para analizar, según el número de Reynolds, estas influencias en el valor de las fuerzas comprendidas en el Efecto Magnus Efecto Magnus a bajos números de Reynolds (Re<300). Los estudios a bajos números de Reynolds tienen la particularidad de que las fuerzas viscosas no son despreciables frente a las fuerzas de
24 Introducción al efecto Magnus 15 inercia. Es por esto por lo que el coeficiente de drag ( ) y en menor medida el coeficiente de lift ( ) son dependientes del número de Reynolds, tal y como observan Kang y Choi [KANG99] en simulaciones en dos dimensiones. Cuanto menor sea el número de Reynolds más importantes serán las fuerzas viscosas y mayor dependencia presentarán dichos coeficientes. Un ejemplo claro de esto se aprecia en la Figura 8. donde se muestran los valores de recopilados por Stojkovic [STOJ02], para el caso particular de muy bajos números de Reynolds (Re<1). Estos flujos se denominan creeping flow (Re<1), y su característica principal es que las fuerzas viscosas tienen gran importancia frente a las de inercia. En ella se puede ver cómo para un mismo φ, al disminuir el Re, se acerca a su valor ideal según la teoría de flujo potencial (sin desprendimiento de capa límite). Esto se debe a que las fuerzas de inercia son menores por lo que tarda más en desprenderse la capa límite.
25 Introducción al efecto Magnus 16 Figura 8. Coeficiente de lift ( ) en función de Re y de [STOJ02]. Para estudiar detenidamente estos coeficiente conviene diferenciar en la parte del lift y drag producidas por la presión ( y ) y la debida a las fuerzas de fricción y turbulencias o inducida ( y ). Por tanto: ; Ec. 1. Componentes de y debidas a la presión y las inducidas. Coeficiente de lift (C : como se observa en la Figura 10 la parte del lift debida a la fricción (C ) es menor que la producida por la presión (C ). Ambas aumentan con el aumento del parámetro φ (relación entre la velocidad de rotación del cilindro y la del fluido incidente), pudiéndose aproximar esta dependencia por una función lineal:
26 Introducción al efecto Magnus 17 Re 40; Re 60; Re 100; Re 160; C~2.57 φ; C~2.5 φ; C~2.48 φ; C~2.46 φ; C~2.31 φ; C~0.26 φ; C~2.29 φ; C~0.21 φ; C~2.31 φ; C~0.17 φ; C~2.33 φ; C~0.13 φ; Obsérvese que en las ecuaciones están escritas para los coeficientes promediados (C, C y. C Esto se debe a que, como ya se ha visto, por debajo de φ 2 los coeficientes se ven afectados por las turbulencias tras el cilindro y presentan oscilaciones. En la Figura 10 se representa la amplitud de las oscilaciones de en función del Re y de φ; obsérvese cómo las oscilaciones desaparecen al superar el φ crítico propio de cada Re representado en la Figura 7. Figura 9. Valor medio del coeficiente de lift ( ) en función de φ., Re=160;, Re=100;, Re=60;, Re=40;, teoría de flujo potencial. También se muestran las contribuciones debidas a la fricción ( ), la presión ( ) y el valor total (), [KANG99].
27 Introducción al efecto Magnus 18 Figura 10. Amplitud de las fluctuaciones del coeficiente de lift(c_l) en función de φ., Re=160;, Re=100;, Re=60;, Re=40;, teoría de flujo potencial, [KANG99]. Coeficiente de drag (C : este parámetro tiene la particularidad de que su componente producida por la fricción ( ) y la producida por la presión ( ) son del mismo orden de magnitud, por lo que presenta una mayor dependencia del número de Re. Al aumentar el número de Reynolds el y C tienden a disminuir. Respecto del parámetro φ (relación entre la velocidad de rotación del cilindro y la del fluido incidente), un aumento en él implica una disminución en el y un aumento del, resultando una disminución del coeficiente de drag total (C ) (Figura 11). Al igual que el, el C presenta valores fluctuantes asociados a los desprendimientos de capa límite, que cesan al superar el φ crítico (Figura 7).
28 Introducción al efecto Magnus 19 Figura 11. Valor medio del coeficiente de drag ( ) en función de., Re=160;, Re=100;, Re=60;, Re=40;, teoría de flujo potencial. También se muestran las contribuciones debidas a la fricción ( ), la presión ( ) y el valor total ( ), [KANG99]. Figura 12. Amplitud de las fluctuaciones del coeficiente de drag( ) en función de φ., Re=160;, Re=100;, Re=60;, Re=40;, teoría de flujo potencial, [KANG99].
29 Introducción al efecto Magnus Efecto Magnus a altos números de Reynolds ( ). Para números de Reynolds 10 de esta magnitud las fuerzas de inercia son mucho más importantes que las viscosas. Para poder estudiar este fenómeno y dado que no existe una explicación matemática de la física del Efecto Magnus, Bychkov [BYCH05] emplea una aproximación basándose en la teoría de flujo potencial y datos experimentales recogidos por él. En ella relaciona la fuerza de lift (Y) que experimenta un cilindro en rotación con el coeficiente de velocidades φ y el ángulo formado por los puntos en los que se desprende la capa límite respecto del centro del cilindro (Figura 13), por medio de la siguiente expresión. Γ Donde: Y es la fuerza de lift que experimenta el cilindro en rotación. es la densidad del fluido que envuelve al cilindro. V es la velocidad que tiene el fluido lejos del área de influencia del cilindro. Γ es la circulación real creada por el cilindro en rotación. L es la longitud del cilindro. La circulación real (Γ) alrededor del cilindro puede ser modelada a partir de la circulación obtenida a partir de la teoría de flujo potencial (Γ ). Γ Γ 2 Donde: R es el radio del cilindro. es la velocidad de rotación del cilindro. es un parámetro que modela el efecto de la separación de la capa límite. Según [ALDO88] está relacionado con los puntos superior e inferior
30 Introducción al efecto Magnus 21 del cilindro en que se desprende la capa límite ( ) (Figura 13) por medio de la siguiente ecuación., 2 El valor de 1 corresponde al caso 2 que es la situación en la que no hay desprendimiento de capa límite, descrita por la teoría de flujo potencial. Para casos reales, y dentro de los números de Reynolds manejados en este apartado ( 10 ), el valor máximo que puede tomar este parámetro es 0.5. Se pueden conseguir valores mayores acoplando discos a lo largo de la superficie del cilindro para aumentar la fricción entre fluido y disco, pero esto hace aumentar demasiado la potencia necesaria para mantener al cilindro rotando y por tanto deja de ser práctico. En la siguiente figura (Figura 13) se puede apreciar la diferencia entre el flujo que predice la teoría del flujo potencial (a) y el flujo real (b) alrededor de un cilindro. También se muestra el comportamiento que presenta la estela tras superar el 4 10 ; al desprenderse la estela del cilindro más tarde disminuye la fuerza de drag que experimenta el cilindro.
31 Introducción al efecto Magnus 22 Figura 13. Diagramas de flujo alrededor de un cilindro. (a) Teoría de flujo potencial. (b)flujo real con desprendimiento de capa límite. A,B son los puntos de estancamiento. S son los puntos en los que se desprende la capa límite. Los datos experimentales publicados sobre este fenómeno a altos números de Reynolds sugieren un comportamiento distinto para valores de ϕ<1 y ϕ>1; en especial los ensayos llevados a cabo en el ITAM SB RAS, por ser los más exhaustivos y las publicaciones de N.M. Bychkov [BYCH05]. Dado que los rangos de ϕ<1 no son relevantes para este proyecto y sólo se darían en periodos de tiempo breves en que arranca el aerogenerador, sólo se estudiará el caso de ϕ> Datos experimentales para ϕ>1. Para valores de ϕ>1 la capa límite se vuelve totalmente turbulenta y los parámetro de lift ( o ) y de drag ( o ) son prácticamente independientes del número de Reynolds. Sin embargo se ha observado la influencia de varios parámetros:
32 Introducción al efecto Magnus 23 Influencia de ϕ y C: con el fin de minimizar las turbulencias que se producen en los extremos del cilindro es común poner un disco a cada lado del cilindro para dirigir el flujo y evitar que escape por los bordes. C es una constante que relaciona el diámetro de este disco (d ) con el diámetro del cilindro (d). La Figura 14 muestra los datos presentados por N.M. Bychkov [BYCH05] para un cilindro de relación de aspecto 8 (L es la longitud del cilindro y d su diámetro) para números de Reynolds Figura 14.Valores experimentales de y de para un cilindro en rotación con 8 y C=1-3, A es la asíntota para 0.5, [BYCH05].
33 Introducción al efecto Magnus 24 En ella se puede apreciar cómo la incorporación de un disco en los extremos aumenta sustancialmente el, llegando a apróximarlo a su valor asintótico para valores de ϕ=3-4. Obsérvese que valores de C >3 no aportan una mejoría significativa. El también se ve afectado por este parámetro, ya que al utilizar los discos que disminuyen las turbulencias en los extremos del cilindro y con ello el. No obstante unos cilindros con diámetros excesivamente grandes implicarían un aumento excesivo de la fricción y el aumentaría. También destaca el crecimiento de y con el aumento de la velocidad de rotación (ϕ), especialmente la región comprendida entre ϕ=1-4, en donde la relación de y con ϕ puede considerarse lineal. Los valores de ϕ y C óptimos para maximizar el lift y minimizar el drag de un cilindro en rotación serían ϕ=4 y C=2. Influencia de la relación de aspecto del cilindro (λ L d): de acuerdo con los datos experimentales de N.M. Bychkov [BYCH05] representados en la (Figura 15) se deduce que tanto el como el dependen del parámetro λ del cilindro. Se aprecia que un aumento en λ influye en el C, especialmente para valores de ϕ>4, produciendo un incremento en su valor, que tiende a acercarse a su límite asintótico (A).
34 Introducción al efecto Magnus 25 Figura 15. Influencia de λ en el (a) y (b), para C=3 y valores de λ = 3(), 6(), 12( ), [BYCH05]. El C presenta una mayor dependencia de λ que el C y, al contrario que éste, tiende a disminuir conforme crece λ de manera significativa. Esto es debido a que la componente del C inducida (ver Ec. 1), es función de la siguiente fórmula. Puesto que >> se puede asumir que. Según se aprecia en la esta aproximación se corresponde con los datos experimentales obtenidos, puesto que el paso de λ=6 a λ=12 supone una reducción en el del 70% - 90% (Figura 15). Por tanto el valor óptimo de λ no debe ser menor de 12, ya que para este valor se aprecia en la Figura 15 que, en rangos óptimos de ϕ=3-4, aunque el C es prácticamente constante a partir de λ=8, el disminuye considerablemente conforme aumenta dicho parámetro.
35 Introducción al efecto Magnus 26 Figura 16. Coeficientes de (a) y (b) en función de λ y ϕ para valores de C=2 (,) y C=3 (), [BYCH05]. En la referencia WHIT04 se trata igualmente el Efecto Magnus y, aunque no se estudia tan exhaustivamente como en BYCH05, tiene la particularidad de aportar la gráfica mostrada en la Figura 17. Su principal particularidad es que los datos que recoge están evaluados para un modelo en dos dimensiones, en el que por tanto no se tienen en cuenta turbulencias tridimensionales. Debido a esto los valores de que recoge son ligeramente superiores. Esta gráfica será de gran utilidad para compararla con los resultados numéricos que se obtengan en las simulaciones numéricas de este proyecto.
36 Introducción al efecto Magnus 27 Figura 17.
37 Capítulo 3 INTRODUCCIÓN A LA DINÁMICA DE FLUIDOS COMPUTACIONAL
38 Introducción a la dinámica de fluidos computacional Historia de la dinámica de fluidos computacional A lo largo de la historia el hombre siempre se ha interesado por conocer el comportamiento de los fluidos. Los primeros estudios acerca de la mecánica de fluidos fueron llevados a cabo por Arquímedes (Figura 18), que sentó las bases de la hidrostática enunciando el principio que lleva su nombre. Con la llegada de la civilización romana tomo importancia el estudio de fluidos, en especial aplicados a la construcción de obras hidráulicas civiles como acueductos, canales, puertos y baños públicos. Con la llegada del Renacimiento se dio un nuevo impulso a esta Figura 18. Arquímedes de Siracusa (287B.C B.C.) ciencia con los estudios de ingenieros como Leonardo Da Vinci, quien en su tratado Del moto e misura dell'acqua expone estudios sobre el movimiento de agua (olas, remolinos, cascadas, chorros de agua, etc.). En el siglo XVII continuaron las aportaciones a la hidráulica con los estudios de científicos tan importantes como Torricelli (inventor del barómetro, Teorema de Torricelli), Pascal (Principio de Pascal) y Newton, quien aplicó sus tres leyes de la mecánica al estudio de los fluidos y analizó la fricción y la viscosidad de los mismos.
39 Introducción a la dinámica de fluidos computacional 30 Durante los siglos XVIII y XIX se sentarían las ecuaciones que hoy en día rigen la mecánica de fluidos gracias a la contribución de numerosos científicos. Algunos de los más destacados son O. Reynolds, S.-D. Poisson, J. W. Rayleigh, J. L. Poiseuille, M. M. Couette o D. Bernoulli que introdujo la famosa ecuación de Bernoulli. En esta época surgirían dos importantes aportaciones que son las bases de la actual mecánica de fluidos computacional (CFD): L. Euler enunció matemáticamente, en sus ecuaciones de Euler, la conservación de masa y momento para un cilindro no viscoso y aportó la teoría de flujo potencial. C. L. Navier (Figura 19) y G. G. Stokes (Figura 20) realizaron un importante avance incluyendo los efectos de la viscosidad en las ecuaciones de Euler desarrollando así las ecuaciones de Navier-Stokes (Ec. 2). Figura 19. Claude Louis Marie Henry Navier Figura 20. George Gabriel Stokes 2 3 Ec. 2. Ecuaciones de Navier-Stokes en coord. cartesianas y notación de subíndices. Debido a la complejidad de estas ecuaciones no se dispone de una solución analítica general y se desarrolla la mecánica de fluidos
40 Introducción a la dinámica de fluidos computacional 31 computacional aplicando métodos para transformar problemas continuos en modelos discretos aproximados. Tras las ecuaciones de Navier-Stokes, durante el siglo XX los esfuerzos se enfocaron en el estudio de la capa límite y el flujo turbulento. L. Prandtl ( ) enunció la teoría de la capa límite, realizó estudios sobre flujos compresibles y desarrolló las primeras teorías sobre ondas de choque supersónicas. T. Von Karman ( ) descubrió la calle de vórtices de Von Karman. G. I. Taylor ( ) avanzó en el estudio del flujo turbulento, propuso una teoría estadística para explicarlo desarrollando la microescala de Taylor que establece la mayor escala de longitud a la que la viscosidad del fluido afecta a la dinámica de los vórtices en los flujos turbulentos. A. N. Kolmogorov ( ) y G. K. Batchelor ( ) investigaron el campo del flujo turbulento aportando interesantes avances en el área. La primera solución numérica obtenida con métodos de CFD fue obtenida en 1933, por A. Thom y consistía en el cálculo de flujos a bajas velocidades que incidían en cilindros. El estudio lo dio a conocer en su publicación [THOM33]. No obstante el CFD era aún un método muy lento debido a la gran cantidad de cálculos que había que realizar. Muestra de ello es que veinte años más tarde, en 1953, el japonés Kawaguti obtuvo una solución para un problema similar utilizando una calculadora mecánica y tardando dieciocho meses trabajando veinte horas semanales. En la década de los sesenta la división teórica de la NASA situada en Los Álamos investigó la dinámica de fluidos computacional, aportando numerosos modelos numéricos que hoy en día se siguen usando, como el conocido modelo de k-.
41 Introducción a la dinámica de fluidos computacional 32 Durante los años setenta las principales investigaciones corrieron a cargo de un grupo de trabajo del Imperial College de Londres que bajo la dirección de D. Brian Spalding contribuyeron a la mejora de los códigos ya existentes, como el modelo k- y al desarrollo de nuevos modelos numéricos. La década de los ochenta supondría un gran salto para el CFD. En 1980 se publicó uno de los libros más importantes hasta la fecha sobre la materia, que daría pie a la creación de multitud de códigos. Se titula Numerical Heat Transfer and Fluid Flow y su autor es Suhas V. Patankar. A principios de esta década surgieron los primeros códigos de CFD comerciales, que tuvieron una gran aceptación y las grandes empresas los empezaron a utilizar en lugar de seguir creando sus propios códigos como habían hecho hasta entonces. Los avances en CFD y el aumento de la potencia de los ordenadores actuales han hecho de este método una herramienta imprescindible a la hora de resolver cualquier problema de dinámica de fluidos actual. La dinámica de fluidos computacional actual permite resolver numéricamente problemas que hasta hace poco eran prácticamente imposibles con los métodos matemáticos tradicionales. Esto implicaba la construcción de numerosos prototipos para probarlos y corregir errores de forma iterativa. Hoy en día el CFD permite realizar modelos matemáticos para estudiar la realidad, mucho más baratos de modificar que los prototipos. Esto ha supuesto una reducción de costes y tiempo de diseño. No obstante estos siguen siendo necesarios, pues los resultados que se obtienen no dejan de ser una aproximación matemática de la realidad. No obstante, debido a la complejidad del proceso, para obtener resultados fiables es imprescindible conocer el programa y el método de discretización que se está utilizando, los pasos a seguir para definir el
42 Introducción a la dinámica de fluidos computacional 33 modelo correctamente, y elegir existentes adecuados para modelar los problemas. el mallado y modelos de ecuaciones 3.2 Métodos de discretización del problema. Como se ha mencionado antes el comportamiento de los fluidos se rige por un conjunto de ecuaciones (ecuaciones de Navier-Stokes) en derivadas parciales que no tienen solución analítica para la mayoría de los problemas de ingeniería. En un problema continuo las variables a estudiar están definidas para cada punto del dominio. El objetivo del CFD es pasar de ese dominio continuo a uno discreto, mediante un proceso de mallado, y así resolver las variables solo en los puntos de la malla. Esto significa pasar a resolver un número finito de incógnitas. Adaptando las ecuaciones del problema continuo al discreto se obtiene un sistema de ecuaciones algebraicas y n incógnitas, que son los valores de las variables discretizadas en cada punto de la malla. Los valores de las variables en los demás puntos se obtienen por interpolación de la malla. Por ejemplo en un problema continuo en el que se quiere hallar el valor de la velocidad en cada punto del dominio, ante la imposibilidad de resolver las ecuaciones, se divide el dominio (discretizado) para así tener un número finito de variables y pasar a resolver ecuaciones algebraicas en lugar de las ecuaciones en derivadas parciales iniciales. Dominio continuo Dominio discreto Figura 22. Sistema de ecuaciones en derivadas parciales + variables continuas [BHASK]. Figura 21. Sistema de ecuaciones algebraicas + variables discretas [BHASK].
43 Introducción a la dinámica de fluidos computacional Métodos de discretización de las diferencias finitas. Para comprender el funcionamiento de este método se resolverá un problema sencillo pero ilustrativo: Considérese un dominio unidimensional en el que se desea resolver el valor de una variable u a lo largo de él. El comportamiento de la variable viene definido por la siguiente ecuación: Ec. 3. Dependiendo del valor del exponente m el problema es lineal ( 1) o no lineal ( 1). Primero se resolverá el problema lineal y posteriormente se analizarán las diferencias que implica la no linealidad. El primer paso es discretizar el dominio, para este problema se ha dividido en 3 partes, teniendo así cuatro nodos equidistantes. Figura 23. dominio discretizado La ecuación que rige el comportamiento de nuestra variable sigue siendo la misma, por lo que para cada nodo será: 0 Siendo i el número de nodo. La principal característica de este método es la obtención de una expresión algebraica a partir de la ecuación diferencial aproximando por las series de Taylor el término diferencial: ΟΔ
44 Introducción a la dinámica de fluidos computacional 35 Despejando: Ec. 4. Aproximación de la derivada parcial. Ο es el error de truncamiento de la serie, que depende de, por lo que cuanto menor sea el tamaño de las celdas de la malla menor error se cometerá. Sustituyendo la Ec. 4 en la Ec. 3 se obtiene: 0 Con lo que se consigue un sistema algebraico para resolver el problema: Para resolver el problema basta con hallar la solución al sistema, invirtiendo una matriz de 4 x 4. No obstante en problemas reales suele haber cientos de miles de nodos, que implican un sistema con multitud de ecuaciones e incógnitas. Resolver estos sistemas implica realizar operaciones de inversión de matrices de muchísimo tamaño, lo que supondría una gran cantidad de memoria. Para evitar este problema las ecuaciones se reordenan para obtener los valores de las variables de un nodo en función de sus valores en los nodos adyacentes: 1
45 Introducción a la dinámica de fluidos computacional 36 En la primera iteración sólo se tienen valores para los nodos adyacentes, a los demás se les considera 0. En sucesivas iteraciones los valores irán convergiendo a la solución exacta. Esto no deja de ser una aproximación del resultado que se obtendría de invertir la matriz del sistema, por lo que se iterará hasta que la diferencia entre los resultados de dos iteraciones consecutivas no exceda del error admitido en la solución. lineal: Se considera ahora el mismo problema, pero con una ecuación no 0 Ec. 5. Tras sustituir la ecuación la Ec. 4 en la Ec. 5 se obtiene: 0 Ec. 6. Para linealizar la ecuación hay que apoyarse en un valor solución supuesto para cada nodo ( ) y después iterar para ir modificando tanto este valor como el valor de la solución real, hasta que la diferencia entre ambos ( ) sea inferior al error admitido en nuestra solución. Elevando la ecuación al cuadrado y reorganizándola queda: 2 Suponiendo se puede despreciar el término : Sustituyendo esta ecuación en la inicial (Ec. 6) desaparece la incógnita elevada al cuadrado, linealizando la ecuación, que queda: 2 0
46 Introducción a la dinámica de fluidos computacional 37 Para la iteración inicial se da un valor supuesto a. En las siguientes iteraciones este valor será el resultado ( ) obtenido en la iteración anterior. La iteración continuará hasta que ambos valores converjan y sea inferior al error asumido en la solución. Al igual que en el caso lineal, para evitar invertir matrices de gran tamaño se despeja la solución de cada nodo en función de la de los adyacentes y al iterar, converge hacia la solución real. Durante la iteración: 1 2 Las principales características de este método son su formulación sencilla y que las mallas curvas deben pasarse a coordenadas cartesianas para que las ecuaciones diferenciales puedan ajustarse a una malla cuadrática cartesiana Métodos de discretización de los volúmenes finitos. Para aplicar este método las ecuaciones que rigen las incógnitas del problema tienen que estar escritas en su forma integral. La idea principal consiste en discretizar el modelo continuo mediante una malla y aplicar las ecuaciones de conservación (de la masa, del momento, de la energía ) a cada celda por separado. Para ilustrarlo supongamos un flujo estacionario incompresible en un dominio continuo de 2 dimensiones al que le queremos aplicar la ecuación de conservación de la masa (Ec. 7): 0 Ec. 7.Ecuación de la conservación de la masa en su forma integral Esta integral se realiza a lo largo de la superficie de cada celda. El vector es normal a la superficie.
47 Introducción a la dinámica de fluidos computacional 38 El método asigna un valor de la variable a cada celda, y aproxima el flujo que entra por cada cara de una celda como la longitud de esa cara por un valor constante de la variable en ella (Figura 24.). La velocidad en el centro de cada celda será. Aplicando la ecuación de conservación de la masa a la celda se obtiene: 0 Los valores, para un punto cualquiera se obtienen interpolando los valores de centro de celda de las celdas adyacentes. Al igual que en otros métodos partiendo de los valores de las condiciones iniciales del problema se interpola hasta que la diferencia entre los resultados de dos interpolaciones consecutivas sea inferior al error admitido en la solución. Figura 24. Este método es el utilizado por el programa FLUENT para resolver los diferentes modelos de ecuaciones que incluye. Una ventaja de este método es que las geometrías complicadas y las mallas no estructuradas se pueden tratar fácilmente sin necesidad de realizar ninguna transformación de coordenadas.
48 Introducción a la dinámica de fluidos computacional Métodos de discretización de los elementos finitos. Este método tiene la particularidad de ser matemáticamente más complejo que los anteriores y es necesario mayor rigor matemático en su aplicación. El método de los elementos finitos resuelve un problema determinado por ecuaciones diferenciales y condiciones de contorno aproximándolo a un problema en forma matricial. Determina la solución exacta solo en ciertos puntos del dominio (los nodos de la malla con la que se discretiza el problema) y obtiene la solución del resto de los puntos por extrapolación. Cuanto más fina sea la malla más precisa será la solución que proporciona este método. Para ilustrar el procedimiento que sigue el método de los elementos finitos se resolverá el siguiente problema unidimensional: Hay varias formas de implementar el método de los elementos finitos, en este caso utilizaremos el método de Galerkin Standard (SGM). En primer lugar se discretizará el dominio continuo (Figura 25 (a)) en dos elementos, dando así lugar a tres nodos (Figura 25 (b y c)),
49 Introducción a la dinámica de fluidos computacional 40 Figura 25. (a) Dominio continuo con condiciones de contorno ( 0, 1). (b) Nodos globales del dominio discretizado. (c) Elementos locales. (d) Funciones de interpolación locales [CHAN02]. El siguiente paso será asumir que la variable a resolver es función lineal de x (Ec. 8.): Ec. 8. Si se despeja y en función de los valores de u en los nodos ( y ) de un elemento cualquiera (e) resulta la siguiente ecuación. 1 1,2 Ec. 9. Donde representa el valor de u en el nodo local N en el elemento (e). son las denominadas funciones de interpolación o de forma del dominio local (de un elemento). 1 ; 0 1 Ec. 10.
50 Introducción a la dinámica de fluidos computacional 41 Estas funciones son propias de cada elemento. Cada una de ellas está asociada a uno de los nodos, tomando el valor 1 en él y 0 en el otro (Figura 25 (d)). El siguiente paso será minimizar el error que supone la utilización de las funciones de forma. Para ello los errores se proyectan ortogonalmente en el subespacio formado por las funciones de interpolación (Ec. 10.) y se igualan a 0 dichas proyecciones., 2 0 Integrando por partes se obtiene. Ec. 11. Sustituyendo la Ec. 9. en la Ec. 11. se obtiene. Ec Esta ecuación es conocida como la formulación débil del problema, puesto que la segunda derivada de la ecuación diferencial se ha transformado en una derivada primera. Es importante destacar que la derivada / del término ya no es la derivada de la variable, sino que es una condición de contorno de Neumann (constante). Igualmente tampoco es ya una función continua en el espacio, sino una función
51 Introducción a la dinámica de fluidos computacional 42 test de la condición de contorno de Neumann. Está función tomara el valor 1 si la condición de contorno de Neumann se aplica en el nodo N, en caso contrario su valor es 0. Además es importante resaltar que, dependiendo de si nos acercamos al punto con condición de contorno de Neumann por la izquierda o por la derecha, el valor de la variable en dicho punto estará aumentando (derivada con signo positivo) o disminuyendo (derivada con signo negativo). Al dotarle de signo a la condición de contorno de Neumann obtenemos lo siguiente. cos, cos Utilizando una notación simplificada la Ec. 12. queda así. 0, 1,2 Esta ecuación es un sistema algebraico de ecuaciones locales, es decir, propias de cada elemento de la malla, y su resolución proporciona los valores de la variable (que en este problema es ) en los nodos de dicho elemento (e). Los términos que componen el sistema son los siguientes (se particularizan sus valores para los datos del problema que se está tratando): Matriz de rigidez, de difusión o de viscosidad ( ): está relacionada con los fenómenos físicos asociados a la derivada segunda de la ecuación diferencial.
52 Introducción a la dinámica de fluidos computacional Vector fuente ( ): Vector de condiciones de contorno de Neumann ( ): cos Estos sistemas locales, propios de cada elemento (e) pueden ser formulados de acuerdo a la notación global de los nodos (, 1,2,3 (Figura 25(d))) y ensamblados formando un sistema global cuya solución será el valor de la variable en cada uno de los nodos de la malla. Para ensamblar los sistemas basta con sumar las contribuciones al nodo global de los elementos adyacentes a él. Para el problema que se está resolviendo los términos del sistema global son
53 Introducción a la dinámica de fluidos computacional En el problema que se está resolviendo 0 cos 0 cos 0,, 0, pues en el enunciado no se da ninguna condición de contorno de Neumann. En el caso que nos ocupa el sistema algebraico global es Resolviéndolo se obtienen los valores de, y, que representan el valor de u en los nodos globales, y con la Ec. 9. se extrapola la solución al resto del dominio. Las geometrías complicadas y las mallas no estructuradas tampoco requieren de transformaciones de coordenadas en este método. 3.3 Mallado. Un aspecto muy importante de la dinámica de fluidos es la discretización del dominio fluido, puesto que las propiedades que tenga la malla influirán tanto en el tiempo de preparación de la simulación, como en la duración del proceso de resolución y la fiabilidad de la solución.
54 Introducción a la dinámica de fluidos computacional Tipos de malla: Las características más importantes de una malla son la forma que adquieren sus celdas y la manera en que estas se reparten. Atendiendo a la forma de las celdas, se distinguen diversas opciones, cada una apropiada para un tipo de geometría y problema: Problemas en 2D: se utilizan celdas de forma de triángulo y cuadrilátero. Problemas en 3D: los volúmenes del mallado pueden tener forma de tetraedro, hexaedro, prisma, pirámide y diversos poliedros. Figura 26. Tipos de celdas [FLUE09] No obstante también es corriente utilizar mallados híbridos (Figura 27) realizando combinaciones de triángulos y cuadriláteros para los problemas de 2D y combinando hexaedros y tetraedros en los problemas de 3D, lo que permite aprovechar las ventajas de cada tipo de celda a lo largo de la malla.
55 Introducción a la dinámica de fluidos computacional 46 Figura 27. Malla híbrida triángulos / cuadriláteros. Otra propiedad destacada de la malla es la distribución que siguen las celdas a lo largo de la misma. Se pueden ordenar de manera estructurada siguiendo un patrón o geometría (Figura 28) o, por el contrario, optar por una malla no estructurada (Figura 29). Figura 28. Malla estructurada a partir de cuadriláteros. Figura 29. Malla no estructurada. Además de mallar el dominio fluido como un solo conjunto también es común dividir el problema en varios bloques (Figura 30) para poder definir en cada uno de ellos propiedades diferentes como el patrón que sigue la malla, la forma de las celdas, o los parámetros de la simulación.
56 Introducción a la dinámica de fluidos computacional 47 Figura 30. Malla estructurada en varios bloques. Mallado cuadrado [FLUE09] Elección de malla adecuada: Es importante tener presente durante este proceso las ventajas y desventajas de cada una de las características expuestas anteriormente y sus implicaciones en la resolución del problema. A la hora de elegir una malla hay que tener en cuenta diversos aspectos: Tiempo de preparación del modelo. Muchos de los problemas en ingeniería implican geometrías complicadas o grandes. Esto hace que realizar una malla estructurada (Figura 28) sea muy trabajoso y prácticamente imposible debido al tiempo que se tardaría en su construcción. Esta es una de las principales causas de optar por mallas no estructuradas. Un error en el que se puede caer con la utilización de mallas estructuradas es intentar simplificar la geometría en demasía para poder ahorrar tiempo en la fabricación de la malla Gasto computacional. Está directamente relacionado con el número de celdas de la malla. A la hora de discretizar un problema las mallas no estructuradas necesitan menos celdas que las que sí lo son. Esto se debe a que al no estar
57 Introducción a la dinámica de fluidos computacional 48 dispuestas siguiendo un patrón se favorece el poder agrupar celdas en las regiones del fluido en que sea necesario más precisión y disminuir su densidad en las áreas menos interesantes. En las mallas estructuradas la disposición de las celdas está condicionada al patrón que deben seguir, lo que hace que en ocasiones se coloquen gran cantidad de celdas en partes del dominio en que no son necesarias. Respecto al tipo de celdas cabe destacar que las de forma de cuadrilátero (2D)/hexaedro (3D) permiten una relación de aspecto mayor que los triángulos/tetraedros puesto que si se alargase en exceso un triángulo, resultaría una celda muy asimétrica y deformada, que afectaría a la precisión de la solución. Las mallas basadas en poliedros necesitan aún menos celdas, favorecen la convergencia y suponen menor coste computacional, pero son más groseras y, por tanto, más imprecisas Difusión numérica: Es un término matemático que hace referencia a los errores de truncamiento debido a los redondeos realizados en las aproximaciones y cálculos. Como consecuencia de éste fenómeno aparece una falsa difusión de los parámetros de la simulación, equivalente a aumentar su coeficiente de difusión real. Dicho error hace más importante cuanto menor es la difusión real en el problema a estudiar. Se sabe que este fallo es inversamente proporcional al refinamiento de la malla, por lo que para minimizarlo se puede aumentar el número de celdas en la malla. También se consigue disminuir la falsa difusión cuando el fluido está alineado con la malla. Por esto conviene emplear mallas estructuradas a base de cuadrados en problemas de 2D y hexaedros para 3D Comprobación de la malla. Dado el importante papel que juega la malla en la precisión y estabilidad del modelo es importante fijarse en algunos de sus atributos,
58 Introducción a la dinámica de fluidos computacional 49 que definen su calidad: distribución de los nodos, tamaño de celdas, suavidad en la transición de un tipo de celdas a otras, asimetría y ubicación de las celdas cercanas a paredes Transición. Al pasar del dominio continuo al discreto se introducen errores de truncado al aproximar las derivadas parciales de las ecuaciones. Estos errores se acentúan cuando hay un cambio brusco de área (2D) o volumen (3D) entre celdas adyacentes. Una posible solución a este inconveniente es refinar la malla en los puntos en que se dé esta situación Forma de la celda: Asimetría: se define como la diferencia de forma entre una celda de la malla y la celda equilateral de volumen equivalente (un cuadrilátero idealmente tendría ángulos de 90º). ó ó í ó ó Este parámetro debe ser siempre inferior a 0.95 con una media de 0.3. Valores superiores podrían dar problemas de convergencia de la solución y afectar a la precisión de la misma. Coeficiente de aspecto: define lo estirada que está una celda. Como ya se ha mencionado anteriormente, para determinados casos como flujos a través de tubos largos y finos, utilizar celdas con un coeficiente de aspecto alto permite resolver un problema necesitando menor cantidad de celdas sin sacrificar la precisión en la solución. Para no afectar a la estabilidad este coeficiente no debe ser superior a 35:1, aunque conviene evitar coeficientes por encima de 5:1 en las partes del fluido alejadas de la capa límite y 10:1 para celdas cuadriláteras/hexaédricas en celdas dentro de la capa límite.
59 Introducción a la dinámica de fluidos computacional 50 Coeficiente de deformación: la deformación de las celdas indica la deformación lateral de la celda (Ej.: La que se produce en un volumen cuando se le somete a un esfuerzo cortante). Se mide mediante un coeficiente adimensional. Para problemas de 2 dimensiones este coeficiente es el valor resultante de la multiplicación escalar del vector que une el centroide de la celda con el centro de cada segmento que delimita su perímetro, por el vector perpendicular a cada uno de estos segmentos. En los casos de 3D el coeficiente se obtendrá multiplicando escalarmente el vector que une el centroide de la celda con el centro de cada una de las caras laterales del volumen (celda), y el correspondiente vector área de cada cara. 1 La celda ideal tendrá un coeficiente igual a 0, mientras que para las peores será cercano a 1. Como norma general este coeficiente no deberá ser superior a 0.95 para celdas en forma de tetraedro y 0.99 para el resto Densidad de nodos y su agrupamiento: Con el mallado del problema se está pasando de un medio continuo a uno discreto, por lo que la validez de la solución que se obtenga, dependerá mucho de la cantidad de nodos en la malla y su distribución. Una baja densidad en partes críticas del dominio del fluido, como pueden ser aquellas en que haya que predecir el desprendimiento de la capa límite en situaciones con gradientes de presión adversos, introducirá errores en la solución final. Nótese que esta situación se dará en algunos puntos del problema que se va a simular en este proyecto, por lo que se prestará especial atención a la distribución de los nodos cercanos a la pared del cilindro en rotación.
60 Introducción a la dinámica de fluidos computacional 51 Es necesario prestar especial atención en las regiones cercanas a paredes, donde se estudian las capas límite, puesto que las celdas en estas regiones desempeñan un papel muy importante en el cálculo. La correcta distribución de las celdas será clave para obtener una solución óptima. Dependiendo del tipo de flujo que se dé en el problema a tratar las mallas deberán tener, en dichas zonas, unas determinadas propiedades. Estudio de capas límite en flujos laminares: para controlar, en estos casos, que el número de celdas cercano a una pared es suficiente para el cálculo de capas límite de forma fiable se emplea, de forma general, el siguiente criterio: Donde: y u flujo. capa límite. 1 Distancia a la pared desde el centroide de la celda adyacente. Velocidad del flujo no perturbado por el objeto inmerso en el Viscosidad cinemática del fluido. Distancia a lo largo de la pared desde el punto de inicio de la Estudio de capas límite en flujos turbulentos: en estos casos, para comprobar la correcta distribución de las celdas, se utiliza el parámetro adimensional y (Ec. 13.), que mide la distancia entre la capa límite y el centroide de la celda más cercana a ésta. Ec. 13. Donde es la viscosidad de fricción.
61 Introducción a la dinámica de fluidos computacional 52 La ley logarítmica que sigue una capa límite en equilibrio en un flujo totalmente desarrollado proporciona los límites superior e inferior para dicho parámetro. Para estudiar la capa límite en flujos turbulentos se puede optar por calcularla mediante modelos de turbulencia (Figura 32), o bien por modelarla mediante unas funciones que describen el comportamiento del fluido dentro de capas límite, denominadas funciones de pared, y que hoy en día son comúnmente utilizadas en los software comerciales de CFD (Figura 31). Figura 32. Cálculo de capa límite mediante modelos de turbulencia [FLUT09]. Figura 31. Modelado de capa límite mediante funciones de forma [FLUT09]. Dependiendo del método elegido la malla deberá tener unas características determinadas. Para resolver la capa límite mediante modelos de turbulencia es necesario acumular un número suficiente de celdas dentro de la capa límite como para poder calcularla con suficiente precisión. Para ello el parámetro y debera tener valores de y 5.
62 Introducción a la dinámica de fluidos computacional 53 En el modelado de la capa límite mediante funciones de pared es importante tener un espacio suficiente entre la pared y la primera celda de la malla para que dichas funciones puedan modelar correctamente la capa límite. De otra manera se empezarían a utilizar los modelos de turbulencia demasiado cerca de la pared, sin haber dado tiempo a las funciones de pared a modelar la totalidad del desarrollo de la capa límite. No obstante es conveniente tener algunas celdas dentro de la capa límite para asegurar la continuidad entre ésta y el flujo calculado con los modelos de turbulencia. Para cumplir con esto requisitos, el y de la malla debe estar situado idealmente entre 30 y 300, prefiriéndose valores cercanos a 30. Se pueden usar funciones de pared incluso con y 11.2 aunque se pierde precisión en las capas límite. Por todo esto, como norma general los valores del de cualquier malla serán de 5, para los casos en que se resuelva la capa límite mediante modelos de turbulencia, y de para los casos en los que se modele mediante funciones de pared. Es conveniente evitar por tanto la región de comprendida entre En flujos turbulentos la densidad de los nodos y su agrupación son también parámetros muy importantes para la correcta simulación de la totalidad del problema, no sólo de las capas límite. Debido a la fuerte interacción entre el flujo medio y sus turbulencias es importante tener una buena resolución de malla para obtener soluciones fiables, ya que este tipo de flujo es todavía más sensible a la distribución de los nodos que el flujo laminar. La densidad de nodos, sobre todo en simulaciones 3D, estará limitada por la memoria y capacidad de cálculo disponibles y la máxima duración admisible en la simulación.
63 Introducción a la dinámica de fluidos computacional Modelos de turbulencia. Los flujos turbulentos se caracterizan por tener campos de velocidad inestables. Estas fluctuaciones a su vez están divididas en escalas de turbulencia (Figura 33), según sea el tamaño de los vórtices que las forman. El cálculo de estas turbulencias a lo largo de todas sus escalas de manera exacta supondría, en los casos prácticos de ingeniería, un gasto computacional muy alto. Figura 33. Escalas de turbulencia. Existen varios planteamientos para simplificar el problema y poder disminuir el gasto computacional, que difieren entre sí en las aproximaciones o suposiciones que se asuman en el tratamiento de las ecuaciones que definen la física a estudiar. Todos estos métodos introducen nuevas variables a las ecuaciones que rigen el comportamiento de los fluidos (excepto el modelo DNS, apartado ), que necesitan ser modelados para llegar a una solución correcta. Las diferentes formas de tratar estas nuevas variables y resolver las ecuaciones son los llamados modelos de turbulencia.
64 Introducción a la dinámica de fluidos computacional Reynolds-averaged Navier-Stokes (RANS): La estrategia de este método se basa en dividir las variables turbulentas implicadas en el problema en un término promediado y otra parte fluctuante (Figura 34). Figura 34. Valor en función del tiempo de la velocidad en un punto, y su valor promediado (línea punteada), [[WILC98]. En el caso de la velocidad, cada una de sus componentes se divide de la siguiente manera. Ec. 14. Parte promediada ( ) y parte fluctuante ( ) de la componente de la velocidad. Igualmente cualquier variable escalar se divide del mismo modo. Ec. 15. Parte promediada ( ) y parte fluctuante ( ) de la componente de la velocidad. La solución al problema de flujo turbulento se obtiene resolviendo sólo la parte promediada de las variables mediante las ecuaciones de Navier-Stokes (Ec. 2) y modelando la totalidad de las turbulencias. Para ello se introducen las ecuaciones Ec. 14 y Ec. 15 en la ecuación de la conservación de la masa y en la de la conservación del momento, que
65 Introducción a la dinámica de fluidos computacional 56 forman las ecuaciones de Navier-Stokes. Tras esto se promedian las ecuaciones respecto al tiempo, llegando a las ecuaciones promediadas de Navier-Stokes (Reynolds Averaged Navier-Stokes ecuations, RANS) Ec. 16. Ecuaciones promediadas de Navier-Stokes. Para una mayor simplicidad en la escritura de las ecuaciones se han omitido el símbolo de promediado en las variables ( ). Estas ecuaciones, que están ahora en función de los valores medios de las variables promediados en el tiempo, presentan ahora un término que es el producto de dos componentes fluctuantes de la velocidad. Aunque el promedio en el tiempo de las componentes fluctuantes por separado ( y ) es cero por definición, el producto de dos fluctuaciones promediado no tiene por qué ser cero y es necesario modelarlo. Una aproximación muy empleada es la llamada hipótesis de Boussinesq. 2 3 Esta hipótesis tiene como principal ventaja el bajo coste computacional que implica utilizarla. Los modelos de turbulencia basados en el método RANS difieren entre sí en la forma en que modelan las turbulencias del problema Modelo Spalart-Allmaras: Fue enunciado por Spalart y Allmaras en 1992 [SPAL92]. Fue pensado para la resolución de aplicaciones aeronáuticas con flujos que involucraran capas límite. Proporciona buenos resultados para capas
66 Introducción a la dinámica de fluidos computacional 57 límite sujetas a gradiente de presión adversos. Esto lo haría útil en el problema que se aborda en este proyecto, sin embargo solo es efectivo a bajos números de Reynolds. Dado que los números de Reynolds manejados son del orden de 10 o mayores se deberá optar por otro modelo. El modelo Spalart-Allmaras es relativamente reciente, por lo que todavía no está contrastada su validez para todos los tipos de flujo estudiados en ingeniería. Últimamente está aumentando su utilización en el estudio de las turbomáquinas. Está basado en el método RANS descrito. Es un modelo de una sola ecuación [Ec. 17] que resuelve un modelo de la ecuación del transporte para la viscosidad cinemática de los torbellinos. 1 Ec. 17. Ecuación del transporte para el modelo Spallart-Allmaras En esta ecuación el término modela la producción de turbulencias y su destrucción debido a la amortiguación viscosa y a la influencia de la pared en las regiones cercanas a ella. Los términos y son constantes propias del modelo y es un parámetro que se define para cada problema a resolver Modelos : Modelo k ε estandar: desde que Launder y Spalding enunciaran este modelo [LAUN72] en 1972, se ha convertido en uno de los más utilizados en el análisis de flujos turbulentos en ingeniería. Se caracteriza por su robustez, economía de recursos computacionales y una razonable precisión para una amplia gama de flujos turbulentos. Se trata de un modelo semi-empírico compuesto por dos ecuaciones de transporte que rigen el comportamiento de la energía cinética de las turbulencias (k) y su tasa de disipación (ε). La ecuación del transporte de
67 Introducción a la dinámica de fluidos computacional 58 k proviene de la original (ecuación del transporte de Reynolds), mientras que la de ε fue obtenida por deducciones físicas y varía ligeramente de su expresión matemática exacta. En la resolución de este modelo se asume que el flujo es totalmente turbulento (altos números de Reynolds) y que los efectos de la viscosidad molecular son despreciables, por lo que sólo será válido para este tipo de flujos. Donde: Ec. 18: Ecuaciones para el modelo estandar. representa la generación de la energía cinética de las turbulencias debida a los gradientes de velocidad. representa la generación de la energía cinética de las turbulencias debido a la flotabilidad. representa la contribución de la fluctuación de la dilatación en flujos turbulentos compresibles debido a la disipación de la energía., y son constantes del modelo. y son los números de Prandtl turbulentos para y respectivamente. y s on parámetros que se definen para cada problema a resolver. El modelo de turbulencia ha sido objeto de múltiples estudios. En ellos se han desarrollado diversas mejoras dando lugar a nuevos modelos que, partiendo de las ecuaciones del modelo, tienen mejor
68 Introducción a la dinámica de fluidos computacional 59 comportamiento. Algunos de ellos son: el modelo RNG [YAKH86] o el modelo realizable [SHIH95]. Modelo RNG k ε: tiene un término adicional en la ecuación que modela la disipación de energía (ecuación que rige a ε) que mejora la precisión en la solución de flujos cuya dirección se modifica rápidamente (rapidly strained flows). El efecto de los vórtices es tenido en cuenta en el modelado de la turbulencia, mejorando la solución para flujos con vórtices. Utiliza una fórmula analítica para el cálculo de los números de Prandtl, mientras que en el modelo k ε estandar se modelan a partir de constantes que hay que introducir. Se han incorporado mejoras que permiten obtener buenas soluciones también para bajos números de Reynolds. Modelo k ε realizable: incorpora una nueva formulación de la viscosidad turbulenta y presenta mejoras en la ecuación que gobierna la disipación de la energía (ε). Con respecto al modelo simple presenta mejores resultados en el cálculo de la dispersión de chorros de fluido y en general problemas que impliquen flujos en rotación y capas límite con gradientes de presión adverso, desprendimiento y recirculación de fluido Modelos : Modelo k ω estandar: es un modelo empírico con dos ecuaciones (Ec. 2 y Ec. 19) basado en el modelo k ω de Wilcox [WILC98], que rigen el transporte de la energía cinética de las turbulencias (k) y su tasa específica de disipación (ω), que puede considerarse como la proporción entre ε y k. Se caracteriza por simular la dispersión de flujos libres (que no están delimitados por paredes), como por ejemplo la propagación del humo al salir de una chimenea. Este tipo de flujos son muy similares al comportamiento de las estelas en la región lejana al objeto que las crea,
69 Introducción a la dinámica de fluidos computacional 60 mezclas de flujos estratificados, chorros de fluido planos, circulares y axiales, por lo que también se utiliza este modelo para su cálculo. También es aplicable a modelos con capas límite. Sus ecuaciones son: Donde: Γ Γ Ec. 19: ecuaciones del modelo k-w estándar. representa la generación de la energía cinética de la turbulencia debido a los gradientes de la velocidad promediada. representa la generación de (tasa específica de disipación de la energía cinética de la turbulencia). Γ y Γ representan la disipación de debido a la turbulencia. y son parámetros que se definen para cada problema a resolver. Modelo SST k ω: a partir del modelo k ω estándar ha surgido este nuevo modelo [MENT94], que incorpora varios cambios. Se han introducido modificaciones en la formulación de la viscosidad turbulenta. En las zonas en que hay influencia de capas límite se emplea el modelo k ω estándar y posteriormente, en las zonas alejadas, cambia a una versión del modelo k ε adaptada para altos números de Reynolds. Estas modificaciones le hacen más preciso para un mayor rango de flujos turbulentos, como por ejemplo perfiles de alas, flujos con gradiente de presión adversos u ondas de choque transónico. Las ecuaciones de este modelo son: t ρk ρku x k Γ x G x Y S
70 Introducción a la dinámica de fluidos computacional 61 Donde: t ρω ρωu x ω Γ x G x Y D S Ec. 20: ecuaciones del modelo SST k-ω representa la generación de la energía cinética de la turbulencia debido a los gradientes de la velocidad promediada. representa la generación de (tasa específica de disipación de la energía cinética de la turbulencia). Γ y Γ representan la disipación de debido a la turbulencia. y son parámetros que se definen para cada problema a resolver. D modela la difusión cruzada (cross diffusion) Reynolds Stress Model (RSM): Es un modelo mucho más elaborado que los anteriores. Su utilización implica la resolución de cinco ecuaciones de transporte adicionales para 2D y siete ecuaciones en 3D. En él se cambia la hipótesis de la viscosidad turbulenta isótropa, por la resolución de las ecuaciones de transporte de las tensiones de Reynolds (tensor de tensiones de Reynolds), y una ecuación para evaluar la tasa de disipación. El modelo predice con más precisión casos como los efectos de curvatura de las líneas de corriente, remolinos y variaciones bruscas de los esfuerzos cortantes. Esto lo convierte en un modelo muy potente a la hora de resolver flujos complejos, especialmente aquellos que presentan algún tipo de anisotropía en el comportamiento de las tensiones de Reynolds. La ecuación de transporte empleada en este modelo es:
71 Introducción a la dinámica de fluidos computacional 62 ó, ó, ó ó ó 2 ó ó 2Ω ó é ó Ec. 21:Ecuación del transporte para el modelo Reynolds stress transport (RSM). En esta ecuación los términos de la ecuación C, D,, P y F se pueden resolver sin necesidad de modelarlos. Sin embargo D,,G, y son términos que necesitan ser modelados. Este hecho es considerado la principal causa de errores de falta de precisión. En la bibliografía [LAUN75] se da información más detallada sobre este modelo Direct numerical simulation (DNS). Los flujos turbulentos se componen de torbellinos de diferentes tamaños (clasificación según una escala de longitud) y duración (clasificación según escala de tiempos). En teoría es posible resolver, sin necesidad de utilizar modelos, todo el espectro de vórtices utilizando un método conocido como direct numerical simulation (DNS). Sin embargo
72 Introducción a la dinámica de fluidos computacional 63 los requerimientos computacionales que se necesitan son proporcionales al cubo del número de Reynolds turbulento ( ), por lo que para números medianamente altos se hace imposible su utilización Modelo LES (Large Eddy Simulation). La filosofía de este método se basa en evitar resolver todas las escalas de turbulencia (Figura 33).Los torbellinos grandes, comparables en tamaño al flujo principal, son resueltos directamente mientras que sólo los de menor tamaño se modelan. De esta manera disminuye el error introducido con los modelos de turbulencia. Los principios básicos de este modelo son: Momento, masa, energía y demás variables son transportadas por los torbellinos grandes. Estos además son más dependientes de la geometría y condiciones de borde del problema. Los torbellinos más pequeños tienden a ser isentrópicos e independientes de la geometría y demás aspectos macroscópicos del problema y por tanto tienen un comportamiento más predecible. Por esto las turbulencias a pequeña escala se creen más fáciles de gobernar mediante una ley universal para todas ellas que los torbellinos grandes y se modelan. La resolución de los torbellinos grandes únicamente permite utilizar mallas menos refinadas y pasos de tiempo mayores en la resolución. Esto reduce considerablemente el coste computacional frente al modelo DNS, aunque sigue siendo mucho mayor que los modelos basados en el método Reynolds-averaged Navier-Stokes (RANS). Se siguen necesitando máquinas potentes para utilizarlo, especialmente para altos números de Reynolds debido a la necesidad de resolver la energía de los remolinos con precisión.
73 Introducción a la dinámica de fluidos computacional Detached Eddy Simulation (DES). Este modelo es considerado frecuentemente como un híbrido entre los modelos RANS y LES, ya que los combina con el fin de reducir la cantidad de recursos necesarios para utilizar el modelo LES. Este modelo utiliza modelos basados en RANS, como el Spalart-Allmaras, el modelo realizable o el modelo SST para resolver las regiones del problema afectadas por capas límite. Las regiones alejadas de la influencia de capas límite, donde las escalas más grandes de torbellinos cobran más importancia, son resueltas con el método LES. Con esta combinación se consigue reducir el coste computacional respecto al método LES, pero todavia sigue siendo significativamente mayor que el necesario para modelos RANS.
74 Capítulo 4 SIMULACIONES PROPIAS DEL EFECTO MAGNUS.
75 Simulaciones propias del Efecto Magnus. 66 Para poder desarrollar la parte analítica de este proyecto es importante conocer el comportamiento de los coeficientes y propios del efecto Magnus y contar con datos fiables que determinen los valores de arrastre y sustentación en función de los parámetros que influyen en ellos. En este apartado se llevó a cabo una batería de simulaciones, con el fin de validar los datos experimentales recogidos en ELLI24, BYCH05 y WHIT04. Se pudo así determinar qué curvas escoger para elaborar el modelo matemático del funcionamiento de un aerogenerador basado en el efecto Magnus, explicado en apartados posteriores. [WHIT04] recoge valores de y obtenidos para un ensayo de un cilindro diferencial (bidimensional), mientras que [ELLI24] y [BYCH05] muestran datos adquiridos mediante ensayos en túneles de viento (tridimensional). Debido a esto se decidió realizar simulaciones de dicho efecto tanto en 2D como en 3D, con el fin de comparar cada tipo de simulación con los correspondientes datos experimentales. Para ello se utilizó el programa ANSYS , y en concreto su herramienta de modelado de dinámica de fluidos CFX Simulación bidimensional del efecto Magnus. Este primer grupo de simulaciones se llevó a cabo para determinar los valores de los coeficientes de y de manera ideal en 2D, modelando el comportamiento ideal de un fluido viscoso alrededor de un elemento diferencial de cilindro. Cabe destacar que puesto que no se tienen en cuenta las turbulencias tridimensionales ni los efectos de borde en los cilindros, se está despreciando la disipación de energía que estos fenómenos implica y por
76 Simulaciones propias del Efecto Magnus. 67 consiguiente los datos que se obtuvieron son diferentes de los coeficientes reales. Dichos valores son meramente teóricos, por lo que no se han utilizado en cálculos posteriores, aunque sí pueden considerarse como el límite máximo que pueden alcanzar los coeficientes y. Para poder comparar los resultados con los datos experimentales recogidos en ELLI24 y BYCH05 todas las simulaciones se realizaron para un , coincidiendo con los órdenes de magnitud de número Reynolds manejados en estos ensayos Geometría: Para llevar a cabo las simulaciones, se definió el dominio por el que circula el fluido como un rectángulo con un orificio que representa el cilindro. Primeramente se realizó una prueba con un dominio rectangular y el cilindro situado en el centro, observándose que el flujo tras el cilindro tardaba mucho más en estabilizarse que en la región anterior a este (Figura 35.). Por ello se decidió disminuir la región de fluido delante del cilindro y aumentar la región tras él. Figura 35. Dominio del fluido inicial.
77 Simulaciones propias del Efecto Magnus. 68 También se observó en estas primeras simulaciones que el tamaño del cilindro era demasiado grande en comparación con el dominio del fluido, lo que hacía que los límites del dominio influyeran excesivamente en la solución final del problema. Una posterior modificación fue desplazar ligeramente el centro del cilindro respecto al eje de simetría del rectángulo (Figura 36.). Con ello se buscó evitar posibles errores numéricos derivados de esta simetría. Tras estos cambios se determinaron las dimensiones finales del dominio utilizado en las simulaciones (Figura 37): Figura 36.Detalle del desplazamiento del cilindro respecto al eje de simetría del rectángulo. El diámetro del cilindro es de 0.2 m, y está separado 0.1 m del eje de simetría del rectángulo. La distancia desde la entrada de fluido en el dominio hasta el centro del cilindro es de 12.5 veces su diámetro, es decir, 2.5 m.
78 Simulaciones propias del Efecto Magnus. 69 La longitud comprendida entre el centro del cilindro hasta la salida del fluido en el dominio es de 3.5 m,17.5 veces su diámetro. Al ancho del rectángulo se le dio un valor de 4,8m, que equivale a 24 veces el diámetro. Figura 37. Geometría y dimensiones empleadas en las simulaciones bidimensionales. ANSYS CFX no permite realizar simulaciones en 2D como tal, sino partiendo de una simulación en 3D con un dominio de espesor despreciable. En el caso analizado, el rectángulo descrito anteriormente tiene un espesor de metros con el fin de conseguir que la simulación fuera en 2 dimensiones.
79 Simulaciones propias del Efecto Magnus Mallado: Para realizar el mallado se optó por utilizar celdas en forma de hexaedros a lo largo de todo el dominio (Figura 38). Esto implica un mayor gasto computacional que si se utilizaran tetraedros, pero tiene como ventajas evitar la difusión numérica ( ) y la obtención una solución más precisa. Para mejorar la distribución de las celdas a lo largo de la geometría se utilizó una malla estructurada, dividiendo la región rectangular mediante un cuadrado concéntrico con el cilindro (Figura 37). Estas dos regiones se subdividieron posteriormente con el fin de poder definir con mayor exactitud la disposición de los elementos en la malla. Esto permitió concentrar mayor número de nodos alrededor de la pared del cilindro y controlar tanto su disposición como su tamaño y forma (Figura 38).
80 Simulaciones propias del Efecto Magnus. 71 Figura 38. Detalle de la malla alrededor de la pared del cilindro. Para realizar estas simulaciones se optó por modelar la capa límite mediante funciones de forma, puesto que su cálculo supondría un coste computacional excesivo. Por esto se prestó atención a la variable y+ ( ) alrededor de la pared del cilindro, obteniendo unos valores comprendidos entre y 83.9 con una media de Esto garantiza una malla adecuada para la simulación del problema con el modelo de turbulencia escogido, el SST, que utiliza el modelo k ω estándar para resolver capas límite en las zonas cercanas a la pared y cambia a una versión del modelo k ε adaptada para altos números de Reynolds en las zonas alejadas. Para llegar al mallado final se buscó un equilibrio entre la precisión del resultado y el gasto computacional. Para ello se realizaron sucesivas mallas y aumentando el número de elementos de cada una de ellas, hasta que la diferencia en la solución fue menor del 5%.La malla final está compuesta por nodos y elementos.
81 Simulaciones propias del Efecto Magnus Condiciones de contorno: Tras obtener el mallado final se buscó aplicar al problema las condiciones de contorno apropiadas para simular en 2D las curvas recogidas en la bibliografía analizada (ELLI24 y BYCH05). Como se indicó anteriormente es necesario partir de una simulación 3D para obtener los resultados del problema simulado en 2D. Para ello se eligió un espesor muy pequeño para el dominio fluido, que permitió que en el ancho la malla sólo tuviese 2 nodos a resolver (Figura 39). Figura 39. Para conseguir que en ambos nodos las variables tuvieran igual valor se aplicó una condición de simetría a ambos lados del dominio fluido. El resto de las condiciones de contorno se establecieron buscando simular los ensayos recogidos en la bibliografía estudiada: Entrada de fluido: manteniendo constante el número de , utilizando como fluido aire a 25ºC y tomando como longitud característica el diámetro del cilindro, se determinó un velocidad de incidencia sobre el cilindro de 5. Paredes libres: se aplicó en los límites superior e inferior de la región. Inicialmente se impuso la condición de pared deslizante, lo que
82 Simulaciones propias del Efecto Magnus. 73 implicaba la eliminación de componente vertical de la velocidad del fluido no tiene componente vertical. Esta limitación resultó influir significativamente en el comportamiento del flujo, evitando que se curvase correctamente tras el cilindro. Para evitar dicho inconveniente se asignó la propiedad de pared libre, permitiendo al fluido entrar y salir libremente por dichas superficies de contorno. Salida del fluido: El aire sale por la superficie de contorno posterior de la geometría del problema, a presión atmosférica. Modelo de turbulencia: en el modelado del comportamiento del fluido se utilizó el modelo de turbulencia SST k ω, debido a su buen comportamiento en los casos que implican altos números de Reynolds, flujos turbulentos y flujos con gradiente de presión adversos, lo que lo hace apropiado para la resolución de este problema (Apéndice A3.4.1 ) Análisis de los resultados obtenidos: Tras imponer las condiciones de contorno y realizar las simulaciones se observaron las líneas de corriente y campos de presión obtenidos, con el fin de evaluar la coherencia de las simulaciones con los resultados experimentales recogidos en la biografía analizada. Se observó que el flujo tras el cilindro se curvaba correctamente, sin verse influido por ninguna condición de contorno y conforme aumentaba la rotación del cilindro el efecto aumentaba (Figura 40, a-d). Igualmente se analizó el campo de presiones alrededor del cilindro (Figura 41), para comprobar que tampoco estuviera influido negativamente por ninguna condición de contorno.
83 Simulaciones propias del Efecto Magnus. 74 a) ϕ=1. b) ϕ=2.
84 Simulaciones propias del Efecto Magnus. 75 c) ϕ=3. d) ϕ=4. Figura 40. Simulación de las líneas de corriente para diferentes valores de ϕ (velocidad de rotación periférica del cilindro/velocidad del viento incidente).
85 Simulaciones propias del Efecto Magnus. 76 a) ϕ=1. b) ϕ=3. Figura 41. Simulación de gradientes de presión para diferentes valores de ϕ (velocidad de rotación periférica del cilindro/velocidad del viento incidente).
86 Simulaciones propias del Efecto Magnus. 77 Con el fin de reproducir las curvas recopiladas en la bibliografía, se realizó una batería de simulaciones que recorriera las curvas representadas. En ésta, para cada problema resuelto, se mantuvo constante la velocidad de incidencia del viento sobre el cilindro y se varió su velocidad de rotación. En las simulaciones se recogieron, para cada caso, los valores de, y el par que es necesario aplicar al cilindro para mantener su rotación constante. En la Tabla 1 se resumen los datos más relevantes obtenidos de las simulaciones teniendo como condiciones iniciales comunes un número de , aire a 25ºC como fluido y una velocidad de incidencia de éste sobre el cilindro 5. Velocidad rotación (rad/s) Condiciones iniciales velocida d rotación (m/s) Resultados de la simulación ϕ Cl Cd Par/área (N m ) Datos de la simulación Iteracio nes Tiempo (min) ,1875 0,1769 0, ,375 0,5012 0, ,75 1,6172 0, ,125 3,1224 0, ,5 4,27 0, ,875 5,7349 0, ,25 7,0557 0, ,625 8,2094 0, ,2787 1, ,375 10,323 1, ,75 11,287 1, ,125 12,001 2, Tabla 1. Resultados de las simulaciones en 2D. Para realizar las simulaciones se utilizó un procesador Intel Centrino Duo de 1.73 GHz y 2.5Gb de memoria RAM. En la Tabla 1 se adjuntan el tiempo y número de iteraciones necesarias para que el problema convergiese.
87 Simulaciones propias del Efecto Magnus Coeficiente de sustentación ( ): Se observa que el coeficiente, representado en la Figura 42, experimenta un crecimiento con el aumento del valor de la relación de velocidades (ϕ). Para rangos en que la relación de velocidades es inferior a 1, se aprecia que deja de comportarse linealmente. Esto es coherente con lo expuesto en BYCH05, por lo que se puede deducir que el Efecto Magnus, para rangos de ϕ menores que 1 tiene un comportamiento dependiente no sólo de ϕ, sino también de otros parámetros como número de Reynolds o la rugosidad del cilindro Coeficiente de lift ϕ(relación de velocidades) Figura 42. Curva de obtenida en las simulaciones 2D. También cabe destacar que la curva que sigue el coeficiente obtenida en las simulaciones está comprendida entre la recta ideal establecida por la teoría de flujo potencial y la establecida por Bychkov como límite máximo real, teniendo una pendiente muy similar a esta última curva (Figura 43). Estos resultados entran dentro de lo esperado y respaldan la validez de las simulaciones. Al tratarse de una simulación bidimensional de un flujo real alrededor de un cilindro infinito se dan las
88 Simulaciones propias del Efecto Magnus. 79 condiciones ideales para las que el coeficiente alcance su valor real máximo. 30 Rectas de coeficiente de Lift Coeficiente de Lift Bychkov Simulaciones Teoría potencial ϕ(relación de velocidades) Figura 43. Comparación de los coeficientes máximos de sustentación obtenidos a partir de la teoría de flujo potencial, experimentalmente por Bychkov y mediante simulaciones numéricas Coeficiente de arrastre ( ): Los valores del coeficiente de drag ( ) obtenidos se muestran en la Figura 44: 2,5 Coeficiente de drag 2 1,5 1 0, ϕ(relación de velocidades) Figura 44. Curva de obtenido en las simulaciones 2D.
89 Simulaciones propias del Efecto Magnus. 80 Se puede apreciar una gran similitud con la gráfica recogida en BYCH05, salvo para los rangos de ϕ mayor de 4, que Bychkov no incluye en su estudio, y que en las simulaciones se observa un aumento drástico del. Esto refuerza la suposición de que el valor máximo de ϕ sea 4, ya que uno mayor implicaría aumentar en exceso el drag comparado con el aumento en el lift. El valor óptimo exacto de ϕ se estudiará en el Capítulo Par por unidad de superficie necesario para la rotación de cilindros: El tercer parámetro recogido en las simulaciones es el par por unidad de superficie de cilindro necesario para mantenerlo rotando a dicha velocidad (Figura 45). 25 Par (Nm) por unidad de sup ϕ(relación de velocidades) Figura 45. Par por unidad de superficie de cilindro. Simulaciones 2D. Los valores de par obtenidos son menores que los que se necesitarán realmente, puesto que las turbulencias despreciadas en las simulaciones disipan energía que es necesario aportar para mantener la rotación del cilindro. No se encontró bibliografía disponible para poder contrastar estos datos, no obstante, teniendo en cuenta lo dicho anteriormente, pueden ser un punto de partida para análisis posteriores.
90 Simulaciones propias del Efecto Magnus. 81 Tras realizar estas simulaciones se intentó simular el comportamiento del cilindro en el buje de un aerogenerador. Para ello se tomó un diferencial de cilindro y se le dio como condiciones de contorno una velocidad de rotación sobre su eje, más un desplazamiento vertical del cilindro, para simular el movimiento del cilindro alrededor del buje. Al intentar resolver dicha simulación se presentó el inconveniente de que dichas condiciones de contorno, en concreto el desplazamiento lineal del cilindro, hacían que la malla se deformase. Este caso no se puede resolver a través de FLUENT. Tras esto se estudió la posibilidad de utilizar ANSYS para resolver el problema como si se tratase de una turbomáquina, para poder simular la rotación del buje, pero debido la excesiva cantidad de recursos computacionales que resolver este caso implicaba, no se pudo abordar dicha simulación. 4.2 Simulación tridimensional del cilindro infinito. Con el objetivo de evaluar el grado de turbulencias que se desprecian en las simulaciones bidimensionales realizadas y calcular con mayor precisión la curva que describe el comportamiento de los coeficientes y, se decidió realizar una simulación tridimensional de un cilindro infinito. Un aspecto importante a tener en cuenta es el tipo de efectos que se tuvieron en cuenta en cada uno de los ensayos tridimensionales recogidos en la bibliografía. En los ensayos recogido en [ELLI24] se puede observar que el cilindro está dispuesto en el túnel de viento de modo que sus extremo sobresalen de él (Figura 46.).
91 Simulaciones propias del Efecto Magnus. 82 Figura 46.Esquema de la disposición del cilindro en el túnel de viento utilizada en los ensayos de [ELLI24]. Consecuentemente en este caso no se están recogiendo las turbulencias producidas en los bordes del cilindro. Debido a esto los valores de recogidos serán superiores a los reales, mientras que los de serán menores, al despreciar el drag inducido por las turbulencias en los bordes. Para evaluar estos datos se decidió realizar una simulación tridimensional, excluyendo en ella los efectos de borde, de la siguiente manera Geometría: El primer paso que se dio para realizar estas simulaciones fue tomar la geometría empleada en las simulaciones en 2D y aumentar el ancho de la región de fluido para, posteriormente proceder al mallado de la región.
92 Simulaciones propias del Efecto Magnus. 83 Durante el proceso de refinamiento de la malla se chocó con la limitación de que el ordenador no disponía de suficiente memoria de almacenamiento para manejar los datos necesarios en el cálculo. Con el fin de solucionar este inconveniente se tomaron varias decisiones: Se mantuvo el desfase del centro del cilindro respecto a la línea de simetría del volumen fluido, desplazándolo ligeramente hacia arriba 0.01m. Se redujo el tamaño del dominio fluido, resultando un cilindro mayor en comparación con las dimensiones del prisma por donde circula el fluido. Con esta medida se corre el riesgo que el desarrollo del flujo se vea afectado en exceso por los límites del problema, por lo que se redujo sólo lo imprescindible. Las nuevas dimensiones son: El diámetro del cilindro es de 0.1 m La distancia desde la entrada de fluido en el dominio hasta el centro del cilindro es de 8 veces su diámetro, es decir, 0.8 m. La longitud comprendida entre la salida del fluido en el dominio hasta el centro del cilindro es de 14 veces su diámetro, 1.4 m. Al ancho del rectángulo se le dio un valor de 1.6m, que equivale a 16 veces el diámetro. Con el fin de ahorrar elementos en la malla, se dividió la región del fluido en dos volúmenes (Figura 47.). Esto permitió emplear elementos hexaédricos en el volumen cercano al cilindro (color verde en la figura) para lograr una mayor precisión y tetraedros en el resto del volumen para poder realizar transiciones de tamaño entre volúmenes de elementos más rápidas, utilizando para ello menos elementos.
93 Simulaciones propias del Efecto Magnus. 84 Figura 47. Geometría empleada en las simulaciones bidimensionales Mallado: Tras el problema de memoria que surgió inicialmente se prestó especial atención a este apartado. Para poder realizar una malla con el mayor número de elementos que la memoria del ordenador permitiera se partió de una malla muy gruesa y se fue refinando hasta que volvió a superarse la memoria disponible. Durante este proceso se jugó con las posibilidades de refinar el volumen de la malla compuesta por hexaedros, y el volumen formado por tetraedros. Se observó que refinando el volumen de hexaedros consumía más memoria que haciendo lo propio en
94 Simulaciones propias del Efecto Magnus. 85 la región de los tetraedros. Como a pesar de ello una concentración de nodos grande alrededor del cilindro proporciona una solución más exacta, se optó por aumentar el número de hexaedros sacrificando el número de elementos tetraédricos especialmente en las zonas más distantes del cilindro. Se realizaron dos mallados finales, uno con mayor cantidad de elementos tetraédricos en el exterior y menor cantidad de elementos hexaédricos en el volumen que contiene al cilindro (Figura 48.), y otra malla donde la mayor concentración de nodos se daba alrededor del cilindro en el volumen compuesto por hexaedros (Figura 49.). El número de elementos resultantes para la malla primera (Figura 48.) es de 854,141, mientras que la segunda malla (Figura 49.) la forman 1,071,765 elementos. Figura 48. Malla con mayor cantidad de elementos tetraédricos en el exterior y menor cantidad de elementos hexaédricos en el volumen que contiene al cilindro.
95 Simulaciones propias del Efecto Magnus. 86 Figura 49. Malla con menor cantidad de elementos tetraédricos en el exterior y mayor cantidad de elementos hexaédricos en el volumen que contiene al cilindro En las simulaciones en 2D el principal inconveniente que se apreciaba en los mallados gruesos es que el efecto Magnus no se simulaba correctamente y se obtenían valores de y muy bajos. Idealmente para validar la malla tridimensional habría que haber ido refinándola hasta estabilizar los valores de las variables estudiadas. Ante las limitaciones encontradas en el proceso de refinamiento de la malla se tomó como criterio para elegir la malla adecuada, de acuerdo con lo observado en las simulaciones bidimensionales, tomar aquella que proporcionara mayores valores de los parámetros estudiados. Antes de validar completamente la malla se estudió el valor de la variable y+ a lo largo de la pared del cilindro. En la Figura 50 se puede apreciar que los valores de esta variable están comprendidos entre y 191.6, un rango adecuado para modelar la capa límite mediante funciones de forma.
96 Simulaciones propias del Efecto Magnus. 87 Figura Condiciones de contorno: En las condiciones de contorno de este problema cabe destacar la impuesta en las paredes laterales de la región problema. Con el fin de simular un cilindro en tres dimensiones sin efectos de borde se optó por resolver un cilindro infinito. Para ello se aplicó a ambos lados del volumen por donde circula el fluido la condición de simetría. El resto de condiciones de contorno impuestas son: La velocidad de rotación del cilindro: parámetro variable en cada simulación, cuyos valores se indican en el Capítulo 4, apartado Velocidad de entrada del fluido en la región: es el valor de la velocidad de incidencia del fluido sobre el cilindro Condición de salida del fluido de la región: el fluido sale a presión atmosférica.
97 Simulaciones propias del Efecto Magnus. 88 Condición de pared libre: los límites superior e inferior del volumen del problema permiten que el fluido entre y salga libremente por ellos facilitando así la correcta curvatura del fluido tras el cilindro. El modelo de turbulencia empleado en este caso es igualmente el modelo SST k ω, dado que el flujo que se simula en este apartado tiene las mismas características que los casos simulados en dos dimensiones Resultados: Elección de la malla adecuada: En primer lugar se llevaron a cabo las simulaciones del problema con la malla de mayor densidad de nodos en la zona formada por tetraedros (Figura 48). En esta primera aproximación se realizaron las simulaciones para un número de Reynolds de 25000, puesto que el problema converge con mayor rapidez. Las condiciones de contorno impuestas para ello fueron una velocidad del viento incidente de 3.75 para un diámetro del cilindro de 0.1m. Los resultados obtenidos se observan en la Tabla 2. Velocidad rotación (rad/s) Condiciones iniciales velocidad rotación (m/s) Resultados de la simulación Datos de la simulación ϕ Cl Cd Iteraciones Tiempo (min.) Tabla 2. Datos obtenidos en la evaluación de la primera malla. Las gráficas obtenidas de los coeficientes de y para este caso están representadas en la Figura 51 (coeficiente de sustentación) y Figura 52 (coeficiente de arrastre).
98 Simulaciones propias del Efecto Magnus Coeficiente de lift ,5 1 1,5 2 2,5 ϕ(relación de velocidades) Figura 51. Evaluación de la primera malla. Coeficiente de sustentación. 6 5 Coeficiente de drag ,5 1 1,5 2 2,5 ϕ(relación de velocidades) Figura 52. Evaluación de la primera malla. Coeficiente de arrastre. A continuación se procedió a resolver con el segundo mallado (Figura 49.) el problema para una velocidad de rotación de 120, con el objetivo de evaluar cual de las dos mayas proporcionaba los valores más altos para el coeficiente, que será la que arroje resultados más cercanos a la realidad. Los resultados obtenidos se muestran en la Tabla 3.
99 Simulaciones propias del Efecto Magnus. 90 Velocidad rotación (rad/s) Condiciones iniciales velocidad rotación (m/s) Resultados de la simulación Datos de la simulación ϕ Cl Cd Iteraciones Tiempo (min.) Tabla 3. Datos obtenidos en la evaluación de la segunda malla. A la vista de los resultados se observó que con la primera malla se obtuvo un 4.185, mientras que con la segunda malla su valor era de La segunda malla arrojaba un valor 17.52% mayor de donde se dedujo, de acuerdo a lo observado en las simulaciones bidimensionales, que era más apropiada para simular el Efecto Magnus y, por tanto, se eligió para resolver las simulaciones finales de este apartado. Aún así, debido a la falta de refinamiento de la malla los que se obtengan de las simulaciones deberán ser contrastados con los datos experimentales Resultados finales: Las simulaciones finales realizadas con la segunda malla se obtuvieron para un valor de número de Reynolds=50000, con el fin de poder compararlas con la bibliografía recogida y con las simulaciones bidimensionales realizadas previamente. Los resultados obtenidos se detallan en la Tabla 4. Velocidad rotación (rad/s) Condiciones iniciales velocidad rotación (m/s) Resultados de la simulación ϕ Cl Cd Torque (N*m/m long) Datos de la simulación Iteraciones Tiempo (min.) Tabla 4. Resultados de las simulaciones tridimensionales de un cilindro infinito.
100 Simulaciones propias del Efecto Magnus. 91 Los datos se muestran a continuación en la Figura 53, Figura 54 y Figura 55. Coeficiente de lift ϕ(relación de velocidades) Figura 53. Coeficiente de sustentación obtenido en la simulación tridimensional de un cilindro infinito. Coeficiente de drag 1,8 1,6 1,4 1,2 1 0,8 0,6 0,4 0, ϕ(relación de velocidades) Figura 54. Coeficiente de arrastre obtenido en la simulación tridimensional de un cilindro infinito.
101 Simulaciones propias del Efecto Magnus. 92 0,08 0,07 (Nm/long. cilindro 0,06 0,05 0,04 0,03 0,02 0, ϕ(relación de velocidades) Figura 55. Par por unidad de longitud obtenido en la simulación tridimensional de un cilindro infinito Análisis de los resultados obtenidos: Con el fin de analizar más detenidamente los resultados conviene observar el desarrollo del flujo obtenido en las simulaciones. En la Figura 56 se pueden observar las líneas de corriente del fluido para las condiciones generales simuladas (Re=50000 y velocidad incidente=7.5 ) y una ϕ 4. En ella se observa que no se da ningún tipo de turbulencia tridimensional a pesar de que la simulación es en tres dimensiones. En la Figura 57 se muestra, para las mismas condiciones, una perspectiva de la trayectoria que siguen unas partículas de control que se liberan en la superficie del cilindro. Igualmente aquí se observa que las turbulencias son bidimensionales, a lo largo de los ejes x e y. Esto indica que las turbulencias tridimensionales son producidas en su totalidad por las inestabilidades introducidas por las turbulencias en los bordes del cilindro. Observando los datos arrojados por estas simulaciones se pudo observar que los coeficientes analizados siguen una tendencia parecida a
102 Simulaciones propias del Efecto Magnus. 93 los resultados obtenidos en las simulaciones en 2 dimensiones. Este hecho reforzó la validez de los resultados obtenidos en las simulaciones bidimensionales, puesto que, al no darse turbulencias tridimensionales en flujo real alrededor del cilindro, el coeficiente de sustentación alcanza en este caso su valor máximo, tal y como se recoge en ellas. Figura 56. Vista lateral de las líneas de corriente obtenidas en la simulación tridimensional de un cilindro infinito. Figura 57. Trayectoria de partículas desprendidas de la pared de un cilindro infinito en rotación obtenidas para una simulación tridimensional.
103 Simulaciones propias del Efecto Magnus. 94 Se observó que el coeficiente presenta un comportamiento similar al recogido en las simulaciones en dos dimensiones tal como se aprecia en la Figura 58. Los datos arrojados en estas simulaciones, fueron similares a los obtenidos en las simulaciones bidimensionales, de lo que se deduce que a pesar de no haberse podido estabilizar las variables estimadas durante el proceso de refinamiento de la malla, la malla empleada proporciona soluciones con una exactitud razonable. Dichos valores se encuentran ligeramente por encima del valor máximo real que pueden alcanzar los coeficientes estudiados [BYCH05], por lo que dichos resultados sólo pudieron considerarse como una aproximación de la realidad y no como una solución válida. Esto es debido a que la malla utilizada en la simulación tridimensional tenía menor densidad de elementos que la malla empleada en las simulaciones bidimensionales y unos valores de la variable y+ que, a pesar de estar comprendidos en un rango óptimo para la simulación, eran mayores que los valores obtenidos en la malla de 2D. Esto hizo que la simulación tridimensional presentase mayor dificultad para capturar correctamente el Efecto Magnus. En relación a los datos recogidos en BYCH05, no pueden compararse con los obtenidos en estas simulaciones, puesto que en la bibliografía se tuvieron en cuenta turbulencias tridimensionales, que en estas simulaciones no aparecen. Por ello los valores obtenidos siguen siendo superiores a los obtenidos experimentalmente (Figura 58).
104 Simulaciones propias del Efecto Magnus COeficiente de lift Coeficiente de lift 2D Coeficiente de lift 3D Bychkov Relación de velocidades Figura 58. Comparación de los valores del coeficiente de sustentación obtenidos en las simulaciones bidimensionales y en las simulaciones tridimensionales. El comportamiento del coeficiente obtenido en estas simulaciones presenta una gran similitud con los resultados obtenidos en las simulaciones 2D (Figura 59.), aunque presenta valores ligeramente mayores. La diferencia en los valores de es debido a que la malla tridimensional es más grosera que la bidimensional, no pudiendo así modelar con la misma precisión el flujo turbulento alrededor del cilindro, responsable de la disminución de la fuerza de arrastre en régimen turbulento. Por otra parte es lógico que tengan igual comportamiento, puesto que siguen sin observarse turbulencias tridimensionales, que son las que incrementarían el valor de la parte inducida del drag. Por esta misma razón los valores obtenidos siguen siendo menores que los recogidos en [BYCH05].
105 Simulaciones propias del Efecto Magnus. 96 2,5 Coeficiente de drag 2 1,5 1 0, ϕ(relación de velocidades) Figura 59. Comparación de los coeficientes de arrastre obtenidos en las simulaciones bidimensionales (curva azul) y en las simulaciones tridimensionales (curva granate). Respecto al par que hay que proporcionar al cilindro por unidad de longitud para que este rote cabe destacar que es similar que el observado en las simulaciones 2D, excepto a altos ϕ. Esto es debido a las causas, anteriormente comentadas, por las que en esta simulación no se recoge el Efecto Magnus con la misma exactitud que en la bidimensional. El par necesario guarda una estrecha relación con el coeficiente de lift, por lo que al disminuir el coeficiente de lift obtenido resultan valores para el par. 4.3 Simulación de cilindro con efecto de borde en 3D. Con el fin de evaluar la influencia de las turbulencias debidas al efecto de borde del cilindro se decidió simular un cilindro de longitud finita. Esta simulación, además de determinar la influencia de las turbulencias tridimensionales que se derivan del efecto de borde, permitiría obtener unas curvas del comportamiento de los coeficientes y similares a los recogidos en BYCH05 (Figura 14).
106 Simulaciones propias del Efecto Magnus Geometría: En este apartado se busca simular el Efecto Magnus en un cilindro finito, cuya longitud se estimó que fuera 10 veces el diámetro, para tener una longitud lo suficientemente grande como para que el flujo se estabilice en la zona alejada del borde. Este hecho permite simular solamente medio cilindro imponiendo la condición de simetría en la sección por la que se secciona dicho cilindro. A la hora de elaborar la geometría para esta simulación se tuvieron en cuenta las limitaciones encontradas a la hora de realizar el mallado en el apartado anterior. Se optó por partir de un cilindro de igual diámetro que el elegido en el apartado y separado de los límites del fluido las mismas distancias. A la sección del cilindro se optó por darle una longitud de 0.5 m, correspondiente a la mitad del cilindro que se iba a simular, manteniendo así la relación de aspecto á 10. Posteriormente se añadió un volumen de fluido en el borde del cilindro de 0.2m de espesor (dos veces el diámetro del cilindro) para poder observar las turbulencias en el borde del cilindro (volumen de color verde en la Figura 60.). Figura 60. Geometría empleada en la simulación tridimensional de un cilindro finito.
107 Simulaciones propias del Efecto Magnus. 98 El volumen donde se encuentra alojado el cilindro se subdividió en dos partes, permitiendo así diferenciar entre una región más cercana al cilindro en la que se refinó más la malla y otra exterior en que el mallado hizo más grueso (Figura 60). A su vez al volumen que no contiene al cilindro (de color verde en la Figura 60) se le realizó un corte cilíndrico, coincidiendo con el cilindro (cilindro marrón en Figura 61). El objetivo de esta acción era crear la superficie correspondiente a la sección lateral del final del cilindro (sección verde en la Figura 62), para así poder imponer en ella la condición de rotación requerida en la simulación. Figura 61. Figura 62. Esta división también permitió refinar la malla en la región cercana al borde del cilindro y capturar así las turbulencias con mayor precisión Mallado: Para realizar el mallado se decidió emplear elementos hexaédricos en el volumen que contiene al cilindro y la región cilíndrica adyacente al borde de éste, y tetraedros para los volúmenes exteriores más distantes
108 Simulaciones propias del Efecto Magnus. 99 Teniendo presente el inconveniente surgido durante el proceso de mallado en la primera simulación tridimensional, se buscó obtener la malla más precisa que las limitaciones de memoria del ordenador permitiesen. Anteriormente se observó que el refinado de los elementos hexaédricos consumía más memoria que el refinado de los tetraedros. Por esto primero se procedió a refinar lo más posible la región compuesta por hexaedros hasta alcanzar el límite permitido por la memoria del ordenador y, posteriormente refinar las regiones formadas por tetraedros. La malla resultante tiene 921,365 elementos, su aspecto puede verse en la Figura 63. Figura 63. Malla empleada en la simulación tridimensional de un cilindro finito en rotación. Finalmente se comprobaron los valores que la variable y+ adquirió a lo largo de la pared del cilindro, para comprobar que la capa límite fue modelada correctamente. En la Figura 64 se puede observar que los
109 Simulaciones propias del Efecto Magnus. 100 valores están comprendidos entre 26.3 y 101.5, lo que supone un rango válido para modelar la capa límite mediante funciones de forma. Figura Condiciones de contorno: Las simulaciones realizadas se han llevado a cabo para un menor número de Reynolds=250000, con el fin de facilitar la convergencia de la simulación. Las condiciones de contorno impuestas son: Velocidad de incidencia del fluido de Velocidad de rotación en el cilindro de 120. Salida del fluido del volumen simulado por la superficie de contorno posterior a presión atmosférica. Condición de simetría en el extremo del volumen en que el cilindro se encuentra seccionado.
110 Simulaciones propias del Efecto Magnus. 101 Condición de pared libre en las superficies superior, inferior y lateral del volumen, que permite el paso libre de fluido, para evitar interferencias de las condiciones de contorno en el correcto desarrollo del flujo Resultados obtenidos: Tras resolver esta simulación para las condiciones anteriormente indicadas y una velocidad de rotación de 120, se obtuvieron los siguientes resultados: Velocidad rotación (rad/s) Condiciones iniciales velocidad rotación (m/s) Resultados de la simulación Datos de la simulación ϕ Cl Cd Iteraciones Tiempo (min.) Se puede apreciar que el coeficiente es ligeramente superior que el obtenido en las simulaciones anteriores para un ϕ 1.6 ( en las anteriores simulaciones en tres dimensiones para el mismo Re). Esto se debe a que las turbulencias introducidas debido a los efectos de borde implican un incremento en la componente inducida del drag). Por otra parte el coeficiente muestra un valor muy inferior al obtenido en las simulaciones anteriores y los recogidos en la bibliografía, tan sólo es aproximadamente dos veces mayor que el. Esto, como ya se ha comentado, es debido a la falta de densidad de nodos en la malla que impide simular correctamente el Efecto Magnus. Sin embargo, la observación del flujo que se obtiene arrojó conclusiones interesantes.
111 Simulaciones propias del Efecto Magnus. 102 Figura 65. Simulación de un flujo alrededor de un cilindro con efecto de borde. Figura 66. Vista superior de la simulación de un flujo alrededor de un cilindro con efecto de borde. En las Figura 65 y Figura 66 se puede observar la trayectoria que siguen las partículas de control que parten desde la pared del cilindro en rotación. En esta imagen se aprecia que las turbulencias causadas en el borde del cilindro tienen una gran influencia en el resto del flujo que rodea al cilindro y por tanto es un efecto que no se puede despreciar.
112 Simulaciones propias del Efecto Magnus. 103 Figura 67. Vista posterior en detalle de la estela creada tras el cilindro. En la Figura 67 se aprecia cómo el flujo de aire forma un vórtice, similar al que se forman en los bordes de las alas de los aviones debido a la diferencia de presión entre la parte superior e inferior. Cabe esperar por tanto que en una simulación más precisa, en la que el Efecto Magnus se represente más fielmente y el lift sea mayor, las diferencias de presión entre la parte superior e inferior del cilindro sean mayores, y por tanto estos vórtices sean mayores. En estudios posteriores sería interesante realizar un estudio pormenorizado sobre las turbulencias que se producen en los bordes del cilindro y el grado de influencia que estas tienen en los valores de y. Una vez conocido el funcionamiento de este fenómeno se podría estudiar posible soluciones para minimizar el incremento en el debido a estas turbulencias, como por ejemplo, la colocación de discos en los bordes del cilindro para separar el flujo alrededor del cilindro y el resto de fluido (BYCH05).

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución