Antonio Pascau Área de Mecánica de Fluidos-LITEC Universidad de Zaragoza

1. Método de volúmenes finitos en Aerodinámica Antonio Pascau Área de Mecánica de Fluidos-LITEC Universidad de Zaragoza

2. Indice (I) • Recordatorio de los diferentes términos de las ecuaciones • Bases del Método de Volúmenes finitos • Disposición de nodos y variables • Ecuación de difusión unidimensional • Condiciones de contorno • Errores de discretización • Ecuación de convección-difusión unidimensional • Diferencias centradas y alternativas • Difusión numérica • Diagrama de variable normalizada • Implementación esquemas alto orden • Relajación MVF en Aerodinámica

3. Indice (II) • Ecuación de convección-difusión con término fuente • Flujo no estacionario • Estabilidad • Método de Von Newmann • Mallas colocalizadas • Velocidad en la entrefase • Desacoplamiento presión-velocidad • Método PWIM • Método SIMPLE • Método PISO • Resolución del sistema de ecuaciones • Mallas no estructuradas MVF en Aerodinámica

4. Ecuaciones de la Mecánica de Fluidos (I) • Cambio temporal + Cambio convectivo = Causas que producen dicho cambio • Cambio temporal integrado en el volumen + Flujo convectivo en caras saliente - Flujo convectivo en caras entrante = Flujo difusivo en caras entrante - Flujo difusivo en caras saliente + Otras fuentes integradas en el volumen MVF en Aerodinámica

5. Ecuaciones de la Mecánica de Fluidos (II) Ecuación integral de cantidad de movimiento Ecuación integral de energía total MVF en Aerodinámica

6. Bases delmétodo de volúmenes finitos • Troceado del dominio de estudio en volúmenes 3D (superficies 2D o intervalos 1D). Cada volumen lleva asociado un nodo, habitualmente en su centroide. • Integración de las ecuaciones correspondientes para las variables de cálculo en cada uno de los volúmenes unitarios. • Estimación de los flujos en las caras de los volúmenes unitarios. Este proceso conlleva la incorporación de relaciones algebraicas entre las variables y sus derivadas en las caras, y los valores en los nodos que comparten la cara. • Estimación de las fuentes de cambio dentro del volumen en función de los valores nodales. • Ensamblaje de todas las relaciones algebraicas de los valores nodales en un sistema de ecuaciones. PROCESO DE DISCRETIZACIÓN MVF en Aerodinámica

7. Disposición de nodos y variables (I) Según la posición relativa de escalares y velocidades • Mallas decaladas (staggered grid): escalares en el centroide de las celdas y componentes de velocidad en las caras. • Mallas colocalizadas (colocated grid): escalares y componentes de velocidad en la misma localización. Según la posición de los escalares en la celda • Variables escalares localizadas en el centro de la celda (cell centred approach) o localizadas en los vértices (cell vertex approach) • Tanto en un caso como en el otro la conservación de la masa se establece en los volúmenes del escalar. MVF en Aerodinámica

8. Disposición de nodos y variables (II) Malla colocalizada Malla decalada MVF en Aerodinámica

9. Discretización de la ecuación de difusión unidimensional (I) Ecuación para una variable genérica f Primer paso: Integración en el volumen finito 0 = -Flujo difusivo saliente + Flujo difusivo entrante En ausencia de fuentes internas el flujo neto difusivo ha de ser cero MVF en Aerodinámica

10. Discretización de la ecuación de difusión unidimensional (II) Segundo paso: Establecimiento de relaciones algebraicas entre flujos y valores nodales Tercer paso: Ensamblaje de las ecuaciones en un sistema Cuarto y último paso: Resolución del sistema MVF en Aerodinámica

11. Imposición de condiciones de contorno Ecuación del nodo 1 Obsérvese la presencia de un nodo ficticio externo (ghost node), que no debe aparecer en el ensamblaje final. Para quitarlo de en medio asumimos una variación lineal de la solución en el entorno del contorno izquierdo. Del mismo modo en el contorno derecho MVF en Aerodinámica

12. Errores de la discretización Errores en los términos discretizados Ecuación modificada (Warming & Hyett 1984) Ecuación del error MVF en Aerodinámica

13. Ecuación de convección-difusión unidimensional (I) Cambio convectivo de f = Causas que producen el cambio (difusión molecular) Primer paso: Integración en el volumen finito MVF en Aerodinámica

