• Author / Uploaded
• Diego Rodríguez

• Categories
• Integral
• Viscosidad
• Unidad Central de procesamiento
• Ecuaciones de Navier-Stokes
• Ecuaciones

Citation preview

SIMULACIÓN DEL FLUJO EN EL INTERIOR DE UN TÚNEL DE VIENTO CON OPENFOAM Diego Rodríguez Rodríguez 21 de diciembre de 2018 Profesora: Elena B. Martín

Mecánica de Fluidos II y CFD1

1

Del inglés Computational Fluid Dynamics, ‘dinámica de fluidos computacional’, parte de la asignatura en la que se engloba este proyecto.

1

2

Índice Índice

3

Índice de figuras

4

Índice de cuadros

4

1. Objetivo 1.1. Definiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 6

2. Datos para el análisis 2.1. Caso Re = 10 . . . 2.2. Caso Re = 100 . . . 2.3. Caso Re = 300 . . . 2.4. Caso Re = 1000 . . 2.5. Caso Re = 200000 .

. . . . .

11 11 11 11 11 12

. . . . . . . . . .

12 12 12 13 13 14 14 16 16 16 16

4. Condiciones de contorno 4.1. Generalidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Variables turbulentas y funciones de pared . . . . . . . . . . . . . . . . . . .

16 16 17

5. Sensibilidad de dominio computacional y malla 5.1. Dominio computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Malla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18 18 19

6. Resultados 6.1. Caso Re = 10 . . . . . . . 6.2. Caso Re = 100 . . . . . . . 6.3. Caso Re = 300 . . . . . . . 6.4. Caso Re = 1000 . . . . . . 6.5. Casos Re = 200 mil . . . . 6.5.1. Modelo K − ω . . 6.5.2. Modelo K − ω SST

22 22 25 30 34 38 38 41

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

3. Entorno de OpenFOAM 3.1. Carpeta system . . . . . . . . 3.1.1. blockMeshDict . . . . 3.1.2. controlDict . . . . . 3.1.3. fvSchemes . . . . . . . 3.1.4. fvSolution . . . . . . 3.1.5. Archivos adicionales . 3.2. Carpeta constant . . . . . . 3.2.1. transportProperties 3.2.2. transportProperties 3.3. Carpeta 0 . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

Anexo A. Código gnuplot medias

46

Anexo B. Código FFT MATLAB

47

Referencias

48

3

Índice de figuras 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37.

Niveles de refinamiento. . . . . . . . . . . . . . . Esquema dominio. . . . . . . . . . . . . . . . . . Residuos Re 10. . . . . . . . . . . . . . . . . . . . Residuos finales Re 10. . . . . . . . . . . . . . . . Coeficientes Re 10. . . . . . . . . . . . . . . . . . Coeficientes ampliados Re 10. . . . . . . . . . . . Animación flujo Re 10. . . . . . . . . . . . . . . . Residuos Re 100. . . . . . . . . . . . . . . . . . . Residuos finales Re 100. . . . . . . . . . . . . . . Coeficientes Re 100. . . . . . . . . . . . . . . . . Señal CL Re 100. . . . . . . . . . . . . . . . . . . Señal CL transformada Re 100. . . . . . . . . . . Animación flujo Re 100. . . . . . . . . . . . . . . Residuos Re 300. . . . . . . . . . . . . . . . . . . Residuos finales Re 300. . . . . . . . . . . . . . . Coeficientes Re 300. . . . . . . . . . . . . . . . . Señal CL Re 300. . . . . . . . . . . . . . . . . . . Señal CL transformada Re 300. . . . . . . . . . . Animación flujo Re 300. . . . . . . . . . . . . . . Residuos Re 1000. . . . . . . . . . . . . . . . . . Residuos finales Re 1000. . . . . . . . . . . . . . . Coeficientes Re 1000. . . . . . . . . . . . . . . . . Señal CL Re 1000. . . . . . . . . . . . . . . . . . Señal CL transformada Re 1000. . . . . . . . . . . Animación flujo Re 1000. . . . . . . . . . . . . . . Residuos Re 200 mil (K − ω). . . . . . . . . . . . Residuos finales Re 200 mil (K − ω). . . . . . . . Coeficientes Re 200 mil (K − ω). . . . . . . . . . Coeficientes ampliados Re 200 mil (K − ω). . . . Animación flujo Re 200 mil (K − ω). . . . . . . . Residuos Re 200 mil (K − ω SST). . . . . . . . . Residuos finales Re 200 mil (K − ω SST). . . . . Coeficientes Re 200 mil (K − ω SST). . . . . . . Señal CL Re 200 mil (K − ω SST). . . . . . . . . Señal CL transformada Re 200 mil (K − ω SST). Animación flujo Re 200 mil (K − ω SST). . . . . Animación K Re 200 mil (K − ω SST). . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15 18 22 23 23 24 24 25 26 27 28 28 29 30 31 31 32 33 33 34 35 35 36 37 37 39 39 40 40 41 42 42 43 43 44 44 45

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

13 17 18 19 19 20 20 20 20

Índice de cuadros 1. 2. 3. 4. 5. 6. 7. 8. 9.

Funciones. . . . . . . . . . . . . . . . . . Lista cond. contorno. . . . . . . . . . . . Funciones de pared. . . . . . . . . . . . . Valores medios de CD dominios. . . . . . Variación % CD . . . . . . . . . . . . . . . Valores medios de CD mallas. . . . . . . % sin refinementBox. . . . . . . . . . . % refinementBox 1, refinementSurface % refinementBox 1, refinementSurface 4

. . . . . . . . . . . . . . 2. 3.

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

10. 11.

% refinementSurface 2/3. . . . . . . . . . . . . . . . . . . . . . . . . . . . Valores medios de y + . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

21 38

NOTA PRELIMINAR El enlace con los archivos de la simulación se puede encontrar aquí debido a su tamaño. Además, se recomienda mantener las animaciones en la misma carpeta del informe para una visualización correcta desde el mismo.

1. Objetivo Calcular el flujo alrededor de un prisma cuadrado de sección 0.06 m x 0.06 m situado en el interior de un túnel de viento. La sección transversal inicial de trabajo del túnel de viento es de 0.6 m x 0.6 m y su longitud de 1.25 m, con el prisma en el centro de la misma. Para simplificar, la simulación está hecha en 2 dimensiones. Obtener el coeficiente de resistencia (o arrastre) y sustentación del prisma bajo estas condiciones de operación del túnel: Re = 10, 100, 300, 1000, 200000. Realizar un análisis de sensibilidad de los resultados tanto a la malla como al dominio computacional escogido. Comparar los resultados de los modelos K −ω y K −ω SST de las simulaciones con flujo turbulento. Hallar el número de Strouhal en los casos con calle de Kármán.

1.1. Definiciones [1] Número de Reynolds: Número adimensional principal en la simulación que caracteriza el movimiento de un fluido relacionando las fuerzas inerciales con las viscosas. Re =

ρv0 lc v0 l c = , µ ν

(1)

donde v0 es la velocidad del fluido incidente, lc la longitud característica, ν la viscosidad cinemática. Número de Strouhal: Número adimensional que relaciona la velocidad característica de oscilación con la velocidad media del fluido. St =

f lc , v0

(2)

donde f es la frecuencia de oscilación del fluido. Coeficiente de resistencia: CD =

FD , 1 2 ρSv 0 2

(3)

donde FD es la fuerza de resistencia, S la superficie donde se mide el coeficiente. Coeficiente de sustentación: CL =

FL , 1 2 ρSv 0 2

(4)

donde FL es la fuerza de sustentación, S la superficie donde se mide el coeficiente. Calle de Kármán: Torbellinos que aparecen tras un objeto inmerso en un fluido cuando, a una velocidad suficientemente alta, la capa límite se desprende del cuerpo formando un patrón periódico. Esto ocurrirá normalmente a Re > 100, y es relacionado con el número de Strouhal. 6

Dominio computacional: Modelo de la realidad donde se va a simular estableciendo unas ecuaciones de gobierno. Malla: Discretización del dominio computacional para aplicar el método de volúmenes finitos con las ecuaciones de gobierno las de Navier-Stokes. Hay diferentes métodos de mallado, estructurado, con celdas de geometría uniforme que demandan menos recursos pero son difíciles de adaptar a dominios complejos; y no estructurado, con celdas más irregulares que siendo menos eficientes en recursos computacionales se adaptan mejor a dominios de forma complicada. Método de Volúmenes Finitos: Aplicación de las ecuaciones de gobierno al dominio discretizado, con cada celda como volumen de control. OpenFOAM hace una discretización temporal implícita de las ecuaciones de Navier-Stokes en forma integral, Ecuación de conservación del momento lineal: ∫ I d ρ dV + ρ(⃗v − ⃗vC )⃗η dS = 0 dt VC (t)

(5)

SC (t)

Ecuación de conservación del momento lineal: ∫ ∫ I I d τ · ⃗η dS + ρf⃗m dV ρ⃗v dV + ρ⃗v (⃗v − ⃗vC )⃗η dS = dt VC (t)

SC (t)

SC (t)

(6)

VC (t)

Ecuación de conservación de la energía interna: ∫ I I ∫ I d ⃗ ·⃗η dS + (Q˙ r + Q˙ q )dV (7) ρ · e dV + ρ · e(⃗v −⃗vC )⃗η dS = ⃗v (τ ·⃗η )dS + k ∇T dt VC (t)

SC (t)

SC (t)

SC (t)

VC (t)

que se pueden representar en notación compacta para ⃗vC = 0 (volumen de control cte.) como ∫ ∫ ∫ ∫ ∂ρϕ ⃗ · ⃗η dS + Sϕ dV, dV + ρϕ⃗v⃗η dS = Γϕ ∇ϕ (8) ∂t VC

SC

SC

VC

con ϕ, Γϕ , Sϕ como variables escalares o vectoriales en cada ecuación. Ahora evalúa temporalmente el residuo de la ecuación (8) –sin el término no estacionario– de manera implícita (en t + ∆t) nf ∑ f =1

ϕf (t + ∆t)F =

nf ∑

(Γϕ )f (∆ϕ(t + ∆t))f S + VP Qu + VP QP ϕP (t + ∆t),

(9)

f =1

con los subíndices f de cara de la celda, y P del centroide de la celda. Tras evaluar también el residuo de la parte no estacionaria (método Euler implícito, por ejemplo) y discretizar los términos ϕf (t + ∆t) y ∆ϕf (t + ∆t), estos serán dependientes de ϕP (t + ∆t), quedando una ecuación matricial para la malla del estilo AΦ(t + ∆t) = b(t), siendo A la matriz, de orden n2P .

7

(10)

Número de Courant: Número adimensional muy importante en las simulaciones que determina el paso de tiempo ∆t a considerar respecto de la velocidad y el tamaño de celda. v0 ∆t 0 Subcapa logarítmica (30 < y + < 300): Zona de la capa límite turbulenta más alejada de la pared, donde los esfuerzos turbulentos y viscosos comparten gran importancia. Aquí se define u+ como u+ = e

1 ln y + + C , κ

uτ y , ν ( ) ∂⟨u⟩ uτ ∼ κy . ∂y y+ =

siendo ahora

(28) (29)

(30)

La constante C y la de von Kármán κ están en los siguientes rangos de valores según Hoffmann [5] dependiendo de la rugosidad de las paredes: κ ∈ (0.38, 0.43) , C ∈ (5, 5.5).

2. Datos para el análisis Los mismos para todos los casos: Viscosidad cinemática del fluido newtoniano Altura del túnel (cte.) Longitud del túnel (véase fig. 2) Lado del prisma cuadrangular Longitud característica (diagonal del cuadrado)

ν = 1·10−5 m2 s−1 0.6 m 2.5 m 0.06 m lc ≃ 0.08485 m

2.1. Caso Re = 10 Velocidad incidente dirección OX v0 ≡ u0 ≃ 0.00118 m s−1

2.2. Caso Re = 100 Velocidad incidente dirección OX v0 ≃ 0.0118 m s−1

2.3. Caso Re = 300 Velocidad incidente dirección OX v0 ≡ u0 ≃ 0.0354 m s−1

2.4. Caso Re = 1000 Velocidad incidente dirección OX v0 ≡ u0 ≃ 0.118 m s−1

11

2.5. Caso Re = 200000 (Valores iniciales, no aplicables necesariamente a todos los contornos, ver subsección 4.2.)

Velocidad incidente dirección OX Viscosidad cinemática turbulenta Energía cinética turbulenta Tasa de disipación de energía turbulenta

v0 ≡ u0 ≃ 23.57 m s−1 νT = 0 m² s−1 K ≃ 2.0833 m² s−2 ω ≃ 34.366 s−1

3. Entorno de OpenFOAM 6 El funcionamiento del programa se basa en código en terminal, que llama a archivos que crean el dominio computacional y la malla, la refinan, recortan una figura concreta en el dominio, simulan, visualizan los resultados, etc. Estos archivos se encuentran en diferentes carpetas del caso a simular, siendo las tres imprescindibles, system, constant, y la carpeta que contiene las condiciones iniciales y de contorno de nuestro problema, que tendrá de nombre el tiempo donde comience la simulación. Se detallan a continuación su contenido e importancia a fin de aclarar terminología y conceptos que se usaren en siguientes secciones.

3.1. Carpeta system Aquí se encuentran archivos relacionados directamente con la simulación, y a los que se llama o usa para realizarla. Hay 4 archivos imprescindible para comenzar un análisis, blockMeshDict, controlDict, fvSolutions, y fvSchemes. 3.1.1. blockMeshDict Con este archivo de código se construye nuestro dominio computacional y se define su discretización indicando el número de celdas en cada eje. Para crear el dominio, en nuestro caso se da valores a 6 variables, que indican la coordenada menor y mayor en cada eje –con el grosor en OZ ≪ que en OX y OY para simular en 2D (véase apartado 3.1.5)–. A continuación pasamos a definir los vértices, en este caso que compondrían la geometría de un túnel de viento cuadrangular, con las coordenadas anteriormente determinadas de manera que se indica la posición de los 4 vértices de cada base del hexaedro que sería el dominio. El apartado blocks fija una numeración a los vértices del hexaedro (hex) que se usarán para definir las caras del contorno; se ajusta también el número de celdas en cada dirección, y si se quiere que estas adquieran un tamaño progresivo en alguna dirección (simpleGrading), que en este proyecto no usaríamos esta funcionalidad. El aparatado edges queda vacío dado que no nos es necesario crear geometrías curvas. En boundary es donde se dan nombre a las caras del hexaedro y el tipo de contorno que serían. En este problema están definidos 6 contornos, 4 como patch y 2 como empty: Los 4 primeros tienen el tipo general patch que nos permitiría aplicar unas condiciones de contorno en esas caras, y los 2 empty indica al programa que esas caras, perpendiculares al 3er eje (OZ) y como la simulación es en 2D se obviaría, serían irrelevantes en la simulación no siendo necesario resolver ahí nuestro problema, quedan “vacías”. Los nombre a las caras son de suma importancia a la hora de ajustar correctamente las condiciones de contorno, y se indican encima de cada contorno antes del entorno entre llaves {}. En nuestro caso contamos con los contornos patch (inlet, oulet, top, bottom), y empty (front, back). Por último, la numeración de los cuatro vértices conforma la cara en el contexto faces. 6

Open Field Operation And Manipulation.

12