14. Ecuación de convección-difusión unidimensional (II) Segundo paso: Establecimiento de relaciones algebraicas entre flujos y valores nodales Tercer paso: Ensamblaje de las ecuaciones nodales Cuarto paso: Resolución del sistema de ecuaciones MVF en Aerodinámica

15. Ecuación de convección-difusión unidimensional (III) Cuando se introduce el término convectivo alguno de los coeficientes puede ser negativo. Para estabilidad es necesario que La condición para que todos los coeficientes sean positivos Se dice que el número de Pèclet ha de ser menor que 2 para que el esquema de diferencias centradas sea estable MVF en Aerodinámica

16. Solución de la ecuación MVF en Aerodinámica

17. Alternativas a diferencias centradas (I) Diferencias desplazadas contraflujo de primer orden (upwind differencing, FOUD) Tercer paso: Ensamblaje de las ecuaciones nodales MVF en Aerodinámica

18. Difusión numérica Cuando se utilizan diferencias desplazadas contraflujo se comete un error de discretización que puede ser interpretado como una difusión numérica representada por un coeficiente dado por Y el número de Pèclet de la discretización se puede definir como Cuando Pèclet real es infinito Pèclet de la discretización es 2 MVF en Aerodinámica

19. Alternativas a diferencias centradas (II) Diferencias desplazadas contraflujo de segundo orden, SOUD Tercer paso: ensamblaje de las ecuaciones nodales Matriz de coeficientes pentadiagonal MVF en Aerodinámica

20. Alternativas a diferencias centradas (III) Quadratic Upstream Interpolation for Convection Kinematics (QUICK) Ajusta una parabóla (polinomio de orden 2) a la solución tomando dos puntos aguas arriba y un punto aguas abajo. Es el esquema más preciso dentro del grupo de los lineales MVF en Aerodinámica

21. Diagrama de variable normalizada (I) Todos los esquemas presentados son susceptibles de ser representados en función de una variable normalizada entre el nodo i+1 y el nodo i-1 Cada uno de los esquemas se puede escribir en términos de una función MVF en Aerodinámica

22. Diagrama de variable normalizada (II) MVF en Aerodinámica

23. Diagrama de variable normalizada (III) Existe un punto dentro del diagrama por el que pasan todos los esquemas de orden 2 (aproximan perfectamente hasta la primera derivada en el entorno i+1, i-1) Si queremos que el esquema aproxime perfectamente hasta la primera derivada el término que la contiene debe desaparecer al hacer el desarrollo en serie de Taylor Sustituyendo b en la primera expresión en función de a y g, y normalizando Todos estos esquemas pasan por el punto (0.5, 0.75) MVF en Aerodinámica

24. Diagrama de variable normalizada (IV) En el diagrama de variable normalizada se pueden proponer esquemas no lineales que pasen por (0.5,0.75) HLPA (Hybrid Linear Parabolic Approximation) NOTABLE (New Option for the Treatment of Advection in the Boundary Layer Equations) MVF en Aerodinámica

25. Implementación de esquemas de alto orden La implementación práctica de estos esquemas se realiza por medio del Factor de Ponderación Direccional FPD (Downwind Weighting Factor, DWF) Para los esquemas de alto orden no lineales este factor vale Dado que la matriz de coeficientes del esquema upwind es la más estable los esquemas de alto orden se introducen como upwind (a matriz de coeficientes) + corrección explícita (a término fuente) MVF en Aerodinámica

26. Relajación Muchas de las ecuaciones de la Mecánica de Fluidos, y en particular las de Navier-Stokes, son no lineales, lo cual obliga a controlar las variaciones de iteración en iteración. Esto se consigue con el denominado factor de relajación (o subrelajación) (underrelaxation factor UF) Expresión resultante MVF en Aerodinámica

27. Dependencia con la relajación y los valores iniciales MVF en Aerodinámica

28. Ecuación de convección difusión con término fuente El tratamiento del término fuente es sencillo Habitualmente el término fuente se valora explícitamente. En ocasiones si existe dependencia con respecto a la propia variable f se linealiza del modo siguiente Para que todo funcione bien SP ha de ser negativo MVF en Aerodinámica

29. Flujo no estacionario (I) La ecuación de convección-difusión no estacionaria es Ahora existe un término adicional y una nueva coordenada sobre la que integrar Se ha asumido implícitamente que el valor nodal es adecuadamente representativo del promedio en el intervalo MVF en Aerodinámica