3.1.2. controlDict Este archivo controla los parámetros temporales de la simulación así como el intervalo de escritura de los resultados y el modo. Comienza especificando el solver a utilizar, en nuestro caso icoFoam o pisoFoam, y el tiempo inicial, que no siempre tiene por qué ser 0 ya que se puede comenzar una simulación en el tiempo deseado con tal de tener la carpeta con esas condiciones iniciales. A continuación se indica el tiempo final en el que debe concluir la simulación, y el paso de tiempo entre iteraciones, muy importante para nuestros solvers dado que en flujo no estacionario dependemos del número de Courant (visto en 1.1) y ∆t es algo a definir cautelosamente para que no nos explote la simulación. Las dos siguientes líneas controlan la escritura de resultados, con writeControl se indica el intervalo de tiempo que se escribirán resultados, siendo los comandos más comunes timeStep, escribiendo datos cada deltaT, y adjustableRunTime para que lo haga cada segundo de tiempo de simulación ajustando los pasos de tiempo si fuera necesario para que coincidan. En la línea debajo podemos definir cada cuánto nos muestre los valores calculados en el directorio principal del caso como carpetas del intervalo de tiempo escogido, una cada uno. Para que además de calcular los campos de presiones y velocidades del problema, el entorno functions sirve para incluir funciones definidas en OpenFOAM y que habrá que configurar con archivos específicos para cada función a llamar. Los datos se guardarán en una carpeta postProcessing que se creará al efecto. Comando función residuals probes forceCoeffs forcesIncompressible

Definición Recoge los residuos mayores de las soluciones para cada componente de los campos definidos en su archivo específico. Saca los datos de las variables escogidas en las celdas más cercanas a un punto especificado. Calcula los coeficientes aerodinámicos (CD y CL ) además de un tercero de momento7 , sobre contorno/s definido/s. Calcula la presión y los esfuerzos viscosos en uno o varios contornos seleccionados para casos incompresibles.

Cuadro 1: Funciones comunes y su uso. (Para mejor entendimiento ver 3.1.5)

3.1.3. fvSchemes Este archivo define las operaciones y el modo de discretización que el solver usaría, los modelos de cálculo de los operadores de las ecuaciones de gobierno discretizadas. En ddtSchemes se indica el tipo de discretización del término no estacionario de la ecuación 8, eligiendo ⃗ , ∇· ⃗ , ∇ ⃗ 2 se determiel método de Euler implícito. La discretización de los operadores ∇ na en los siguientes tres entornos (gradSchemes, divShemes, laplacianSchemes) con ϕ la variable que acompaña al operador (p, U, k, omega...) y el método de cálculo Gauss-linear que es el modo estándar de discretización del MVF para integración de Gauss (teorema de Gauss-Ostrogradsky), aun habiendo más que se podrán escoger en base a ese. 7

Coeficiente de momento: CM =

M

, 1 2 2 ρv0 S ·c

(31)

donde las nuevas variables serían M como el momento que provocan la fuerzas aerodinámicas (cabeceo en un avión), y c que podría extrapolarse al perímetro del prisma en este proyecto (normalmente es la longitud de la cuerda de un ala en aeronáutica).

13

3.1.4. fvSolution Aquí se definen los solvers y tolerancias que se usará para la simulación. Hay que especificar para cada variable el tipo de solver, el preacondicionador de matriz –para solvers de gradiente conjugado–, tolerancia de soluciones y la relativa del solver. Estos dos últimos parámetros son muy importantes ya que configuran el orden de magnitud de los residuos de las iteraciones internas –para resolver la matriz inversa– y las soluciones de la ecuación (32), indicando primero en tolerance el valor del residuo inicial menor con el cual se consideraría que la simulación habría convergido y pararía automáticamente. En relTol se fija el ratio entre el residuo inicial y final en cada iteración, cuanto más se acerca a 0 más próximos estarán, obligando con 0 a obtener el valor de tolerance en cada ∆t. Cuando se use un solver con suavizante para reducir las iteraciones de la resolución del sistema de ecuaciones (smoothSolver) se deberá indicar el modo, con preferencia por el Gauss-Seidel y el simétrico (GaussSeidel y symGaussSeidel). Para el caso de usar código PISO, SIMPLE, o PIMPLE, relacionados con la resolución del campo de presiones de manera iterativa, se deben indicar unos correctores como las veces que se resuelve el caso y la ecuación de la presión8 en cada paso de tiempo para actualizar resultados. Adicionalmente también se indica una presión de referencia en las celdas (pRefCell y pRefValue) para que el problema dé un rango de presiones relativas y no la absoluta (caso incompresible). 3.1.5. Archivos adicionales Incluyendo otros archivos en la carpeta nos abre el abanico a más funcionalidades. En este proyecto se hace uso primordial del archivo snappyHexMeshDict y extrudeMeshDict para el conformado final del dominio computacional y el refinado de la malla; además también se usa la ejecución en paralelo por lo que el archivo decomposeParDict será con lo que controlemos esa función, y luego están los códigos para la configuración del cálculo de las funciones indicadas en controlDict (3.1.2). snappyHexMeshDict: Dada la profundidad de esta inestimable función se describirán solo los apartados que a nuestro caso les hayan sido necesarios. Para modelar el dominio computacional de este caso, se extruyó el prima central del túnel en Solidworks exportándolo en formato .stl (estereolitografía). Con ello en la carpeta triSurfaces (ver sección 3.2), se especifica el nombre del archivo en el entorno geometry junto con el tipo de geometría y su nombre que será el del contorno. Añadiendo más {} definimos más geometrías, y como en nuestro dominio se decidió realizar un refinado a una zona de la malla (ver subsección 5.2), usamos el código implícito searchableBox para acotar una caja con las coordenadas del vértices mínimo y máximo como está en la fig. 2. Dentro del apartado castellatedMeshControls, en refinementSurfaces, se define las superficies del dominio que se quiera tener con un tamaño de celda menor al inicial del archivo blockMeshDict. Aquí tenemos varios niveles de refinamiento a escoger, siendo 0 el tamaño inicial de celda (no refina) y no tiene valor máximo (véase fig. 1). El primer número del paréntesis significa el nivel de refinamiento en esas celdas del tipo de contorno especificado en la línea siguiente, mientras que el segundo aplica ese nivel –si fuera distinto– a celdas que tengan mayor ángulo en su intersección entre ellas que el descrito en el comando resolveFeatureAngle posterior al entorno refinementSurfaces. 8

Ecuación Poisson para la presión: ⃗ 2 P = −∇·(ρ⃗ ⃗ ⃗ v ). ∇ v ∇⃗

14

(32)

En refinementRegions, si hemos especificado previamente en el apartado geometry una geometría distinta del .stl, podremos elegir el nivel y lugar de refinamiento respecto a la geometría. Existen 3 modos: inside refina las celdas del interior de la geometría con el nivel tomado del segundo número del paréntesis, outside es igual pero refinando el exterior, y distance refina de acuerdo a una distancia de la geometría indicada en el primer valor del paréntesis con un nivel del segundo valor. El primer número es ignorado en los otros 2 modos. El comando locationInMesh da una coordenada que es importante que esté fuera de la figura cerrada a recortar y refinar, dado que indica la región que mantendrá las celdas inalteradas. Los demás entornos definen el recorte (snap) y demás ajustes avanzados de añadidura de celdas y alisado y otros ajustes avanzados.

Figura 1: Niveles de refinamiento de las celdas. Extraído de simscale.com. extrudeMeshDict: Este código brinda la funcionalidad de, usando la malla resultado de las operaciones del snappyHexMesh y dado que nuestro caso se simulará en 2D, convertir el dominio a un estado bidimensional que cumple ahora las condiciones del tipo de contorno empty definidas en las paredes frontal y trasera del dominio en el archivo blockMeshDict. decomposeParDict: Con este archivo, en las simulaciones corren en paralelo, configuramos el peso de la simulación que cada procesador ejecutará y el método. En numberOffSubdomains indicamos la cantidad de descomposiciones del dominio que se quieran hacer, que puede o no coincidir con el numero de procesadores/hilos de la CPU. Esto creará tantas carpetas como subdominios e irá escribiendo los resultados temporales (copia las condiciones iniciales en cada una) de sus parte descompuesta que luego habrá que unir con el comando en terminal reconstructPar. La línea preserveFaceZones es opcional y lo que hace es mantener en un mismo procesador celdas completas y vecinas, no descomponer mitad de una. En cuanto a los métodos, el que se ha usado ha sido simple, que mediante unos coeficientes ajustamos la carga de celdas que soportan los procesadores en cada dirección como el número de descomposiciones por eje, en nuestro caso (3 2 1) en (OX OY OZ). Para una carga igual se usaría el método scotch, por ejemplo. probes: Aquí configuramos la función del mismo nombre llamada en controlDict, asignando un nombre a la carpeta donde se guardarán los datos, la frecuencia de escritura y su intervalo, en este caso timeStep y 1, escribe cada ∆t. Por último hay que especificar los campo que quieras que calcule y dónde, asignando uno o varios puntos con su coordenada.

15

forceCoeffs: Al igual que probes, es un archivo de configuración en este caso de la función que calcula los coeficientes de las fuerzas. El control del tiempo de escritura es el mismo que el anterior, y se indica el contorno donde se quiera calcular los coeficientes, la incompresibilidad o no del fluido, la dirección de FL y FD , así como el módulo de la velocidad incidente, lc , y el área de referencia (grosor del dominio en OZ·altura figura).

3.2. Carpeta constant En este directorio se almacena el resultado del dominio y su discretización y el resultado de las operaciones de snappyHexMesh en la carpeta polyMesh como archivos con las posiciones de las celdas, sus caras y las vecinas, los vértices, y los contornos. Además, como se vio en el apartado 3.1.5, para que snappyHexMesh recorte una figura ha que almacenar el archivo en formato .stl en la carpeta triSurface, de donde lo cogerá y usará en sus operaciones. 3.2.1. transportProperties Archivo breve pero imprescindible, pues define el valor de la viscosidad cinemática ν del fluido de trabajo, entre corchetes9 las dimensiones (m2 s−1 ) seguidas de su valor numérico. 3.2.2. transportProperties Solo necesario para simulaciones en régimen turbulento, configura el modelo computacional a usar. Primero se especifica el tipo de modelo, RAS o LES, y dentro de este el modelo según el orden (visto en mayor profundidad en 1.1), en nuestro caso kOmega y kOmegaSST. Dejando turbulence on, y printCoeffs on, pudiendo definir los coeficientes de cierre distintos de los implícitos con el entorno kOmegaCoeffs {} o kOmegaSSTCoeffs {} (de igual modo para otros modelos).

3.3. Carpeta 0 (cual sea el tiempo inicial) Esta carpeta, que tendrá por nombre el valor numérico del tiempo en el que inicie la simulación, contiene los archivos de las variables de las ecuaciones con las condiciones iniciales y de contorno (detalladas en la sección 4), campo de velocidades y presiones en laminar/transitorio, y demás variables turbulentas del modelo a usar para condiciones a ese régimen. En cuanto a las condiciones iniciales, en nuestro problema se han dado valores para los términos de la ecuación (1), en los que se despeja v0 y el código la calcula previa indicación en la línea U #calc "$Re*$ni/$D" (v0 = Relc· ν ). Previamente se han debido definir las dimensiones de la variable fluida que estemos configurando, así como el valor que tendría en en el dominio computacional inicialmente, en cada celda antes que sean afectadas por el fluido de trabajo (en nuestro caso al túnel no llega aire hasta que no comienza la simulación, quedan en 0 la/s componente/s).

4. Condiciones de contorno 4.1. Generalidades Imprescindibles para la resolución de las ecuaciones de gobierno, estas condiciones restringen el comportamiento del fluido a lo largo del dominio computacional definiendo unas condiciones en los contornos que atraviesa. Estas, que están resumidas en la tabla 2, se 9

Dimensiones: [kg m s K mol A cd].

16

indican para cada variable fluida en sendos archivos en la carpeta constant en el entorno boundaryField con el nombre del contorno encabezando cada tipo de condición. Comando condición

Definición

noSlip

Condición de no deslizamiento, la velocidad en dichas paredes sería 0.

fixedValued

Se especifica un valor fijo de la variable en ese contorno con la línea value y el escalar/vector uniforme o no.

zeroGradient

El gradiente de la variable es 0 en esa pared

∂ϕ ∂xi

= 0.

empty

No se impone condición al ser un contorno vacío.

symmetryPlane

Indica que el contorno actúa como plano de simetría.

calculated

Asume un valor calculado con las demás de la carpeta 0 y/o con el modelo de turbulencia implícito; es opcional especificar un valor.

Cuadro 2: Tipos de condiciones de contorno más usadas y su interpretación.

4.2. Variables turbulentas y funciones de pared Para cuando el fluido de trabajo se encuentra en régimen turbulento hay unas condiciones de contorno específicas que prevalecen para las variables fluidas del modelo de turbulencia a usar. Estas son las leyes de pared, que definen el comportamiento del fluido en la zona más próxima a las paredes (ver apartado 1.1). La resolución de las ecuaciones de los modelos de turbulencia en la subcapa viscosa con la condición de contorno de no deslizamiento dan resultados poco fiables debido a su incapacidad de predecir el valor de C . Para ello hay que imponer unas condiciones de contorno específicas llamadas funciones de pared para las variables a resolver en cada modelo de turbulencia, en nuestro caso K y ω. Para determinar sus condiciones iniciales se opera de la siguiente manera para cada variable: 3 3 K = (⟨u⟩I)2 = (0.05·23.57 m s−1 )2 ≃ 2.0833 m² s−2 , 2 2

(33)

con I la escala turbulenta característica para cada túnel pero que aquí se ha supuesto del 5 %. Una vez obtenida la energía cinética turbulenta, esta se puede relacionar con la tasa de disipación de esta manera: √ √ K 2.0833 m² s−2 = ≃ 34.366 s−1 , (34) ω= ℓ 0.07·0.6 m siendo ℓ la escale que describe el tamaño de los torbellinos con gran carga energética en un flujo turbulento, que se estima un 7 % del diámetro del túnel de viento (altura en nuestro caso). En el caso de νT , su valor inicial lo definiremos 0 m² s−1 .

17

Comando función kqRWallFunction omegaWallFunction nutkWallFunction

Interpretación Condición de contorno para K como cubierta sobre la de zeroGradient para implementar las funciones de pared con alto Re. Condición en la pared para la tasa de disipación, √ siendo mezcla de la 2 2 subcapa viscosa y logarítmica, ω = ωvisc. + ωlog. . Esta condición para νT usando funciones de pared se basa en la energía cinética turbulenta.

Cuadro 3: Funciones de pared como condiciones de contorno para las variables turbulentas.

5. Sensibilidad de dominio computacional y malla 5.1. Dominio computacional Para estimar la longitud adecuada del túnel de viento escogí una malla de 100 x 50 (OX OY) celdas iniciales con la superficie del cuadrado con un nivel de refinamiento 3, y una caja de refinamiento interior 1, saliendo un total de 42994 celdas. Y 0.6 m

máx X

mín 2.5 m

0.625 m