30. Flujo no estacionario (II) El proceso se repite para los dos términos restantes pero ahora hay que especificar en qué momento temporal se evalúan las derivadas espaciales • Casos: • q=1 Euler explícito • q=0 Euler implícito • q=1/2 Crank-Nicholson El método explícito es condicionalmente estable para cualquier tipo de esquema convectivo, esto significa que el paso temporal ha de estar por debajo de un determinado valor para que el esquema funcione MVF en Aerodinámica

31. Flujo no estacionario (III) El método Euler implícito tiene el rango de estabilidad del esquema convectivo, exactamente igual que el método Crank-Nicholson Consideremos el caso de una ecuación con parámetros constantes, y llamemos a=G/r. Utilizando diferencias desplazadas para el término convectivo y centradas para el difusivo obtenemos El problema puede aparecer si el primer coeficiente en el lado derecho es negativo MVF en Aerodinámica

32. Flujo no estacionario (IV) La condición de estabilidad se puede concretar en O lo que es lo mismo MVF en Aerodinámica

33. Método de Von Newmann de análisis de estabilidad (I) El error se representa como una serie de armónicos de Fourier con una amplitud (compleja) dependiente del tiempo. La condición de estabilidad es equivalente a satisfacer una condición para que la amplitud no se amplifique en el tiempo El armónico m del error en un instante posterior es Insertando en la ecuación algebraica entre valores nodales MVF en Aerodinámica

34. Método de Von Newmann de análisis de estabilidad (II) Llamando El factor de amplificación resulta ser Y su módulo La condición de que el módulo sea menor que la unidad coincide con la de positividad de los coeficientes MVF en Aerodinámica

35. Comentarios generales sobre estabilidad lineal • La estabilidad lineal va indisolublemente ligada al método de discretización espacial y al temporal. No se pueden valorar por separado. • Tanto la condición de positividad de los coeficientes como de amplificación del error según Von Newmann, que dan los mismos resultados, son condición suficiente pero no necesaria para estabilidad lineal. • Ninguno de los métodos anteriores tiene en cuenta las condiciones de contorno. • Las condiciones de contorno aparecen en el ensamblaje final en la matriz de coeficientes y en otra derivada de la anterior: la de iteración. • Si se desea saber si un determinado problema va a converger hay que conocer un parámetro de la matriz de iteración (a definir), el radio espectral. Éste ha de ser menor que la unidad • Pero de nuevo, el radio espectral menor que la unidad es condición suficiente pero no necesaria para que el problema converja. MVF en Aerodinámica

36. Velocidad en la entrefase • Al estar las velocidades en la misma localización que los escalares el balance convectivo de entrada-salida en cada volumen se ha de realizar con velocidades que no están directamente disponibles. • El modo de determinar las velocidades en las caras en mallas colocalizadas es muy peculiar y viene influenciado por la imposibilidad de utilizar un promedio geométrico. • Si se utiliza como velocidad en la entrefase un promedio de los valores nodales se produce un desacoplamiento entre el campo de velocidad y el de presión. El gradiente de presión en la velocidad nodal tiene una diferencia2d-nodal y la divergencia de la velocidad tiene una diferencia 2d-nodal en la ecuación de la presión deducida de la de continuidad. MVF en Aerodinámica

37. Presión en la ecuación de cantidad de movimiento En mallas decaladas la localización de las velocidades y la presión es diferente. Al integrar el gradiente de presión en la ecuación de cantidad de movimiento aparece una diferencia 1d entre nodos de presión Por el contrario si se utilizan mallas colocalizadas al realizar la misma operación de integración en el volumen, aparece una diferencia 2d Si se realizara un promediado geométrico para el cálculo de la velocidad en la cara el desacoplamiento sería total MVF en Aerodinámica

38. Método PWIM (I) PWIM (Pressure Weighted Interpolation Method) es un método desarrollado para mantener el acoplamiento entre presión y velocidad Expresión general de las ecuaciones ensambladas para los nodos i e i+1 Ecuación ficticia del nodo de la entrefase Esta ecuación no existe, ya que no hay localizada ninguna velocidad en la entrefase. Sin embargo formalmente su estructura ha de ser idéntica a las ecuaciones reales para ui y ui+1 MVF en Aerodinámica

39. Método PWIM (II) Si se despejan las velocidades se obtiene La hipótesis fundamental en PWIM es MVF en Aerodinámica

40. Método PWIM (III) La expresión final para la velocidad en la entrefase es • Tiene dos partes: • ume es la media entre valores nodales • uce es la corrección debida al campo de presiones en el entorno del punto considerado La corrección ha de ser cero en determinadas circunstancias, por ejemplo, flujo completamente desarrollado MVF en Aerodinámica