Figura 2: Esquema del dominio computacional. (No se muestra el refinamiento del cuadrado). Con ello, fui definiendo distintas longitudes del eje OX positivo del dominio en el archivo blockMeshDict, haciendo simulaciones para túneles de 0.9375, 1.25, 1.5625, 1.875, 2, 2.25, 2.5625, 2.8125 metros hacia atrás (a sumar 0.625 m de la parte anterior del túnel). Para decidir cuál sería la configuración más óptima, calculé el CD y saqué la media con cada longitud (véase código en el anexo A). Ahora fui comparando el % de variación de dichas medias con respecto a las longitudes de 1.5625, 1.875, 2, y 2.25 metros, y llegué a la conclusión que, partiendo un umbral de 2 % de diferencia, el de 1.875 tenía una media de los porcentajes de variación menor de entre las 4 longitudes, no haciendo más datos debido a que con los dos túneles más largos analizados la variación es muy pequeña, así que por economizar opté por el túnel de 2.5 m en total. A continuación se detallan los valores calculados y usados.

18

Longitud parte posterior (m)

CD

0.9375 1.25 1.5625 1.875 2 2.25 2.5625 2.8125

168,3626698 152,0802478 146,5984022 144,764155 147,6985533 142,6960067 140,8252618 143,7543066

Cuadro 4: Valores medios de CD para las diferentes longitudes. Longitud (m) 1.5625 1.875 2 2.25 2.5625 2.8125

% sobre 1.5625 m 0 1,251205432 -0,75045233 2,661963196 1,311000171 1,940059036

% sobre 1.875 m 1,267058943 0 -2,027019945 1,428632897 2,720903688 0,697581785

% sobre 2 m -0,744862491 -1,986748164 0 3,386997722 4,653594347 2,670470755

% sobre 2.25 m 2,73476154 1,44933862 -3,50573695 0 1,31100017 -0,74164652

Medias % 1,282755101 0,817431473 1,595890434 0,95069076 − −

Cuadro 5: Variaciones porcentuales11 de los valores de CD y media de las 4 longitudes referenciadas. Como se ve, en verde está indicado el porcentaje medio menor de variación del CD para las 4 primeras longitudes (se lee horizontalmente), interpretado como que para el túnel de 1.875 m de longitud posterior, las variaciones porcentuales medias del coeficiente para distintos tamaños del dominio serían menores con esa longitud, habiendo más saltos con las demás analizadas.

5.2. Malla Para el estudio de la sensibilidad de malla, la relación de aspecto de las celdas es crucial, por lo que se calculó la relación del dominio decidido anteriormente y se aplicó a cada celda para que fueran lo más cuadradas posible. Se decidió comenzar por un tamaño inicial de celda de 1 cm, aumentando dos veces y disminuyendo una el número de celdas en cada dirección en grupos de 20, con lo que quedarían tamaños iniciales de 230x40, 250x60, 270x80 y 290x100 celdas en sendos archivos blockMeshDict. Además, para cada uno de estos se realizaron distintas operaciones de refinado a la malla inicial con niveles 0 (sin refinar) y 1 para una caja de refinamiento (refinementBox), y 2 y 3 para la superficie del prisma (véase el apartado 3.1.5 para mejor entendimiento). Eventualmente se decidió también configuarar un nivel 2 de refinamiento en la caja, pero se descartó por la cantidad de recursos computacionales y tiempo necesarios al aumentar la velocidad, y también por ser un problema bidimensional que consideramos las mallas obtenidas como suficientemente fiables. Con ello se presentan los siguientes datos extraídos como medias de los valores del coeficiente de resistencia calculado en la superficie del prisma, y sus variaciones en % para cada malla. 11

Cálculo de la variación porcentual:

B % = 100 1 − , A

con A la variable sobre la que se calcula la variación, la de referencia, y B la variable que varió.

19

(35)

Tamaño malla

CD

Malla 230x40 rb 0 rs 2 Malla 230x40 rb 1 rs 2 Malla 230x40 rb 1 rs 3

136,234681150800 145,146465686001 172,421297139202

Malla 250x60 rb 0 rs 2 Malla 250x60 rb 1 rs 2 Malla 250x60 rb 1 rs 3

161,358762184321 177,908563021761 190,129225123040

Malla 270x80 rb 0 rs 2 Malla 270x80 rb 1 rs 2 Malla 270x80 rb 1 rs 3

174,788575214201 192,429907767201 185,169297404481

Malla 290x100 rb 0 rs 2 Malla 290x100 rb 1 rs 2 Malla 290x100 rb 1 rs 3

164,202959367399 182,068090377002 185,772561136896

Cuadro 6: Valores medios de CD para diferente número de celdas y niveles de refinamiento (rb ≡ refinementBox, rs ≡ refinementSurface).

Malla 230x40 Malla 250x60 Malla 270x80 Malla 290x100

Sobre la de 230x40

Sobre la de 250x60

Sobre la de 270x80

Sobre la de 290x100

0 18.4418 28.2996 20.5295

15.5703 0 8.32295 1.7627

22.0574 7.6835 0 6.0562

17.0327 1.7321 6.4467 0

Cuadro 7: Variaciones porcentuales de las medias de CD en mallas sin refinementBox.

Malla 230x40 Malla 250x60 Malla 270x80 Malla 290x100

Sobre la de 230x40

Sobre la de 250x60

Sobre la de 270x80

Sobre la de 290x100

0 22.57175 32.5764 25.4375

18.4151 0 8.1623 2.33801

24.5718 7.5463 0 5.3847

20.2790 2.2846 5.6912 0

Cuadro 8: Variaciones porcentuales de las medias de CD en mallas con refinementBox nivel 1 y refinementSurface nivel 2.

Malla 230x40 Malla 250x60 Malla 270x80 Malla 290x100

Sobre la de 230x40

Sobre la de 250x60

Sobre la de 270x80

Sobre la de 290x100

0 10.2701 7.3935 7.7434

9.3136 0 2.6087 2.2914

6.8845 2.6786 0 0.3258

7.1869 2.3457 0.3247 0

Cuadro 9: Variaciones porcentuales de las medias de CD en mallas con refinementBox nivel 1 y refinementSurface nivel 3.

20

Malla 230x40 rb 1 rs 2 rb 1 rs 3

0 18.7912

15.8187 0

Malla 250x60 0 6.8691

6.4276 0

Malla 270x80 0 3.7731

3.9211 0

Malla 290x100 0 2.0347

1.9941 0

Cuadro 10: Diferencias porcentuales de las medias de CD en las mismas mallas pero con distintos refinementSurface (niveles 2 y 3). Una vez presentados los datos en la tabla 6, se caracterizan para cada nivel de refinamiento12 con el fin de hallar valores medios de CD del orden del 2 % de diferencia o menos en las distintas mallas. Como se ve, en la tabla 7 se alcanza ese umbral, pero decidimos que es una malla bastante gruesa que no tiene caja de refinamiento que daría mayor fiabilidad a los resultados obtenidos en la estela (calle de Kármán). Donde se vuelven a tener valores en el rango aceptado es en la tabla 9, con valores francamente próximos para las mallas de 270x80 y 290x100 celdas con ese refinamiento, por lo que ese refinamiento será el adecuado para nuestro proyecto. Dado que hay dos opciones, se decidió realizar una última comparación para decidirnos exponiendo los porcentajes de variación del coeficiente de resistencia para, sobre el mismo número de celdas iniciales, los dos casos con caja de refinamiento (tabla 10), como modo de ver la estabilidad en los resultados con la variación del refinamiento. La malla que se dicidió como más fiable fue la de 290x100 con un nivel de refinamiento 1 en refinementBox y 3 en refinementSurface, con un total de 225890 celdas. 12

Está disponible una hoja de cálculo con valores más detallados aquí .

21