41. Método PWIM (IV) La resolución de la componente de velocidad requiere relajación La ecuación final resultante es idéntica a la no relajada con un término extra que representa la contribución de la iteración anterior MVF en Aerodinámica

42. Método PWIM (V) La expresión final requiere el almacenamiento de la velocidad en las caras, del mismo modo que se almacena en los nodos, si no la solución final depende del factor de relajación La convección en caras (phi en OpenFoam) es una variable más del cálculo que debe ser almacenada MVF en Aerodinámica

43. Acoplamiento velocidad-presión. Método SIMPLE La presión no tiene una ecuación que nos sirva para su determinación. Es necesario obtenerla de la de continuidad, que en un dominio unidimensional sería la segunda ecuación que nos permitiría obtener el campo fluido, es decir, la velocidad y la presión. Llamemos u* a la velocidad obtenida después de resolver cantidad de movimiento y supongamos que existe un+1 que satisface continuidad MVF en Aerodinámica

44. Acoplamiento velocidad-presión. Método SIMPLE (II) El método SIMPLE empieza restando ambas ecuaciones, obteniendo una ecuación para las correcciones de velocidad • SIMPLE asume una serie de simplificaciones: • Se desprecia H’ • El producto de dos correcciones se desprecia frente a términos que sólo tienen una • A’ se desprecia frente a los gradientes de P’ MVF en Aerodinámica

45. Acoplamiento velocidad-presión. Método SIMPLE (III) Como en mallas colocalizadas existen tanto velocidades en las caras como en los nodos es necesario repetir la obtención de las correcciones de velocidad en las caras Si se sustituyen las expresiones obtenidas anteriormente para las correcciones de los valores nodales se tiene La relación entre correcciones de velocidad y gradiente de correcciones de presión en las caras es formalmente idéntica a la existente en los nodos MVF en Aerodinámica

46. Acoplamiento velocidad-presión. Método SIMPLE (IV) Del mismo modo que se escribe la expresión correspondiente a la corrección de velocidad en la entrefase ‘este’ de un volumen, se pueden escribir las correspondientes a las otras entrefases de un dominio bidimensional, la oeste, sur y norte. Los subíndices refieren a posiciones relativas de la variable con respecto al punto donde se evalúa la velocidad MVF en Aerodinámica

47. Acoplamiento velocidad-presión. Método SIMPLE (V) Supongamos que, como ya se ha mencionado antes, el campo un+1 satisface continuidad, eso significa que La sustitución de un+1 = u* + u’ y la expresión de u’ en términos de gradientes de P’ conduce a MVF en Aerodinámica

48. Algoritmo del método SIMPLE • Suponer un campo de velocidades y de presiones asignado a cada uno de los nodos del dominio • Iterar en la ecuación de cantidad de movimiento con un número fijado de iteraciones internas (inner iterations) • Ensamblar la ecuación de corrección de presión, determinando el término fuente (imbalance de masa) • Corregir la presión con P’ y la velocidad con el gradiente de P’ • Con las nuevas velocidades y presiones volver a 2 a lo largo de un número fijado de iteraciones externas (outer iterations) o hasta que el error sea menor que una tolerancia dada. MVF en Aerodinámica

49. Acoplamiento velocidad-presión. Método PISO La única diferencia entre el método SIMPLE y PISO estriba en una segunda corrección introducida por PISO que mejora la satisfacción de conservación de masa por el campo de velocidades en las caras. Recordemos que SIMPLE despreciaba H’. PISO también lo desprecia en el primer paso, y en ese sentido es idéntico a SIMPLE hasta ese momento, pero lo evalúa explícitamente en el segundo paso, lo cual representa una corrección adicional. MVF en Aerodinámica

50. Algoritmo del método PISO • Suponer un campo de velocidades y de presiones asignado a cada uno de los nodos del dominio • Iterar en la ecuación de cantidad de movimiento con un número fijado de iteraciones internas (inner iterations) • Ensamblar la ecuación de corrección de presión, determinando el término fuente (imbalance de masa) • Corregir la presión con P’ y la velocidad con el gradiente de P’ • Con las nuevas velocidades evaluar explícitamente H’ • Ensamblar una nueva ecuación para la nueva corrección de presión P’’ en la que el término fuente no sólo contiene el imbalance de masa • Corregir la presión con P’’ y la velocidad con el gradiente de P’’ • Con las nuevas velocidades y presiones volver a 2 a lo largo de un número fijado de iteraciones externas (outer iterations) o hasta que el error sea menor que una tolerancia dada MVF en Aerodinámica