6. Resultados 6.1. Caso Re = 10 Para el primer caso, con una velocidad de 0.00118 m s−1 , una partícula fluida en movimiento rectilíneo uniforme tarda en recorrer la longitud total del túnel 2118.64 segundos ), por lo que ajustamos el tiempo final en el controlDict a 2200 con un ∆t = 0.1 s (t = l(túnel) v0 (22000 iteraciones), a sabiendas que debería dejar más tiempo pero con la finalidad de observar la evolución de los residuos iniciales hasta ese tiempo para ver si necesitaba continuar simulando. Como en fvSolution se ajustó el valor de tolerance en 1·10−8 (para todos los casos, además de relTol a 0, salvo para este que es 0.001 menos en la iteración final de la presión), cuando el residuo llega a esos órdenes de magnitud da por convergidos los cálculos, y se ve como la fig. 3 no llega a las 22000 iteraciones, parando sobre la 11 mil. Por tanto, se decide obviamente no proseguir la simulación más allá de 2200 s y se extraen los datos.

Figura 3: Gráfico con los valores de los residuos iniciales de esas tres variables. También se representarán los valores de los residuos iniciales (fig. 4), lo cuales responden de la resolución de las ecuaciones de las variables fluidas por el algoritmo PISO y la ecuación de Piosson. Estos vemos como, una vez superada la figura por el flujo, se estabilizan en valores por debajo del umbral de tolerance de 1·10−8 .

22

Figura 4: Gráfico con los valores de los residuos finales de esas tres variables.

Figura 5: Gráfico con los valores de los coeficientes aerodinámicos calculados en la superficie del prisma. 23

Dado que los valores de CL y CM están solapados debido al dominio en OY, se muestra a continuación la zona ampliada, del orden de 10−4 .

Figura 6: Gráfico ampliado para ver sin solapamiento la evolución de CL y CM .

Figura 7: Campo de velocidades visualizado en el instante t = 1100 s. La animación completa se puede ver aquí . Como se ve, no se produce calle de Kármán y el flujo se estabiliza, por lo que no habría que calcular el número de Strouhal según el enunciado, pero ya se aprecia que a partir de los 1000 segundos el flujo podría considerarse estacionario por ser f ∼ 0. 24

6.2. Caso Re = 100 Subiendo 10 veces el valor de la velocidad a 0.0118 m s−1 , el tiempo que demora la partícula en recorrer la longitud se reduce en ese orden de magnitud hasta un t = 211.86 s, por lo que, y al contrario que en el caso anterior, sabiendo que la calle de Kármán comienza a aparecer en el entorno de un Re de 90 a 100 se ajusta primeramente un tiempo de simulación de 500 segundos esperando que para esos instantes ya haya comenzado a notarse el efecto de la estela. Y en efecto esto comienza sobre los 250 s, creciendo la amplitud hasta cerca de los 400 s cuando se estable el movimiento armónico simple. En cuanto al paso temporal, este será de 0.02 s, por lo que inicialmente nos quedarían 25000 iteraciones temporales. En cuanto a los residuos, observamos en la fig. 8 como reflejan estos también la periodicidad que adoptan los valores de las variables fluidas debido a dicha estela en el flujo. Además, vemos como sus valores quedan demasiado elevados respecto a lo configurado en fvSolution del orden de 10−8 , aun siendo residuos iniciales que reflejan operaciones de la matriz y no en la resolución de la ecuación de Poisson, que serían los residuos finales del archivo log.icoFoam y que tienen valores significativamente menores como se ve en la fig. 9 no superando el umbral de las tolerancias impuestas.

Figura 8: Gráfico con los valores de los residuos iniciales en cada iteración temporal (paso de tiempo). Cabe aclarar que el dominio del eje X sobrepasa la iteración 25 mil como habíamos indicado antes, y la razón, si bien viene explicada posteriormente, es que sobre el primer archivo log.icoFoam se trasladaron los datos de un segundo que llegaba hasta las 27500 iteraciones. 25

Figura 9: Gráfico con los valores de los residuos finales de la resolución de la ecuación (32) en cada paso de tiempo. Respecto de los coeficientes aerodinámicos, dado que existe calle de Kármán y hay que calcular la frecuencia del movimiento armónico del fluido, se decidió que el tiempo final de la primera simulación de 500 s era insuficiente ya que arrojaba pocos periodos, se determinó, partiendo del último valor temporal como condición inicial de esta segunda simulación, llegar hasta los 550 s. Como se definió en la ecuación (2), el número de Strouhal depende de la velocidad media de oscilación del fluido –caracterizada como el producto de la frecuencia y la longitud característica– y de su velocidad media –la incidente–. Tanto lc como v0 nos son conocidas, por lo que procederemos a hallar f . Para ello necesitamos unos valores que describan un movimiento armónico, y dado que el cálculo de los coeficientes aerodinámicos es un objetivo de este proyecto serán los datos que usemos. A la hora de elegir con cual de los 3 coeficientes calcular la frecuencia, llegamos a la conclusión que, sabiendo que para un cilindro la frecuencia de CD es el doble que la de los torbellinos (la de oscilación del fluido) y la frecuencia de CL es la misma, se podría usar el mismo razonamiento aunque no sea un cilindro nuestra figura, pero suponemos el concepto el mismo dado que la fuerza de resistencia no es sensible al signo de cada uno de los torbellinos de signo opuesto que se forman en cada periodo. Una vez teniendo los datos, habría que procesarlos con algún programa para obtener la frecuencia por medio del análisis en series de Fourier, un análisis de una señal discreta con un algoritmo basado en la transformada de Fourier como el FFT13 que está implementado 13

Fast Fourier Transform, ‘Transformada rápida de Fourier’. La función fft calcula la transformada de

26

en MATLAB 14 . El código utilizado está en el anexo B y devuelve 2 gráficas, siendo en la segunda (fig. 12) en la que se representa el valor de la frecuencia de oscilación de los valores de CL .

Figura 10: Representación gráfica de los valores de los coeficientes aerodinámicos. Fourier discreta de un conjunto de N datos introducidos en una variable (en este caso CL , ftt(Cl)), transformando el resultado a dominio de frecuencia con la entrada como dominio del tiempo. Esta transformación se define como

F=

N −1 ∑

− 2πi N kn

xn

∀k = 0, . . . , N − 1,

(36)

n=0

que para la implementación informática según el algoritmo de Cooley-Tukey hace uso de la posibilidad de escribirse como una matriz de Vandermonde de la siguiente manera   1 1 1 1 ··· 1 2 3 N −1 1  ω ω ω ··· ω  2 4 6 2(N −1)  2πi 1  1 ω ω ω · · · ω W =√  (37)  , ω = e− N .  .. .. .. .. . N . . . .  . . . . . 1 ω (N −1)

ω 2(N −1)

ω 3(N −1)

···

ω (N −1)(N −1)

De esta manera es como se implementa, en forma de ecuación matricial con

F = W x,

(38)

siendo x los valores a transformar introducidos. Más información en las siguientes referencias: [6], [7], [8], y [9]. 14 En la documentación del programa se puede encontrar su implementación, en la cual está basado el código aquí utilizado. Documentación web.

27

Figura 11: Representación de los valores de CL frente al tiempo. Una vez ejecutado el código y obtenida la frecuencia, procedemos a calcular número de Strouhal según la ecuación (2), St ≃ 0.222.

Figura 12: Representacion en dominio de la frecuencia de la trasformación, indicándose el valor máximo de la frecuencia en la señal de CL . 28

Figura 13: Campo de velocidades visualizado en el instante t = 350 s. La animación completa se puede ver aquí .

29

6.3. Caso Re = 300 Para esta velocidad de 0.0354 m s−1 , una partícula fluida tardaría 70.62 segundos en recorrer el túnel, y por conservadurismo se redujo solo en 100 s el tiempo final de la simulación, por lo que en esos 405 s (debido al ∆t se corrigió añadiendo 5 s más para tener un número natural de iteraciones, 27000) una partícula recorrería más de 5 veces la longitud completa del túnel de viento. Esto se aprecia como excesivo en la fig. 14, ya que los residuos se estabilizan desde la iteración 1300 aproximadamente, pudiendo haberse reducido el tiempo de simulación. El paso de tiempo para este caso se calculó de 0.015 s manteniendo Co < 0.5. Seguimos con el inconveniente de lo elevados que resultan los residuos iniciales, sobre todo los de la presión, caballo de batalla en todo nuestro trabajo. Aun así, los finales dan valores mucho más óptimos, y es que cabe la posibilidad de que, a medida que se aumenta la velocidad del fluido, le sea más difícil al solver lineal resolver la matriz –sobre todo para la variable de presión susceptible a gradientes altos– y por eso nos muestre mayores residuos que implicarían mayores errores. Una posible mejora podría ser aumentar el número de iteraciones para actualizar los valores de la presión mediante el valor de la línea nCorrectors (ver última parte del apartado 3.1.4), ya que en todos los casos se fijó en 2 (2 iteraciones internas en cada una temporal) salvo en el del flujo a Re 10 que está en 3. Es precisamente en ese caso donde vemos que esa iteración extra mejora en un orden de 10 el resultado de la presión, pero no se puede extrapolar a los demás casos ya que en ese no habría calle de Kármán y el fluido es más estable. Habría que probarlo para velocidades mayores, aunque por limitación temporal no se realizó en este proyecto menor.

Figura 14: Gráfico con los valores de los residuos iniciales de esas tres variables. 30

Figura 15: Gráfico con los residuos finales de las variables de la ecuación de Poisson en cada ∆t.

Figura 16: Gráfico con los valores de los coeficientes aerodinámicos. 31

En la fig. 16 se aprecia lo indicado en el caso anterior de cómo la frecuencia de los valores de CD es mayor que la de CL , probablemente relacionado con lo expuesto entonces.

Figura 17: Gráfica con la representación individual de CL . Obtenida la frecuencia como muestra la fig. 18 calculamos el número de Strouhal, que incrementaría algo su valor hasta St ≃ 0.225.

32

Figura 18: Potencia de la señal con el valor máximo de la frecuencia en la señal de CL . El eje X va hasta 33.¯3 s−1 , se muestra reducido para que se aprecien los picos en la señal.

Figura 19: Campo de velocidades visualizado en el instante t = 300 s. La animación completa se puede ver aquí . 33

6.4. Caso Re = 1000 Con una velocidad de 0.118 m s−1 , el tiempo empleado por una partícula fluida para alcalzar la salida en línea recta sería de 21.19 segundos, por lo que aquí se limitó drásticamente el tiempo de simulación a 40 s considerando con los anteriores casos que debería ser suficiente para alcanzar y obtener pasos temporales suficientes con el movimiento oscilante para calcular la frecuencia. El intervalo de tiempo para esta simulación se calculó ya sensiblemente menor que los anteriores, de 0.002 s, ya que se disparaba ya rápidamente el Courant si ajustábamos más por encima. De nuevo la fig. 20 nos muestra como los residuos iniciales siguen siendo elevados, sin embargo los de la velocidad se han reducido algo. En cuanto a los finales, también los de la velocidad han disminuido, apreciándose claramente un amortiguamiento en la dispersión de los mismos, lo que podría indicar que fue una simulación numéricamente más estable.

Figura 20: Residuos iniciales obtenidos en cada iteración temporal.

34

Figura 21: Residuos finales en cada paso de tiempo.

Figura 22: Representación de los tres coeficientes aerodinámicos. 35

Observamos como los resultados de los coeficientes se vuelven más erráticos alejándose de la forma de una onda senoidal a medida que se aumenta la velocidad, incluso teniendo varias frecuencias y desde luego dejando ya de mantener cierta amplitud constante.

Figura 23: Señal del coeficiente de sustentación obtenido. Lo descrito anteriormente se refleja en la salida de la fig. 24, en la cual observamos varios picos en los que la señal de CL sería más potente (varias frecuencias), escogiendo la mayor para determinar la frecuencia de oscilación del fluido. Vemos como ha aumentado bastante, lo que significa un movimiento oscilatorio más rápido que se debería reflejar en la animación de la fig. 25, viendo cómo se producen y disipan los torbellinos en la estela del fluido (calle de Kármán). Con estos datos procedemos a calcular nuevamente el número de Strouhal, St ≃ 0.216. Vemos como disminuye algo, y es que aunque f aumente, no lo hace en el mismo o mayor orden de magnitud que la velocidad lo hace (v0 ̸∝ St).

36

Figura 24: Valores de la transformada según el algoritmo FFT. El eje X va hasta 250 s−1 , se muestra reducido para que se aprecien los picos que aparecen en la señal.

Figura 25: Campo de velocidades visualizado en el instante t = 40 s. La animación completa se puede ver aquí , pero desafortunadamente no ajusté correctamente el intervalo de escritura y paraFoam solo es capaz de reproducir 4 saltos de tiempo quedando una animación muy pobre.

37

6.5. Casos Re = 200 mil En este régimen los efectos del flujo turbulento son evidentes, por lo que predominan las leyes de pared en la capa límite(apartado 1.1). Según estas, el valor máximo de la distancia a la pared de la capa sería de 300 (distancia adimensional, y + ), por lo que es un factor limitante a la hora de configurar las simulaciones. Para los dos modelos de turbulencia se realizó la simulación en las mismas condiciones del dominio computacional y malla de los anteriores casos, ignorando los valores de y + que resultaren. Como se ve en la tabla 11, estos sobrepasan el límite de la subcapa logarítmica (y por tanto la capa límite) en la paredes superior e inferior de nuestro dominio,consecuencia de que la celda menor en esas paredes tiene un tamaño excesivo. Contorno top bottom prisma_proyecto

K −ω 322.765 322.643 9.7339

K − ω ref. 161.949 161.798 9.73767

K − ω SST 308.833 305.988 9.86228

K − ω SST ref. 154.871 154.211 9.85747

Cuadro 11: Valores medios de y + para los diferentes modelos y mallas en el tiempo final. ref ≡ refinamiento paredes. Entonces se decidió refinar un nivel (dividir a la mitad el tamaño de la celda) la celda más cercana a las paredes top y bottom que son las conflictivas. Para ello se crearon 2 geometrías nuevas (refinementBox, véase apartado 3.1.5) con el grosor de una celda calculado teniendo en cuenta la altura del túnel (0.6 m) y el número de celdas en esa dirección (100). Entonces 0.6 el grosor de la zona refinada sería de 100 = 0.006 m, por lo que se fijarían las zonas de refinamiento a 0.294 m del origen. Tras esto ya se podrían considerar válidos los resultados dado que y + < 300. 6.5.1. Modelo K − ω A la velocidad de 23.57 m s−1 , el tiempo necesitado por una partícula fluida para recorrer el túnel sería de 0.105 segundos, configurándose un tiempo de simulación de cerca del doble, 0.2 s. En estos regímenes de velocidades del fluido, el número de Courant es crítico debido a dicha magnitud, por lo que el paso de tiempo tiene que ser muy pequeño, en nuestro caso de 1·10−5 s, y por lo tanto obteniendo un total de 20000 iteraciones temporales. Vemos en la fig. 26 como los residuos iniciales se van reduciendo paulatinamente aun no llegando a los valores de la tolerancia, por lo que la simulación, a pesar de que según ese criterio el programa no la considera convergida, nosotros sí aceptamos esos valores tan bajos como válidos y no continuamos la simulación más allá. Contrastan estos resultados con la tendencia que veníamos teniendo a lo largo de los últimos 3 casos, aquí los residuos no se estabilizan en torno a un valor al cabo de un tiempo tras un valle y ascenso, sino que continúan su descenso. Esto es debido al modelo de turbulencia usado para la resolución de las variables fluidas, K − ω no es de fiable en los efectos de la subcapa viscosa en las paredes del dominio computacional como el modelo K − ω SST (ver apartado 1.1), por lo que, y según muestra la animación 30, no existirá calle de Kármán al no poder simular el desprendimiento de la capa límite del fluido en la pared del prisma. El fluido, tras unos instantes iniciales de establecimiento de la coriente, se estabiliza en una corriente ∼ unidimensional cuasiestacionaria con un zona de mínima presión tras el prisma. Respecto de los residuos iniciales, estos sí siguen la tónica de mantenerse por debajo del valor de tolerance que arrojarían resultados fiables de las ecuaciones fluidas. 38

Figura 26: Residuos iniciales obtenidos en cada iteración temporal.

Figura 27: Residuos finales en cada paso de tiempo. 39

Figura 28: Representación de los tres coeficientes aerodinámicos.

Figura 29: Gráfico ampliado para ver sin solapamiento la evolución de CL y CM . 40

Como en el caso de Re 10, los coeficientes aerodinámicos se estabilizan tras un corto instante, aunque en la fig. 29 se ve como el de sustentación realiza un movimiento armónico con una amplitud muy pequeña, probablemente producida por cierta oscilación del fluido que en la animación se es incapaz de apreciar. De todos modos tendría valores despreciables como se ve.

Figura 30: Campo de velocidades visualizado en el instante t = 0.1 s. La animación completa se puede ver aquí . 6.5.2. Modelo K − ω SST Con este modelo se mantienen las condiciones iniciales del anterior, ya que no añade más variables turbulentas, solo implementa otras ecuaciones promediadas con sus correspondientes coeficientes de cierre, todo ello implícito. Lo que sí cambiamos aquí es el tiempo total de simulación, aumentándolo teniendo en mente la oscilación del fluido y por tanto la necesidad de obtener suficientes periodos para el cálculo de la frecuencia. Por ello, el tiempo final queda establecido en 0.5 segundos, con un número total de iteraciones temporales de 50 mil. Los residuos iniciales ahora sí tienen la forma que venían adoptando los otros casos, y es que la elección del modelo de turbulencia según las condiciones de trabajo se torna fundamental a la hora de obtener los resultados deseados. Aquí, la posibilidad que da el K − ω SST de simular los desprendimientos de la capa límite del fluido en el prisma brindan los efectos reales al mostrar, como consecuencia, la calle de Kármán que aparece. Como siempre, los residuos de la presión son los que adoptan mayores valores, cabiendo destacar los picos que se muestran en los de la tasa de disipación de energía turbulenta ω y que en la fig. 32 se ven como variaciones de los residuos finales, lo que da a entender que es una variable sensible a los cambios en el fluido.

41

Figura 31: Residuos iniciales obtenidos en cada iteración temporal.

Figura 32: Residuos finales en cada paso de tiempo. 42

Figura 33: Representación de los tres coeficientes aerodinámicos. En cuanto a los coeficientes, vemos, como ya se apunto anteriormente, como CD va adquiriendo frecuencias distintas mientras CL mantiene un periodo más o menos constante además de la amplitud, como se mustra más claramente en la fig. 34

Figura 34: Señal del coeficiente de sustentación obtenido. 43

Tenemos una señal de CL bastante estable en periodo, por lo que la fig. 35 no nos mostrá muchos picos de potencia, lo que sí tenemos es, como cabía esperar, una frecuencia que nuevamente ha aumentado en 56 unidades. Esto implica un movimiento oscilatorio más rápido que se debería reflejar en la animación de la fig. 25, viendo cómo se producen y disipan los torbellinos más rápidamente. Con estos datos calculamos el número de Strouhal, St ≃ 0.202. Vuelve disminuir, teniendo el caso mismo que con Re 1000 en el que el orden de magnitud de aumento de la frecuencia es menor que el de la velocidad del fluido incidente.

Figura 35: Valores de la transformada según el algoritmo FFT. El eje X va hasta 5·104 s−1 .

Figura 36: Campo de velocidades visualizado en el instante t = 0.24 s. La animación completa se puede ver aquí . 44

A continuación se muestran la animación del campo energía cinética turbulenta.

Figura 37: Energía cinética turbulenta visualizada en el instante t = 0.21 s. La animación completa se puede ver aquí .

45

ANEXOS A. Código gnuplot para extraer las medias del CD set print "Medias_Cd.dat" stats '230x40/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'D' nooutput stats '230x40_rb1_rs2/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'E' nooutput stats '230x40_rb1_rs3/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'F' nooutput stats '250x60/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'P' nooutput stats '250x60_rb1_rs2/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'Q' nooutput stats '250x60_rb1_rs3/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'R' nooutput stats '270x80/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'N' nooutput stats '270x80_rb1_rs2/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'O' nooutput stats '270x80_rb1_rs3/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'S' nooutput stats '290x100/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'A' nooutput stats '290x100_rb1_rs2/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'B' nooutput stats '290x100_rb1_rs3/postProcessing/forceCoeffs/0/forceCoeffs.dat' using 2 " %lf %*lf %lf %*lf %*lf %*lf" name 'C' nooutput

print "Medias de los valores del coeficiente de arrastre en las distintas mallas (rb: nivel refinamiento caja interior, rs: nivel refinamiento superficie cuadrado)" print "" print "Malla 230x40 rs 2 = ", D_mean print "Malla 230x40 rb 1 rs 2 = ", E_mean print "Malla 230x40 rb 1 rs 3 = ", F_mean print "" print "Malla 250x60 rs 2 = ", P_mean print "Malla 250x60 rb 1 rs 2 = ", Q_mean print "Malla 250x60 rb 1 rs 3 = ", R_mean print "" print "Malla 270x80 rs 2 = ", N_mean print "Malla 270x80 rb 1 rs 2 = ", O_mean print "Malla 270x80 rb 1 rs 3 = ", S_mean print "" print "Malla 290x100 rs 2 = ", A_mean print "Malla 290x100 rb 1 rs 2 = ", B_mean print "Malla 290x100 rb 1 rs 3 = ", C_mean

46

B. Código FFT MATLAB clear all load forceCoeffs100.txt; %Sustituir por el nombre de cada archivo t=forceCoeffs100(:,1); %t como variable de los tiempos (1.ª columna) Cl=forceCoeffs100(:,4); %Cl como variable de los valores de C_L (4.ª columna) figure(1) %Gráfica de C_L grid on plot(t,Cl) grid on axis auto title('Coeficiente de sustentación') xlabel('Tiempo (s)') ylabel('Amplitud de los valores del coeficiente') Fs = 1/(t(2)-t(1)); % Frecuencia de muestreo (inversa del periodo) l = length(t); Y = fft(Cl); P2 = abs(Y/l); P1 = P2(1:l/2+1); P1(2:end-1) = 2*P1(2:end-1); % Saca la potencia de la señal sin ruido f = Fs*(0:(l/2))/l; maxP1 = max(P1) maxP1f = f(P1==maxP1) figure(2) plot(f,P1) grid on title('Espectro de potencia de la señal Cl(t)') xlabel('f (1/s)') ylabel('|P1(f)|') axis([-1 25 -.05 .25]) hold on plot(f(P1==maxP1),maxP1,'+r') text(f(P1==maxP1),maxP1,[' (' num2str(f(P1==maxP1)) ' , ' num2str(maxP1) ')']) hold off

47

Referencias [1] E. B. Martín, “Apuntes teoría cfd y turbulencia.” 2018. [2] C. Nguyen, “Turbulence Modeling,” 2005. [3] D. Wilcox, Turbulence Modeling for CFD. No. v. 1 in Turbulence Modeling for CFD, DCW Industries, 2006. [4] F. R. Menter, “Two-equation eddy-viscosity turbulence models for engineering applications,” AIAA Journal, vol. 32, pp. 1598–1605, Aug. 1994. [5] K. Hoffmann and S. Chiang, “Computational fluid dynamics, vol. 3,” Wichita, KS: Engineering Education System, 2000. [6] Wikipedia, “Transformada rápida de Fourier.” [7] Wikipedia, “Transformada de Fourier discreta.” [8] Wikipedia, “DFT matrix.” [9] Wikipedia, “Cooley–Tukey FFT algorithm.”

48