Exposición y revisión del método SIMPLE con aplicaciones a flujos de geometría sencilla

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster en Ingeniería Aeronáutica Exposición y revisión del método SIMPLE con aplicaciones a flujos de geometría sencilla Autor: Aitor Robles Corpa Tutor: Miguel Pérez-Saborid Sánchez-Pastor Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Máster Máster en Ingeniería Aeronáutica Exposición y revisión del método SIMPLE con aplicaciones a flujos de geometría sencilla Autor: Aitor Robles Corpa Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor Titular Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Máster: Exposición y revisión del método SIMPLE con aplicaciones a flujos de geometría sencilla Autor: Aitor Robles Corpa Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos E ste documento supone para mí el fin de otra etapa muy importante en mi vida, quizás una de las más complicadas y duras. Supone, además, la reafirmación de los valores con los que me educaron. La prueba de que, si se quiere algo, con esfuerzo, ganas y sacrificio, se pueden alcanzar todas las metas que uno se proponga en esta vida. Por supuesto, no habría llegado hasta aquí solo, y son muchas las personas a las que me gustaría agradecer su apoyo, sus palabras y los buenos consejos. Empezando por mi familia, por los que están y por los que se fueron, y en especial a mi madre. Sin su paciencia, cariño, comprensión y valores vitales, no sería quien soy. A mi tutor, Don Miguel Pérez-Saborid Sánchez-Pastor, por haber confiado en mí una vez más para desarrollar este estudio. Gracias por su orientación y dedicación. A mis amigos, que son pocos, pero son mi segunda familia. En especial a Loly y Desirée, que pese a la distancia que nos separa, siempre os he sentido cerca para cualquier cosa. Seguro que nos esperan muchas cosas buenas en el futuro que seguir disfrutando juntos. Y no me puedo olvidar de A.R., sin su consejo, quizás no habría comenzado esta aventura que estoy a punto de terminar. A todos, muchas gracias. Aitor Robles Corpa Sevilla, 2024 I Resumen E l campo de la Dinámica de Fluidos Computacional (CFD) ha experimentado un crecimiento significativo, siendo vital en la ingeniería y en particular en la industria aeronáutica. Este trabajo investiga el método SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), como procedimiento base de las simulaciones de dominios fluidos. Se presenta una revisión exhaustiva de la literatura de CFD, seguida de la formulación matemática de las ecuaciones de Navier-Stokes y la implementación del algoritmo SIMPLE. Se resuelven casos de geometría sencilla y flujo incompresible, evidenciando su eficacia y versatilidad. Se discuten los resultados en comparación con la literatura existente, seguido de conclusiones y sugerencias para futuras investigaciones. Los anexos detallan el método line-by-line, las ecuaciones implementadas y los códigos en Matlab. Este estudio no solo profundiza en el método SIMPLE, sino que también sirve como recurso educativo introductorio en CFD, proporcionando las herramientas básicas necesarias para comprender y aplicar algoritmos de simulación de fluidos. III Índice de Figuras 1.1 Ejemplo de mallas estructuradas [12] 9 1.2 Ejemplo de mallas estructuradas tipo O y C [12] 10 1.3 Ejemplo de malla estructurada multi-bloque [12] 10 1.4 Ejemplo de mallas no estructuradas [12] 11 1.5 Mallado estructurado co-ubicado 12 1.6 Ejemplo de problema 1D 13 1.7 Ilustración del Checkerboard pressure oscillation 13 1.8 Mallado estructurado escalonado 14 2.1 Teorema de Gauss 19 2.2 Caras y normales de la celda de referencia 20 2.3 Algoritmo SIMPLE 29 2.4 Ilustración del sistema de ecuaciones de la fila j32 2.5 Aproximación de la interpolación de ˆue34 2.6 Pared izquierda sin deslizamiento 42 2.7 Capa límite 44 2.8 Adaptación de velocidades ueyvnen condiciones cercanas a la pared 46 2.9 Pared superior con deslizamiento 48 2.10 Entrada del fluido a velocidad uniforme 50 2.11 Salida del fluido a presión uniforme 52 2.12 Celda cerca de la salida del fluido a presión constante 53 3.1 Canal de sección constante 58 3.2 Condiciones de contorno del canal de sección constante 59 3.3 Residuos y campo de presiones Re =100 60 3.4 Residuos y campo de presiones Re =200 60 3.5 Residuos y campo de presiones Re =300 60 3.6 Residuos y campo de presiones Re =400 61 3.7 Residuos y campo de presiones Re =450 61 3.8 Campo de velocidades y perfiles, Re =100 63 3.9 Campo de velocidades y perfiles, Re =200 63 3.10 Campo de velocidades y perfiles, Re =300 63 3.11 Campo de velocidades y perfiles, Re =400 64 3.12 Campo de velocidades y perfiles, Re =450 64 3.13 Comparación de perfiles de velocidad horizontal, x=0.0175m65 3.14 Comparación de perfiles de velocidad horizontal, x=0.025m65 XI XII Índice de Figuras 3.15 Comparación de perfiles de velocidad horizontal, x=0.0375m65 3.16 Problema de cavidad impulsada 66 3.17 Condiciones de contorno del problema de cavidad impulsada 67 3.18 Residuos y líneas de corriente Re =100 69 3.19 Residuos y líneas de corriente Re =200 69 3.20 Residuos y líneas de corriente Re =500 69 3.21 Residuos y líneas de corriente Re =1000 70 3.22 Residuos y líneas de corriente Re =1500 70 3.23 Residuos y líneas de corriente Re =2000 70 3.24 Campo de presiones y velocidades Re =100 72 3.25 Campo de presiones y velocidades Re =200 72 3.26 Campo de presiones y velocidades Re =500 72 3.27 Campo de presiones y velocidades Re =1000 73 3.28 Campo de presiones y velocidades Re =1500 73 3.29 Campo de presiones y velocidades Re =2000 73 3.30 Comparación de perfiles de velocidad horizontal, x=0.05 m 74 3.31 Flujo en canal con cambio de sección 75 3.32 Condiciones de contorno del canal con cambio de sección 76 3.33 Residuos y líneas de corriente Re =50 78 3.34 Residuos y líneas de corriente Re =100 78 3.35 Residuos y líneas de corriente Re =150 78 3.36 Residuos y líneas de corriente Re =200 79 3.37 Residuos y líneas de corriente Re =300 79 3.38 Campo de velocidades y perfiles, Re =50 81 3.39 Campo de velocidades y perfiles, Re =100 81 3.40 Campo de velocidades y perfiles, Re =150 81 3.41 Campo de velocidades y perfiles, Re =200 82 3.42 Campo de velocidades y perfiles, Re =300 82 4.1 Comparación de perfiles de velocidad a la salida con el perfil teórico de Poiseuille 85 4.2 Comparación de longitudes de entrada 86 4.3 Nuevos resultados para Re =300 87 4.4 Determinación de la longitud de entrada para Re =300 88 4.5 Comparación con perfil de Poiseuille para Re =300 89 4.6 Comparación de la longitud de entrada para Re =300 89 4.9 Comparación de resultados entre el método MacCormack y Hou et al [10] 92 4.10 Comparación de perfiles de velocidad para Re =100 93 4.11 Comparación de perfiles de velocidad para Re =400 93 4.12 Esquema de separación del flujo [15] 94 4.13 Validación de datos, Hossain et al. [15] 95 4.14 Resultados para Re =350 96 4.15 Resultados para Re =400 97 4.16 Resultados para Re =500 98 4.17 Resultados de Hossain et al. [15] para Re =100 99 4.18 Resultados de Hossain et al. [15] para Re =350 99 4.19 Resultados de Hossain et al. [15] para Re =500 99 A.1 Mallado para resolución mediante el método line-by-line con barrido vertical 107 B.1 Celda interior de referencia 112 B.2 Pared izquierda sin deslizamiento 115 Índice de Figuras XIII B.3 Pared derecha sin deslizamiento 117 B.4 Celda cerca de la pared derecha 118 B.5 Pared inferior sin deslizamiento 119 B.6 Pared superior sin deslizamiento 121 B.7 Celda cerca de la pared superior 122 B.8 Pared superior con deslizamiento 124 B.9 Celda cerca de la pared superior con deslizamiento 125 B.10 Entrada del fluido a velocidad constante 127 B.11 Entrada del fluido a presión constante 129 B.12 Salida del fluido a presión constante 131 B.13 Celda cerca de la salida del fluido a presión constante 132 Índice de Tablas 3.1 Canal longitudinal constante: iteraciones y residuos en función de Re 61 3.2 Cavidad impulsada: iteraciones y residuos en función de Re 68 3.3 Cambio de sección: iteraciones y residuos en función de Re 77 4.1 Error relativo de la velocidad máxima 85 4.2 Error relativo de la longitud de entrada 86 4.3 Error relativo de la velocidad máxima para Re =300 88 4.4 Error relativo de la longitud de entrada para Re =300 88 4.5 Errores relativos de la posición del centro de rotación 91 XV Introducción E l campo de la Dinámica de Fluidos Computacional (CFD, por sus siglas en inglés) ha experimentado un crecimiento exponencial en su relevancia y aplicación en la industria aeronáutica y en diversos campos de la ingeniería. Este crecimiento se ha visto impulsado por la capacidad de los métodos numéricos para simular y predecir el comportamiento de dominios fluidos con una precisión cada vez mayor, lo que ha permitido optimizar diseños, reducir costos y acelerar el desarrollo de productos. El presente Trabajo Fin de Máster surge de la motivación de explorar y comprender a fondo cómo funcionan códigos comerciales de simulación de dominios fluidos como pueden ser ANSYS Fluent, Autodesk CFD o CFD++. En este sentido uno de los primeros métodos desarrollados fue el método SIMPLE (Semi-Implicit Method for Pressure-Linked Equations). Este método sirve como base para muchos de los programas comerciales disponibles en la actualidad. Por ende, este estudio se centra en analizar en detalle este algoritmo y comprender su implementación en la práctica de la simulación computacional de fluidos. El objetivo principal de este estudio es proporcionar una comprensión detallada de las ecuaciones y algoritmos subyacentes al método SIMPLE, así como su aplicación práctica mediante la implementación de un código computacional. Este código se utiliza para resolver tres casos de estudio representativos: un canal longitudinal de sección constante, un problema de cavidad impulsada y un canal con cambio de sección, con el fin de ilustrar la eficacia y versatilidad del método en diferentes configuraciones. Es importante destacar que, si bien existen numerosos recursos y herramientas comerciales disponibles para la simulación de flujos fluidos, la comprensión profunda de los fundamentos teóricos y la capacidad para desarrollar y modificar códigos personalizados son habilidades esenciales para cualquier estudiante o investigador interesado en este campo. Este trabajo pretende proporcionar una síntesis clara y concisa de la literatura existente sobre los aspectos teóricos y prácticos del CFD, con el propósito de orientar y facilitar el aprendizaje y desarrollo de habilidades en este ámbito. Es de esperar que el estudiante o investigador interesado pueda utilizar este proyecto para desarrollar sus propios códigos sin necesidad de recurrir a cajas negras como son los códigos comerciales. El documento se ha estructurado en seis capítulos que tratan distintos aspectos del estudio realizado. En el primer capítulo, se lleva a cabo una revisión exhaustiva de la literatura sobre Dinámica de Fluidos Computacional (CFD), explorando la historia de su desarrollo, clasificando los métodos utilizados y proporcionando una explicación conceptual del algoritmo SIMPLE y sus variantes, así como una discusión sobre el mallado del dominio. El segundo capítulo se centra en la formulación matemática fundamental, partiendo de la ecuación general de advección-difusión y desarrollando XVII XVIII Introduccion las ecuaciones de Navier-Stokes, junto con su discretización. Este capítulo es crucial para comprender los fundamentos del método SIMPLE, detallando cada paso del algoritmo y adaptando las ecuaciones a diversas condiciones de contorno. En el tercer capítulo, se presentan los tres problemas seleccionados para ser resueltos con el código implementado, describiendo sus características y condiciones de contorno, y analizando los resultados obtenidos. El cuarto capítulo se dedica a la discusión de estos resultados, comparándolos con los hallazgos de la literatura actual. Finalmente, el quinto capítulo recoge las conclusiones derivadas del estudio y propone líneas futuras de investigación. Además, el documento incluye tres anexos: el Anexo A aborda los fundamentos del método line-by-line necesario para resolver los sistemas de ecuaciones, el Anexo B proporciona un resumen de las ecuaciones implementadas del método SIMPLE y un compendio de los ajustes requeridos según las condiciones de contorno, y el Anexo C contiene los códigos desarrollados en Matlab. En conclusión, este estudio no solo se centra en profundizar en el método SIMPLE y su aplicación práctica, sino que también aspira a servir como un recurso introductorio y educativo en el fascinante campo de la dinámica de fluidos computacional (CFD). Más allá de simplemente presentar los fundamentos teóricos y prácticos del método SIMPLE, este trabajo tiene como objetivo proporcionar una base sólida para aquellos que se inician en el campo de la CFD, brindándoles las herramientas necesarias para comprender, implementar y adaptar algoritmos y técnicas de simulación de fluidos. 1 Revisión bibliográfica A ntes de abordar cualquier proyecto en ingeniería, es esencial realizar una investigación preliminar sobre el tema en cuestión, familiarizándose con los conceptos y terminología empleada en dicho campo. El propósito de este primer capítulo es hacer una revisión de la bibliografía actual sobre la Dinámica de Fluidos Computacional (CFD), para proporcionar al lector una base sobre los fundamentos de esta disciplina, asegurando así una mejor comprensión del presente estudio. Se comenzará por un breve repaso histórico sobre el desarrollo de los métodos CFD, para después pasar a clasificarlos según el tipo de metodología empleada. Posteriormente, se centrará la explicación en el método SIMPLE, que es el que se desarrollará más en profundidad durante todo este estudio, y se abordarán sus variantes. Finalmente, se dará una visión sobre los tipos de mallados y se justificará el uso del mallado co-ubicado frente al mallado escalonado. La dinámica de fluidos computacional (CFD) es una herramienta de análisis utilizada por los ingenieros desde hace casi cincuenta años. Como explica Charesworth [ 5 ], en sus inicios, el CFD se aplicó a problemas de aerodinámica, pero pronto se comprendieron sus ventajas y se empezó a emplear en todas las disciplinas de la ingeniería para resolver una gran variedad de problemas. El gran logro del CFD es su capacidad para predecir, con una precisión razonable, los detalles de complejos patrones de flujo que afectan a problemas prácticos de ingeniería. Esto se consigue sin realizar pruebas con modelos en túneles de viento, tanques de agua o instalaciones experimentales similares. En su lugar, el CFD utiliza la capacidad de los ordenadores para resolver las ecuaciones que rigen la dinámica de fluidos mediante técnicas numéricas. Resulta evidente que los cálculos numéricos pueden realizarse en un ordenador más rápidamente y con menos esfuerzo que los tradicionales ensayos experimentales con modelos o maquetas. Sin embargo, la fiabilidad de estos cálculos está condicionada por su validación experimental dentro de ciertos parámetros y condiciones. Aun así, tal como recalca Charesworth [ 5 ], el CFD ofrece una herramienta de análisis mucho más rentable que las pruebas experimentales. Esta eficiencia permite una reducción significativa del tiempo necesario para analizar un problema de flujo de fluidos, lo que facilita un estudio más profundo y una mejor optimización del diseño técnico. Aunque los ensayos con modelos siguen y seguirán siendo una parte importante de la ingeniería, el CFD se ha convertido cada vez más en una herramienta más accesible y esencial que ya forma parte de los procesos normales de diseño y optimización. Como es lógico, el desarrollo de estas herramientas está vinculado al avance de la capacidad computacional de los ordenadores. En los próximos años, es de esperar que, con la incorporación de la inteligencia artificial y la computación cuántica, el CFD experimente un avance significativo tanto en precisión como en reducción de tiempos de cálculo. Este crecimiento proyectado en las capacidades del CFD promete transformar aún más la ingeniería, ofreciendo soluciones más rápidas, precisas y económicas para el análisis y diseño de sistemas complejos de flujo de fluidos. 1 8Capítulo 1. Revisión bibliográfica A continuación, se revisan y explican las diferencias fundamentales de cada una de las variantes. 1.4.1 SIMPLER El método SIMPLER (SIMPLE Revised) es una mejora respecto del método SIMPLE. Fue desarrollado también por Brian Patankar para aumentar la tasa de convergencia de su método predecesor. La principal diferencia es que se comienza por calcular unas pseudos-velocidades iniciales en base a un campo de presiones preliminar calculado de manera directa, en lugar de empezar por una suposición de ambos campos. Esto hace que la estimación inicial sea mucho mejor y, por tanto, el número de iteraciones hasta llegar a la solución, disminuye, consiguiendo una mejor convergencia del método. Los pasos simplificados del método SIMPLER: 1. Cálculo de pseudo-velocidades: 2. Resolución de la ecuación de presión 3. Resolución de ecuaciones de momento 4. Corrección de presión 5. Actualización del Campo de velocidades 6. Evaluación de la convergencia. 7. Si la solución no converge, se repite el proceso desde el paso 2. 1.4.2 SIMPLEC El método SIMPLEC (SIMPLE Consistent) es la segunda variante del método SIMPLE, que fue desarrollada por Van Doormal y Raithby [ 2 ]. Los pasos del algoritmo SIMPLEC son esencialmente los mismos que SIMPLE. La diferencia entre ambos métodos está en la forma que tienen de discretizar las ecuaciones de momento. Como se verá durante el desarrollo matemático, en un momento dado es necesario aproximar las correcciones de velocidades. Es en ese punto, donde ambos métodos difieren. En el caso de SIMPLEC, la expresión de aproximación se simplifica mucho más mediante la eliminación de los términos menos significativos. Con esta simplificación, se consigue un algoritmo más eficiente en ciertos casos, logrando una convergencia más rápida que SIMPLE. 1.4.3 PISO Por último, el método PISO (Pressure Implicit with Splitting of Operators) se considera una extensión del método SIMPLE. Se basa en un método que originalmente estaba desarrollado para el cálculo no iterativo de flujos compresibles no estacionarios. Issa lo adaptó y convirtió en un método iterativo para flujos estacionarios incompresibles. Se diferencia principalmente en que incluye 2 pasos de corrección para la presión, añadiendo pasos adicionales al algoritmo para mejorar la precisión y la convergencia. 1. Estimación Inicial del Campo de Presión y velocidades 2. Solución de las Ecuaciones de Momento 3. Corrección de la Presión 4. Actualización del Campo de velocidades 1.5 Tipos de mallado 9 5. Corrección de la Presión 6. Actualización del Campo de velocidades 7. Evaluación de la convergencia 8. Si no converge la solución, repetir desde el paso 2 con los nuevos valores de presión y velocidad. El algoritmo SIMPLE es muy sencillo conceptualmente y es la base de muchos códigos CDFs actuales. Sus variantes están enfocadas a mejorar la convergencia y ahorrar tiempo de computación. Sin embargo, Charlesworth [ 5 ] señala que las comparaciones han demostrado que el rendimiento de cada algoritmo depende de las condiciones de flujo, del grado de acoplamiento entre la ecuación de momento y las ecuaciones escalares, y de la cantidad de sub-relajación utilizada. Por lo tanto, dadas todas estas variables, es complicado afirmar que un método sea mejor que otro y no queda más remedio que elegir un algoritmo en particular según el caso. 1.5 Tipos de mallado Por último, pero no menos importante, es necesario echar un vistazo a los tipos de mallas usados en CFD. Como se ha explicado en esta revisión bibliográfica, el uso de una malla es fundamental para poder resolver el dominio fluido, pues es la forma en la que se debe subdividir el dominio en una serie de pequeñas celdas o volúmenes. Esta división es la que permite discretizar las ecuaciones que rigen el flujo de fluidos y convertirlas en un gran sistema de ecuaciones algebraicas simultáneas que pueden resolverse numéricamente [ 5 ]. En este sentido, como ya se adelantó en apartados anteriores, las mallas se suelen dividir en mallados estructurados y no estructurados. Figura 1.1 Ejemplo de mallas estructuradas [12]. Por su sencillez, las primeras técnicas CFD que se desarrollaron utilizaban mallas ortogonales estructuradas. Estas mallas se construyen típicamente a partir de cuadrados o rectángulos en dos dimensiones y cubos o prismas en tres dimensiones y, obviamente, son relativamente fáciles de construir. Su principal ventaja es que su estructura lineal facilita mucho su implementación en códigos informáticos para resolver las ecuaciones algebraicas resultantes del proceso de discretización. Sin embargo, aunque las mallas ortogonales estructuradas presentan ventajas significativas, éstas se ven mitigadas por el número limitado de dominios que pueden mallarse utilizando cuadrados o cubos. En la mayoría de las aplicaciones de ingeniería de interés práctico, los fluidos se encuentran en geometrías complejas que requieren técnicas de mallado más sofisticadas. 10 Capítulo 1. Revisión bibliográfica Los primeros intentos de encontrar técnicas de generación de mallas capaces de mallar geometrías más complejas seguían girando en torno a las mallas estructuradas. La clave fue que se eliminó la restricción ortogonal de la malla, lo que dio lugar a técnicas de generación de mallas libres de ajustarse a formas más complicadas. Los tipos de mallas que surgieron entonces se clasifican como mallas O, C y H. Algunos ejemplos se ilustran en la figura 1.2. Otra técnica que surgió basada también en mallas estructuradas, son las mallas multi-bloque, donde hay uno o más niveles de subdivisión. Normalmente en el nivel exterior, hay bloques generalmente grandes y se reserva el nivel más fino para las zonas de acoplamiento entre regiones o de especial interés en el estudio. Figura 1.2 Ejemplo de mallas estructuradas tipo O y C [12]. Figura 1.3 Ejemplo de malla estructurada multi-bloque [12]. Este tipo de mallas son más flexibles que las estructuradas y permite usar mayor resolución en aquellas regiones donde es necesario, aunque por contrapartida son más complejas de programar[ 12 ]. Para lograrlo, fue necesario replantearse por completo la forma de generar las mallas. El resultado fue el desarrollo de técnicas de generación de mallas no estructuradas que permiten mallar casi cualquier configuración geométrica actual. La mayoría de estas técnicas utilizan triángulos y tetraedros, en dos y tres dimensiones respectivamente, como elemento básico para construir las mallas. Así surgieron las mallas no estructuradas, como por ejemplo las mostradas en la figura 1.4. Evidentemente, este cambio de elemento básico obliga a desestructurar la malla, ya que mantener cualquier estructura lineal es imposible, a no ser que usen triángulos muy regulares, por supuesto, no es posible en absoluto con tetraedros. Sin embargo, este cambio es el que permite a las mallas no estructuradas discretizar geometrías muy complejas. Además, las mallas no estructuradas también 1.5 Tipos de mallado 11 facilitan enormemente el refinamiento adaptativo como parte del proceso de solución. Esta técnica consiste en utilizar información de la solución encontrada en una malla para crear una nueva malla con áreas localizadas más refinadas donde existen patrones de flujo complejos. Es evidente, que este método iterativo, permite que las ecuaciones de gobierno controlen lo fina o gruesa que debe ser una malla en una región determinada y, por lo tanto, conducen a soluciones más precisas. Tanto la flexibilidad geométrica, como la posibilidad del refinamiento adaptativo favorecen otra ventaja más sobre las mallas estructuradas, que es una reducción significativa del tiempo y el esfuerzo necesarios para generar las mallas. Según señala Charlesworth [ 5 ], a menudo, en los cálculos CFD basados en mallas estructuradas, el proceso de generación de la malla constituye en más de un cincuenta por ciento al tiempo total necesario para alcanzar una solución. Por lo tanto, cualquier técnica que reduzca el tiempo, supone una ventaja importante, puesto que eso permite bien realizar más simulaciones de un problema determinado o bien, reducir el coste de realización de un estudio específico. Figura 1.4 Ejemplo de mallas no estructuradas [12]. A primera vista, las mallas no estructuradas ofrecen muchas ventajas para una buena estrategia de generación de mallas, pero inevitablemente también presentan desventajas. El principal inconveniente es la falta de estructura lineal en las mallas, lo que suponen una mayor carga computacional por tener que almacenar los datos topológicos de la malla. Debido a esta falta de estructura, la discretización de las ecuaciones se vuelve mucho más compleja y requiere de nuevas técnicas de discretización para poder lograr los mismos sistemas algebraicos. Sin embargo, pese a ello, y gracias a los avances informáticos, las mallas no estructuradas son las que generalmente se emplean en la actualidad en cualquier programa CFD [9]. En conclusión, para cualquier estudio CFD, la elección de un mallado adecuado es fundamental a la hora de obtener una solución lo más exacta posible y, sobre todo, en términos de ahorro de tiempo, ya que no solo importa el tiempo empleado en la construcción del mallado sino también el tiempo de calculo que depende del mallado seleccionado. En este sentido y dados los objetivos del presente estudio, se ha optado por un mallado estructurado, ya que al tratarse de ser un manual introductorio a CFD, la sencillez de la discretización permite mostrar perfectamente el proceso de discretización de las ecuaciones NSE en volúmenes finitos. 12 Capítulo 1. Revisión bibliográfica 1.5.1 Mallado estructurado: escalonado vs co-ubicado En este punto, el lector debe tener claro que este estudio desarrollará el método de volúmenes finitos junto con el método SIMPLE, que permite el acoplamiento de las velocidades y presiones, para resolver las ecuaciones que gobiernan los dominios fluidos. Para este desarrollo, se empleará un mallado estructurado que permitirá discretizar las ecuaciones NSE en los volúmenes de control definidos, para posteriormente, crear sistemas de ecuaciones algebraicos que resolver computacionalmente. Sin embargo, antes de poder realizar esta discretización, y definir los nodos donde se establecerán las variables del dominio fluido, queda por abordar un problema que surge precisamente por desacoplar las velocidades y presiones. Es lo que se conoce como problema de oscilación de presiones en tablero de ajedrez, o en inglés, Checkerboard pressure oscillation. Los métodos de discretización en volúmenes finitos, empiezan por dividir mediante un mallado el dominio fluido en un conjunto de volúmenes discretos. Las variables de flujo como las velocidades, presiones, densidades, temperaturas, se vinculan bien en el centro de estos volúmenes o en sus vértices, y las ecuaciones del problema se expresan en términos de relaciones algebraicas que conectan estos puntos discretos entre sí. Es intuitivo asociar todas las variables de flujo en los mismos nodos en todo el dominio computacional. Este tipo de mallados como el de la figura 1.5 se conocen como mallados co-ubicados. Sin embargo, como apunta Charlesworth [5], aquí es donde surge el problema. Si las velocidades y las presiones se calculan en las mismas ubicaciones, se puede dar el caso de que un campo de presiones no uniforme pueda actuar como si de un campo uniforme se tratase en las ecuaciones de momento discretizadas. Figura 1.5 Mallado estructurado co-ubicado. El origen de este problema está en que, para resolver la ecuación de continuidad, es necesario interpolar las velocidades en las caras de los volúmenes y es esa interpolación la que provoca las oscilaciones de presión en forma de tablero de ajedrez. Para ilustrar este concepto, se empleará el mismo ejemplo que usa Mazumder [ 13 ]: un flujo unidimensional, de densidad constante con una velocidad de entrada igual a 1. 1.5 Tipos de mallado 13 Figura 1.6 Ejemplo de problema 1D. En este caso, la ecuación de continuidad para la celda O quedaría expresada como 1.5, donde se expresa que el flujo que entra por la cara wdebe ser igual al flujo que sale por la cara e. (ρue−ρuw)∆y=0(1.3) Pero las velocidades en las caras uw y ue , son desconocidas, puesto que las velocidades están discretizadas en los nodos centrales. Si se aproximan de manera lineal según las expresiones: uw=uO+uW 2ue=uO+uE 2 si se sustituye en la expresión de continuidad y simplificamos, quedaría como: uE=uW(1.4) Como se observa, el propio valor de la velocidad en O, no aparece en la ecuación, esto indica que matemáticamente, cualquier combinación en forma de tablero de ajedrez como la de la figura 1.7 sería una posible solución de la continuidad. Y dado que la continuidad es la que está asociada a la corrección de las presiones mediante el método SIMPLE, lo que se obtiene en realidad son oscilaciones de presiones para intentar converger a las infinitas soluciones posibles. Por lo tanto, el algoritmo no llegaría a converger y de hacerlo, no se puede garantizar que sea una solución con verdadero sentido físico. Figura 1.7 Ilustración del Checkerboard pressure oscillation . 14 Capítulo 1. Revisión bibliográfica Una de las primeras soluciones para remediar este problema de la interpolación debido al desacoplamiento, fue sugerida por Harlow y Welch, fue el uso de un mallado escalonado. El concepto era sencillo, se trataba de vincular las distintas variables de velocidad, presión y escalares en puntos distintos del dominio fluido, tal como muestra la figura 1.8. Figura 1.8 Mallado estructurado escalonado. Al escalonar estos nodos, se provoca que los volúmenes finitos empleados para discretizar las ecuaciones de cantidad de movimiento y continuidad, sean diferentes. Por lo tanto, si de nuevo se evalúa la continuidad en el volumen O, la expresión resultante es: (ρue−ρuw)∆y=0(1.5) uw=uOue=uE sustituyendo y simplificando: uO=uE(1.6) donde ahora sí se tiene en cuenta el valor de la velocidad de la propia celda O y, se elimina la necesidad de interpolación de las velocidades en las caras, consiguiendo una vinculación de continuidad real entre ambos volúmenes adyacentes. Es cierto que el mallado escalonado, consigue solucionar el problema del acoplamiento de presión y velocidad, pero también introduce una mayor complejidad en el procedimiento de discretización. Deja de ser intuitivo en cuanto a la interpolación de las variables y por lo tanto complica las expresiones y el algoritmo computacional resultante. Además, si se quisiera aplicar a mallas no estructuradas, sería aún más complejo. A raíz de esas complicaciones, surgió un nuevo enfoque para trata de mantener todas las variables en el mismo nodo, y por tanto seguir usando los mallados co-ubicados. Rhie y Chow [ 13 ] desarrollaron en 1983 el método PWIM (Pressure-Weight Interpolation Method). Se trata de un nuevo método de interpolación para obtener las velocidades en las caras de los volúmenes, pero que tiene en cuenta las presiones del resto de volúmenes adyacentes, por lo tanto, deja de ser una ponderación lineal que es la que provoca el problema de oscilaciones de presiones en tablero de 1.5 Tipos de mallado 15 ajedrez. Los detalles de este método serán desarrollados durante el capítulo 2 y se entenderá mejor por qué se tiene en cuenta la presión para la ponderación en la interpolación. En conclusión, PWIM es la última pieza que se necesita para comenzar a discretizar las ecuaciones NSE en una malla estructurada co-ubicada, y desarrollar el algoritmo SIMPLE para la resolución de dominios fluidos en CFD, alcanzando la convergencia hacia la solución. 2 Formulación matemática E n este capítulo, se presenta la base teórica y matemática sobre la cual se fundamentará el desarrollo computacional del proyecto. Se inicia con la ecuación general de advección-difusión, la cual constituye un modelo fundamental para la descripción de procesos de transporte en fluidos de cualquier magnitud que se quiera estudiar. Partiendo de esta ecuación, se realizará su discretización para un volumen finito, lo que permitirá observar las claves necesarias para su implementación en un marco computacional. Posteriormente, se realizará una particularización de esta formulación general para el caso específico de las ecuaciones de Navier-Stokes, que describen el comportamiento de fluidos viscosos e incompresibles. Estas ecuaciones, al ser de suma importancia en la mecánica de fluidos, servirán como base para el desarrollo del método de resolución computacional de dominios fluidos que se utilizará en el proyecto para los diferentes casos que se han propuesto. A continuación, se explicará en detalle el método SIMPLE (Semi-Implicit Method for Pressure Linked Equations), el cual es ampliamente utilizado para la resolución numérica de las ecuaciones de Navier-Stokes en problemas de flujo laminar e incompresible. Se expondrán los fundamentos teóricos de este método, así como cada uno de los pasos necesarios para obtener expresiones matemáticas que puedan ser implementadas en un código de programación. Finalmente, se abordará la revisión de las expresiones obtenidas para las diferentes condiciones de contorno que son comunes en problemas prácticos de flujo de fluidos. Se considerarán condiciones como paredes sin deslizamiento, velocidades y presiones de entrada y salida, con el objetivo de asegurar una adecuada modelización de los fenómenos físicos involucrados en el problema. Por lo tanto, este apartado es de fundamental interés pues es donde se sientan las bases teóricas y matemáticas necesarias para el desarrollo y la implementación computacional del modelo fluidodinámico propuesto en el proyecto. 2.1 Ecuación general de Advección-Difusión La ecuación de advección-difusión (2.1) describe la evolución en el tiempo y el espacio de una magnitud física φ en un medio continuo. En su forma más general, se compone de cuatro términos principales que representan diferentes fenómenos físicos presentes en el sistema: ∂ ∂t(ρφ)+∇·(ρuφ−Γ∇φ) = Sφ(2.1) 17 24 Capítulo 2. Formulación matemática Ec. Continuidad φ=1;Γ=0;Sφ=0(2.27) ∂ρ ∂t+∇(ρu) = 0(2.28) Ec. Cant. Movimiento φ=u;Γ=µ;Sφ= +∇p+ρ fm(2.29) ∂(ρu) ∂t+∇·(ρuu−µ∇u) = −∇p+ρ fm(2.30) A medida que se analiza su formulación, es crucial considerar las simplificaciones y suposiciones necesarias que permiten reducir la complejidad del problema y llegar a una descripción más manejable y aplicable a problemas específicos como los que aplican a esta revisión. En primer lugar, se considerará un flujo bidimensional isotérmico, lo que implica que las propiedades del fluido, como la densidad y la viscosidad, no varían en el tiempo ni en el espacio. Además, es lógico asumir que la gravedad es despreciable en comparación con otros efectos, lo que permite centrarse únicamente en las fuerzas inerciales y viscosas del fluido. Otra hipótesis fundamental típica es que el flujo se encuentra en régimen estacionario, lo que significa que las propiedades del fluido en cualquier punto del dominio no cambian con el tiempo, es decir, se puede despreciar el término transitorio de ambas ecuaciones. Si a esto se le agrega, que solo se van a estudiar flujos incompresibles, donde la densidad del fluido permanece constante en todo el dominio, implica que por lo tanto, las variaciones de velocidad se producen únicamente por cambios en la dirección y magnitud del flujo. Bajo estas simplificaciones, las ecuaciones de Navier-Stokes se reducen a tres ecuaciones principales: una de continuidad (2.31) y dos de cantidad de movimiento (2.32 y 2.33). ∇(ρu) = 0(2.31) ∇·(ρu u −µ∇u) = −∂p ∂x(2.32) ∇·(ρu v −µ∇v) = −∂p ∂y(2.33) La ecuación de continuidad representa la conservación de la masa y establece que la divergencia del campo de velocidad es cero. Por otro lado, las ecuaciones de cantidad de movimiento, derivadas de la segunda ley de Newton para fluidos, relacionan las variaciones de velocidad con las fuerzas que actúan sobre el fluido, como la presión y la viscosidad. Estas ecuaciones forman un sistema acoplado de tres incógnitas: las dos componentes de velocidad del fluido ( u=u i+v j ) y la presión p . Resolver este sistema permite determinar el campo de velocidades y la distribución de presiones en el dominio, lo que es crucial para comprender y predecir el comportamiento del flujo en cualquier dominio fluido de interés para aplicaciones en ingeniería. Al revisar este sistema de ecuaciones, se observa que en las ecuaciones de cantidad de movimiento, las tres incógnitas u , v y p aparecen de manera explícita. Sin embargo, en la ecuación de continuidad, 2.2 Ecuaciones de Navier-Stokes para flujos incompresibles 25 solo se hace referencia a la velocidad. Más adelante, en este estudio, se explicará cómo manipular la ecuación de continuidad para introducir la presión p , lo que permitirá resolver el sistema de ecuaciones utilizando el algoritmo SIMPLE. El próximo paso es discretizar estas ecuaciones en un volumen finito, utilizando las expresiones deducidas previamente de la ecuación general de advección-difusión. 2.2.1 Discretización y linealización de la ecuación de cantidad de movimiento para un volumen finito En primer lugar, se debe abordar la discretización de las ecuaciones de cantidad de movimiento en un volumen finito. Para ello, basta con tomar la expresión general discretizada previamente obtenida (2.26) y sustituir los mismos parámetros que se han utilizado anteriormente en (2.32 y 2.33) para las dos dimensiones. Esto permite expresar las ecuaciones de cantidad de movimiento en términos discretos (2.34 y 2.35 ). Donde como se puede observar, en el sistema de ecuaciones, siguen apareciendo las incognitas ˆu,ˆvyp. |(ρu)e|+(ρu)e 2ˆuO−|(ρu)e|−(ρu)e 2ˆuE−µe ˆuE−ˆuO ∆x∆y− −|(ρu)w|+(ρu)w 2ˆuW−|(ρu)w|−(ρu)w 2ˆuO−µw ˆuO−ˆuW ∆x∆y+ +|(ρv)n|+(ρv)n 2ˆuO−|(ρv)n|−(ρv)n 2ˆuN−µn ˆuN−ˆuO ∆y∆x− −|(ρv)s|+(ρv)s 2ˆuS−|(ρv)s|−(ρv)s 2ˆuO−µs ˆuO−ˆus ∆y∆x=−Zn sZe w ∂p ∂xdx dy (2.34) |(ρu)e|+(ρu)e 2ˆvO−|(ρu)e|−(ρu)e 2ˆvE−µe ˆvE−ˆvO ∆x∆y− −|(ρu)w|+(ρu)w 2ˆvW−|(ρu)w|−(ρu)w 2ˆvO−µw ˆvO−ˆvW ∆x∆y+ +|(ρv)n|+(ρv)n 2ˆvO−|(ρv)n|−(ρv)n 2ˆvN−µn ˆvN−ˆvO ∆y∆x− −|(ρv)s|+(ρv)s 2ˆvS−|(ρv)s|−(ρv)s 2ˆvO−µs ˆvO−ˆvs ∆y∆x=−Zn sZe w ∂p ∂ydx dy (2.35) Si se reorganizan los términos de la ecuación agrupando todos los términos que multiplican a cada una de las incógnitas, se consigue obtener los coeficientes de enlace, y una expresión mucho más compacta. Estos coeficientes AX i y AY i son fundamentales en la resolución de las ecuaciones discretizadas, ya que representan las contribuciones de los nodos vecinos al nodo central, es decir, determinan la forma en que las variables están relacionadas entre sí. Ahora sí, como se indicó anteriormente, es posible integrar el término fuente de estas expresiones. En el caso de las ecuaciones de Navier-Stokes y con las simplificaciones que se han realizado, este término fuente se reduce al gradiente de presiones, que es responsable de impulsar el flujo del fluido en cada dirección y por lo tanto contribuye a cada una de las velocidades. 26 Capítulo 2. Formulación matemática Tomando de nuevo el recinto de la celda de control e integrando: SX O=−Zn sZe w ∂p ∂xdx dy =−Zn s (pe−pw)dy = (pw−pe)∆y(2.36) SY O=−Zn sZe w ∂p ∂ydx dy =−Ze w (pn−ps)dx = (ps−pn)∆x(2.37) Dado que la presión en las caras no se conoce, sino que los valores se están calculando en el centro de cada una de las celdas, es necesario aproximar los valores de las caras mediante interpolación, utilizando una ponderación basada en la distancia de las celdas adyacentes. En principio, este método será válido siempre que las celdas estén posicionadas en el interior del dominio. Sin embargo, este será un punto importante a revisar cuando las celdas estén sujetas a condiciones de contorno, ya que esta aproximación podría dejar de ser válida. Por el momento, sin embargo, se considera válida. SX O= (pw−pe)∆y=pO+pW 2−pE+pO 2∆y=1 2(pW−pE)∆y(2.38) SY O= (ps−pn)∆x=pO+pS 2−pN+pO 2∆x=1 2(pS−pN)∆x(2.39) En última instancia, estas ecuaciones compactas y las expresiones de los coeficientes de enlace constituyen el núcleo implementado en el algoritmo SIMPLE para resolver las ecuaciones de cantidad de movimiento. A través de este proceso, se obtendrán los valores de velocidad horizontal ˆu y vertical ˆv. AX OˆuO+AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS=SX O(2.40) AX O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (2.41) AX E=−|(ρu)e|−(ρu)e 2+µe ∆x∆y(2.42) AX W=−|(ρu)w|+(ρu)w 2+µw ∆x∆y(2.43) AX N=−|(ρv)n|−(ρv)n 2+µn ∆y∆x(2.44) AX S=−|(ρv)s|+(ρv)s 2+µs ∆y∆x(2.45) 2.2 Ecuaciones de Navier-Stokes para flujos incompresibles 27 AY OˆvO+AY WˆvW+AY EˆvE+AY NˆvN+AY SˆvS=SY O(2.46) AY O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (2.47) AY E=−|(ρu)e|−(ρu)e 2+µe ∆x∆y(2.48) AY W=−|(ρu)w|+(ρu)w 2+µw ∆x∆y(2.49) AY N=−|(ρv)n|−(ρv)n 2+µn ∆y∆x(2.50) AY S=−|(ρv)s|+(ρv)s 2+µs ∆y∆x(2.51) 2.2.2 Discretización y linealización de la ecuación de continuidad para un volumen finito Queda por abordar la discretización de la ecuación de continuidad. Al seguir el proceso discutido hasta ahora, se debe tomar la ecuación diferencial (2.28) y llevar a cabo su integración a lo largo y ancho de la celda CFD. No obstante, también se puede abordar la discretización desde el punto de vista físico de esta ecuación. La ecuación de continuidad indica que debe existir un equilibrio físico entre los flujos que entran y salen del dominio fluido, en este caso de la celda. Por lo tanto, al discretizar la ecuación de continuidad, se debe evaluar el flujo que entra o sale por cada una de las caras este E, oeste W, norte Ny sur Sde la celda. ZVO ∇(ρu)dV =ZS (ρu)·n dA =0(2.52) ZS (ρu)·n dA =ρe(u·n)eAe+ρw(u·n)wAw+ρn(u·n)nAn+ρs(u·n)sAs=0(2.53) [(ρu)e−(ρu)w]∆y+[(ρv)n−(ρv)s]∆x=0(2.54) Al revisar la expresión, Se observa que la velocidad sigue apareciendo de forma explícita, sin embargo, no se encuentra rastro alguno de la presión en la ecuación de continuidad. Este asunto será abordado más adelante en este estudio, y se procederá a modificar la expresión para incorporar explícitamente la presión. Llegado este punto, ya se cuenta con las ecuaciones fundamentales de Navier Stokes expresadas en la forma adecuada para poder ser implementadas en cualquier código y poder ser resueltas de forma numérica. 28 Capítulo 2. Formulación matemática 2.3 Desarrollo del método SIMPLE La simulación de la dinámica fluidos computacional (CFD) es una herramienta crucial en ingeniería y ciencia para comprender y predecir el comportamiento de los fluidos en una amplia gama de aplicaciones, desde el diseño de aviones y automóviles hasta la optimización de procesos industriales. En este sentido, uno de los métodos más utilizados en CFD es el método SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), que proporciona una solución aproximada eficiente a las ecuaciones de Navier-Stokes para flujos incompresibles. 2.3.1 Fundamentos del método SIMPLE Resolver directamente las ecuaciones de Navier-Stokes es computacionalmente costoso debido a la complejidad de los términos no lineales y la necesidad de discretizar el dominio de simulación. El método SIMPLE aborda estos desafíos mediante un enfoque iterativo que simplifica la resolución del problema. Tal como indica Charlesworth en su tesis, la dificultad en resolver las ecuaciones de Navier-Stokes incompresibles radica en la falta de una ecuación explícita para el campo de presión. Este desafío se supera con el algoritmo SIMPLE, que emplea un procedimiento de aproximación y corrección para el campo de presión. Dicho de otra forma, se propone inicialmente un campo de presión y velocidades, para luego resolver las ecuaciones de momento y, posteriormente, se actualiza el campo de presión mediante un mecanismo adecuado. Este mecanismo consiste en aprovechar la ecuación de continuidad, ya que dada una presión correcta, una solución de las ecuaciones de momento también cumplirá con la ecuación de continuidad. En el contexto del algoritmo SIMPLE, esto se traduce en formular la ecuación de continuidad como una ecuación de corrección de presión, de ahí el nombre del método, que de forma semiimplicita, consigue una ecuación de continuidad expresada en presiones. Al resolver esta ecuación, se actualizan tanto el campo de presión como las velocidades utilizadas en la solución de las ecuaciones de momento. La derivación detallada de esta ecuación de corrección de presión se mostrará más adelante en el apartado 3.3.2.5. Durante el proceso de iteración, se garantiza la conservación de la masa y la presión en el dominio de simulación mediante la aplicación adecuada de los esquemas de corrección de velocidad y presión. Sin embargo, puede requerir un número significativo de iteraciones para alcanzar la convergencia y la solución puede ser sensible a la elección de parámetros y condiciones iniciales. 2.3.2 Algoritmo SIMPLE Como se puede observar en el diagrama de la figura 2.3, el método SIMPLE sigue un esquema de iteración entre la velocidad y la presión del fluido. Comienza con una suposición inicial para la velocidad y la presión del fluido en el dominio de simulación. Luego, se llevan a cabo tres pasos principales de iteración: predicción de velocidad, corrección de presión y corrección de velocidad. Estos pasos se repiten iterativamente hasta que se cumpla el criterio de convergencia seleccionado. A continuación, se detallarán minuciosamente cada uno de los pasos del algoritmo, explorando en profundidad el desarrollo de las ecuaciones requeridas para su implementación en un código. 2.3 Desarrollo del método SIMPLE 29 Primera iteración 𝑘 = 1 𝑢𝑘𝑣𝑘𝑝𝑘 Cálculo de coeficientes de enlace y fuentes para Ec. Cant. Movimiento 𝐴𝑂 𝑋𝐴𝑊 𝑋𝐴𝐸 𝑋𝐴𝑆 𝑋𝐴𝑁 𝑋𝑆𝑂 𝑋 𝐴𝑂 𝑌𝐴𝑊 𝑌𝐴𝐸 𝑌𝐴𝑆 𝑌𝐴𝑁 𝑌𝑆𝑂 𝑌 Cálculo de coeficientes de enlace y fuentes para Ec. Continuidad 𝐴𝑂 𝑝𝐴𝑊 𝑝𝐴𝐸 𝑝𝐴𝑆 𝑝𝐴𝑁 𝑝𝑆𝑂 𝑝 Resolución de las Ec. Cant. Movimiento en X e Y ො𝑢𝑂ො𝑢𝑊ො𝑢𝐸ො𝑢𝑆ො𝑢𝑁𝑅2 𝑋 ො𝑣𝑂ො𝑣𝑊ො𝑣𝐸ො𝑣𝑆ො𝑣𝑁𝑅2 𝑌 Interpolación de velocidades en las caras mediante PWIM ො𝑢𝑤ො𝑢𝑒ො𝑢𝑠ො𝑢𝑛 ො𝑣𝑤ො𝑣𝑒ො𝑣𝑠ො𝑣𝑛 Resolución de la Ec. Continuidad (Corrección de presiones) 𝑝𝑂 ′𝑝𝑊 ′𝑝𝐸 ′𝑝𝑆 ′𝑝𝑁 ′𝑅2 𝑝 Actualización de velocidades y presiones 𝑢𝑘+1 𝑣𝑘+1 𝑝𝑘+1 ¿Converge? 𝑅2 𝑋𝑅2 𝑌𝑅2 𝑝 Resultados y Postprocesado 𝑢 𝑣 𝑝 𝑘 = 𝑘 + 1 Figura 2.3 Algoritmo SIMPLE. 30 Capítulo 2. Formulación matemática 2.3.2.1 Primera iteración Como se han indicado anteriormente, el algoritmo SIMPLE comienza con la suposición de las velocidades uk y vk y presiones pk en el campo fluido, lo que representa una primera aproximación a la solución buscada para k=1 . Es importante que esta aproximación inicial tenga sentido físico y no esté demasiado alejada de la solución final, ya que una desventaja del método es su sensibilidad a las condiciones iniciales, lo que puede dificultar la convergencia hacia una solución. En función del problema a resolver, es crucial elegir las condiciones iniciales adecuadas. Sin embargo, esta variabilidad también puede ser útil para garantizar la consistencia de la solución. Si diferentes condiciones iniciales conducen a la misma solución del campo fluido, es probable que esta sea una solución sólida y confiable. Por ejemplo, en los casos tratados más adelante, se han considerado diferentes condiciones iniciales. Para un canal longitudinal y para el problema de cambio de sección, se propuso una velocidad uniforme como condición inicial, simulando un flujo uniforme en todo el canal, permitiendo que los efectos viscosos adapten gradualmente el campo fluido. Por otro lado, para el problema de cavitación, se optó por una velocidad inicial nula, permitiendo que la condición de contorno del movimiento de la cubierta superior y los efectos de advección y difusión moldeen el campo fluido con cada iteración. 2.3.2.2 Cálculo de los coeficientes y fuentes para las ecuaciones de cantidad de movimiento Siguiendo el diagrama, el siguiente paso implica el cálculo de los coeficientes de enlace AX i y AY i , y las fuentes SX O y SY O , necesarios para resolver las ecuaciones de cantidad de movimiento en los ejes X e Y . Es importante recordar que los coeficientes de enlace, derivados en secciones anteriores y expresados desde 2.41 hasta 2.45 para el eje X y 2.47 hasta 2.51 para el eje Y , representan las relaciones entre celdas adyacentes. Dado que estas expresiones se dedujeron a través de la discretización de la ecuación de cantidad de movimiento, las expresiones son válidas para celdas interiores. En principio, para celdas interiores, los coeficientes de enlace en los ejes X e Y tienen expresiones equivalentes. Respecto a las fuentes, solo tienen en cuenta los gradientes de presión entre las celdas, lo que impulsa al fluido a través de ellas y cuyas expresiones se encuentran en 2.38 y 2.39. En este punto, es crucial considerar las condiciones de contorno para una representación precisa del dominio fluido, pues no todas las celdas serán interiores. En el proceso de implementación del código, es fundamental tener en cuenta estas condiciones, dado que diferentes áreas del dominio estarán sujetas a distintas condiciones. Por consiguiente, es imperativo reflejar esta variabilidad en el código correspondiente para asegurar una representación precisa del sistema físico en estudio. En este caso, mediante diferentes bucles que recorran las distintas zonas según sus condiciones de contorno, por ejemplo recorriendo las celdas por filas o por columnas. Esta consideración implica la necesidad de tener en cuenta las condiciones de contorno que afectan tanto a las expresiones de los coeficientes como a las fuentes. Las adaptaciones específicas de las expresiones conforme a las condiciones de contorno serán abordadas en secciones posteriores. Un detalle a tener en cuenta en este apartado, específicamente durante la primera iteración, es la necesidad de interpolar la velocidad en las caras. Para este propósito, una aproximación inicial válida podría ser una interpolación ponderada por distancia, similar a la expresión (2.55). Es importante señalar que esto podría plantear interrogantes, ya que como se explicó previamente, esto puede dar origen al problema conocido como oscilaciones de presiones en tablero de ajedrez (checkboard 2.3 Desarrollo del método SIMPLE 31 pressure oscillation). Sin embargo, como se discutirá más adelante, este escenario está contemplado dentro del método de interpolación PWIM para la aproximación de las velocidades en las caras de las celdas. ue=uE+uO 2(2.55) 2.3.2.3 Resolución de las ecuaciones de cantidad de movimiento Llega el momento de la resolución por separado de las ecuaciones de cantidad de movimiento en los ejes X e Y . Aunque la metodología empleada es la misma, una nos proporcionará el campo de velocidades horizontales ˆu y la otra el campo de velocidades verticales ˆv . Sin embargo, es crucial tener en cuenta que el campo de velocidades obtenido de esta resolución no será un campo de velocidades convergente completo. Es decir, cumplirá únicamente con las ecuaciones de cantidad de movimiento, pero aún no satisfará la ecuación de continuidad. Para distinguirlo, estos campos se representan con los caracteres (ˆu)y(ˆv). A partir de la ecuación discretizada (2.34), podría pensarse que bastaría con resolverla para las incógnitas uO , uE , uW , uN y uS . Pero desde el punto de vista computacional, es más eficiente abordar la resolución por partes, aplicando un método ADI. En este estudio, se emplea el método Line-by-line, que resuelve primero las ecuaciones en dirección de las filas y luego en dirección de las columnas. Los detalles del método Line-by-line se pueden encontrar en el anexo A, explicado por Majumdar. Tal como se indica en la siguientes ecuaciones, en el primer barrido, las incógnitas son uO,uEyuW, mientras que en el segundo barrido, las incógnitas son uO,uN, y uS. AX Oˆ uO+AX Wˆ uW+AX Eˆ uE=SX O−AX NˆuN−AX SˆuS(2.56) AX Oˆ uO+AX Nˆ uN+AX Sˆ uS=SX O−AX WˆuW−AX EˆuE(2.57) Para hacer aún más eficiente la resolución del sistema, se realiza en términos correctivos. Esto implica calcular la corrección necesaria para cada velocidad para cumplir con la ecuación, lo cual está relacionado con el concepto de los residuos de la ecuación, que se utilizan como criterio de convergencia. De esta manera, si se introducen los términos correctivos u′ , se tiene a la izquierda las incógnitas a resolver y a la derecha el valor del residuo de dicha ecuación. AX Oˆuk+1 O+AX Wˆuk+1 W+AX Eˆuk+1 E=SX O−AX Nˆuk N−AX Sˆuk S(2.58) ˆuk+1 i=ˆuk i+u′ ii=O,E,W,N,S(2.59) AX O(ˆuk O+u′ O)+AX W(ˆuk W+u′ W)+AX E(ˆuk E+u′ E) = SX O−AX Nˆuk N−AX Sˆuk S(2.60) AX Ou′ O+AX Wu′ W+AX Eu′ E=SX O−AX Nˆuk N−AX Sˆuk S−AX Oˆuk O−AX Wˆuk W−AX Eˆuk E(2.61) RO=SX O−AX Nˆuk N−AX Sˆuk S−AX Oˆuk O−AX Wˆuk W−AX Eˆuk E(2.62) 32 Capítulo 2. Formulación matemática Debido tanto a la discretización como a la linealización de la expresión, la mayoría de los métodos añaden un factor de amortiguamiento inercial ωinercial que modera estas pequeñas correcciones en la velocidad, quedando la expresión donde el término de la derecha constituye los residuos de las ecuaciones. El valor de este amortiguamiento suele oscilar entre 0.2y0.8. (1+ωinercial)AX Ou′ O+AX Wu′ W+AX Eu′ E=RO(2.63) El proceso para obtener el campo de velocidades comienza calculando los residuos de todas las celdas antes de iniciar las iteraciones. Esto se hace para evaluar la convergencia del campo de velocidades mediante la raíz del error cuadrático medio (RMSE), una medida comúnmente utilizada para evaluar la calidad del ajuste del modelo. RX 2=v u u t n−1 ∑ i=2 m−1 ∑ j=2 R2 Oi,j(2.64) Una vez completado este paso, comienza una iteración interna, donde se resuelven los sistemas primero por filas y luego por columnas. La ecuación 2.63, es propia de cada celda, pero forma parte del sistema para una fila j , de manera que como resultado, cada fila k en realidad constituye un sistema tridiagonal, fácil de resolver para cualquier código que tenga implementada una función TDMA. Después de resolver el sistema de la fila k, se actualizan las velocidades con las correcciones obtenidas y se pasa a la siguiente fila. Una vez completadas todas las filas, se repite el proceso, pero esta vez resolviendo los sistemas en columnas. Esta iteración se conoce como iteración interna, ya que se utiliza para resolver una ecuación específica [13]. A x =d(2.65) A=         AO1,jAE1,j0 AW2,jAO2,jAE2,j AW3,jAO3,j ... ......AEn−1,j 0AWn,jAOn,j         ,x=         u′ O1,j u′ O2,j u′ O3,j . . . u′ On,j         ,d=         RO1,j RO2,j RO3,j . . . ROn,j         Figura 2.4 Ilustración del sistema de ecuaciones de la fila j. 2.3 Desarrollo del método SIMPLE 33 2.3.2.4 Método PWIM para la interpolación de las velocidades en las caras Hasta este punto, se han obtenido valores de velocidad en el centro de cada celda que cumplen con la ecuación de cantidad de movimiento, pero no con la ecuación de continuidad. Por lo tanto, la solución solo ha convergido parcialmente. Para lograr una convergencia completa, es necesario satisfacer también la ecuación de continuidad. Sin embargo, al analizar esta ecuación (2.54), se observa que los flujos están expresados en función de las velocidades en las caras de las celdas, no en las velocidades centrales de cada celda. Por lo tanto, se requiere algún tipo de interpolación para aproximar las velocidades en las caras. Aquí radica una de las claves del método. Como se explicó previamente, no cualquier tipo de interpolación es adecuada, ya que el uso de una interpolación basada en la ponderación de distancia puede dar lugar al problema del tablero de ajedrez, como se explicó en secciones anteriores. La clave para interpolar los valores de velocidad en las caras radica en el método desarrollado por Rhie & Chow en 1983 y posteriormente refinado por Miller & Schmidt en 1988 [ 13 ], el cual todavía es ampliamente utilizado en muchos códigos CFD en la actualidad. Este método, conocido como PWIM (Pressure-Weight Interpolation Method) en inglés, utiliza una ponderación por presión en lugar de por distancia. A continuación, se muestra su desarrollo detallado. La ecuación de cantidad de movimiento puede expresarse de manera más concisa si se agrupan todos los términos de las celdas adyacentes en un único sumatorio, como se muestra en 2.66. A continuación, se deriva una expresión para el valor de las velocidades uO y vO en el centro de la celda O. ∑ nb AX nb ˆunb =AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS(2.66) AX OˆuO+∑ nb AX nb ˆunb =SX O=pW−pE 2∆y(2.67) AY OˆvO+∑ nb AY nb ˆvnb =SY O=pS−pN 2∆x(2.68) ˆuO=−1 AX OO ∑ nb AX nb ˆunb +pW−pE 2AX OO ∆y(2.69) ˆvO=−1 AY OO ∑ nb AY nb ˆvnb +pS−pN 2AY OO ∆x(2.70) De manera similar, la expresión (2.72) para la celda adyacente Este seguiría una formulación parecida, adaptando los términos de enlace y las presiones del término fuente. Siguiendo este desarrollo, es posible suponer una expresión similar para la velocidad en la cara este ue como se muestra en (2.73). Hay que destacar que en esta expresión no es necesario incorporar la división por 2 en el término fuente, ya que el gradiente de presión en la cara de la celda es exacto y no aproximado como en las otras dos expresiones del centro de la celda. ˆuO=−1 AX OO ∑ nb AX nbOˆunb|O+pW−pE 2AX OO ∆y(2.71) ˆuE=−1 AX OE ∑ nb AX nbEˆunb|E+pO−pEE 2AX OE ∆y(2.72) ˆue=−1 AX Oe ∑ nb AX nbeˆunb|e+pO−pE AX Oe ∆y(2.73) 40 Capítulo 2. Formulación matemática uk+1 e=ˆue+ωuvu′ e=ˆue+ωuv "1 AX OE +1 AX OO#p′ O−p′ E 2∆y(2.118) vk+1 n=ˆvn+ωuvv′ n=ˆvn+ωuv "1 AY ON +1 AY OO#p′ O−p′ N 2∆x(2.119) uw|O=ue|W(2.120) vs|O=vn|S(2.121) 2.3.2.8 Evaluación de la convergencia de los residuos En la etapa final del proceso iterativo, es crucial evaluar la convergencia de los resultados, determinando si es necesario continuar iterando o si se puede considerar el proceso como finalizado. En el marco del algoritmo SIMPLE, este proceso se lleva a cabo mediante la evaluación de los residuos de las ecuaciones de cantidad de movimiento y continuidad, que se calculan al inicio de cada iteración en su correspondiente bloque. Es importante recordar que estos residuos representan la raíz del error cuadrático medio, una medida estándar de la discrepancia entre los resultados obtenidos y los esperados. Para determinar la convergencia, se comparan los valores de los residuos RX 2 , RY 2 y Rp 2 con una tolerancia predefinida, que en este estudio se establece en 10−10 . Si todos los residuos son inferiores a esta tolerancia, se detiene la iteración, indicando que se ha alcanzado una convergencia satisfactoria en el campo de velocidades y presiones. En caso contrario, se continua el proceso iterativo comenzando por la actualización de los coeficientes de enlace y fuentes de las ecuaciones de cantidad de movimiento. Es recomendable adoptar dos precauciones adicionales durante este proceso. En primer lugar, se sugiere visualizar los residuos a medida que avanza el algoritmo, lo que permite detectar si están convergiendo o divergiendo. En segundo lugar, se debe establecer un límite máximo de iteraciones del algoritmo para evitar posibles divergencias en los residuos. Este límite de iteraciones debe ajustarse según la complejidad del problema que se esté resolviendo, ya que puede variar considerablemente según la configuración específica del dominio fluido y las condiciones de contorno. En resumen, el algoritmo SIMPLE representa un proceso iterativo que gradualmente converge hacia una solución del campo de velocidades y presiones en un dominio fluido. Inicia con condiciones iniciales y se enfoca inicialmente en ajustar las velocidades, las cuales determinan los flujos que ingresan y salen de cada celda del dominio. A partir de los resultados de velocidad, se derivan los gastos no balanceados, que necesitan ajustarse mediante correcciones en las presiones. Por ejemplo, si en una celda el flujo de entrada es elevado respecto al flujo de salida, esto incrementará la presión, lo que a su vez influirá en las velocidades en la siguiente iteración. Este proceso de ajuste continúa iteración tras iteración, con la presión influyendo en las velocidades y viceversa, hasta que se alcanza un equilibrio. Este equilibrio se traduce en la convergencia del algoritmo, donde las velocidades y las presiones alcanzan valores estables iteración tras iteración. Alternativamente, si se establece un límite en el número de iteraciones y este límite se alcanza antes de lograr la convergencia, el proceso se detiene. 2.4 Particularización de las ecuaciones para las condiciones de contorno 41 2.3.3 Implementación del algoritmo SIMPLE en Matlab Para la implementación de este algoritmo, se ha diseñado un programa principal que actúa como punto central, donde se definen y gestionan los parámetros clave del problema a resolver. Estos parámetros incluyen la geometría del dominio fluido, las propiedades del fluido, el mallado co-ubicado estructurado, las tolerancias e iteraciones y los factores de amortiguamiento, entre otros aspectos fundamentales. La estructura principal del algoritmo SIMPLE se encuentra en este programa principal, el cual se encarga de llamar a cada bloque funcional representado en la figura 2.3 como una función independiente. Además de la lógica del algoritmo, el programa principal también incluye una funcionalidad para generar una representación gráfica dinámica de los residuos, a medida que se realizan las iteraciones. Esta visualización en tiempo real proporciona una manera efectiva de monitorear la convergencia del algoritmo y comprender el comportamiento del sistema durante su ejecución. Una vez que se completa la ejecución del algoritmo, el programa principal se encarga de almacenar los resultados obtenidos, incluyendo los campos de velocidades y presiones, así como los datos relacionados con los residuos y el número de Reynolds. Estos resultados son utilizados más adelante para análisis adicionales en procesos de postprocesado o para su visualización en herramientas externas. La integración de todas estas funcionalidades en un único programa ofrece una solución completa y eficiente para la resolución de problemas de dinámica de fluidos mediante el algoritmo SIMPLE. El código implementado en Matlab, que incluye tanto el programa principal como las funciones auxiliares, está disponible para su consulta en el Anexo C adjunto a este documento. Este recurso proporciona una referencia detallada y accesible para aquellos que deseen examinar el código fuente utilizado en la implementación del algoritmo descrito. 2.4 Particularización de las ecuaciones para las condiciones de contorno Una vez deducidas las expresiones generales que deben emplearse para las celdas interiores en el algoritmo SIMPLE, es momento de revisar estas ecuaciones y considerar cómo se deben adaptar según las diversas condiciones de contorno presentes en los dominios fluidos que se estén analizando. Estas adaptaciones deben implementarse en el código dentro del contexto de las celdas afectadas por cada condición de contorno específica. Para lograr esta particularización de las ecuaciones, es crucial comenzar con las expresiones generales no discretizadas y luego aplicar nuevas aproximaciones. No es viable derivarlas directamente de las ecuaciones discretizadas, ya que estas últimas incorporan otras aproximaciones que podrían no ser válidas para estos casos particulares. Se expondrán cuatro condiciones de contorno típicas: •Pared sin deslizamiento •Pared con deslizamiento •Entrada de fluido a velocidad constante •Salida de fluido a presión constante Cada una de estas condiciones requerirá ajustes específicos en las ecuaciones del algoritmo, los cuales se explicarán detalladamente para su implementación correcta en el código. 42 Capítulo 2. Formulación matemática 2.4.1 Pared sin deslizamiento Figura 2.6 Pared izquierda sin deslizamiento. La condición de pared sin deslizamiento, comúnmente denotada como "Wall" en muchos programas comerciales, es esencial en los dominios fluidos, ya que el fluido debe permanecer confinado dentro del recinto. Esta condición se caracteriza porque el fluido en contacto con la pared no presenta movimiento, es decir, tanto su velocidad horizontal como vertical son nulas. Respecto a la presión sobre la pared, aunque es desconocida, nunca será nula, lo que la convierte en una incógnita que debe aproximarse de alguna manera. En este caso particular, se considera una pared situada a la izquierda (figura 2.6), lo que implica que las velocidades en la cara oeste de la celda sean nulas. El proceso de desarrollo que se presentará a continuación puede aplicarse a cualquier otra pared, adaptando los términos según la ubicación específica de la misma. Para adaptar las ecuaciones, se seguirá el diagrama del algoritmo SIMPLE, revisando cada uno de los bloques y ajustando los términos según sea necesario. Ecuación de cantidad de movimiento Se comienza por examinar la adaptación de las ecuaciones de cantidad de movimiento (2.122 y 2.123 ) . Dado que la pared se encuentra a la izquierda w , se deben revisar los términos relacionados con el flujo Jwy la presión pw. JX e−JX w∆y+JX n−JX s∆x= (pw−pe)∆y(2.122) JY e−JY w∆y+JY n−JY s∆x= (ps−pn)∆x(2.123) Al considerar el flujo Jw , se recuerda que este término abarca tanto el flujo convectivo como el flujo difusivo, expresado tanto para el eje X (2.124) como para el eje Y (2.125). Al aplicar la condición de contorno de pared, las velocidades son nulas, lo que implica que no hay flujo convectivo y solo queda el término difusivo, el cual debe aproximarse de alguna manera. JX w=ρwuwuw−µw ∂u ∂xw (2.124) JY w=ρwuwvw−µw ∂v ∂xw (2.125) 2.4 Particularización de las ecuaciones para las condiciones de contorno 43 Una aproximación sencilla es utilizar una aproximación de Taylor de segundo orden basada en dos puntos. Por ejemplo, al desarrollar los términos para el eje X, donde uO y uE se expresan en función de uW , se puede combinar ambos para despejar el valor de la derivada parcial de u con respecto a xen w. uO=uw+∆x 2 ∂u ∂xw +(∆x/2)2 2! ∂2u ∂x2w (2.126) uE=uw+3∆x 2 ∂u ∂xw +(3∆x/2)2 2! ∂2u ∂x2w (2.127) Despejando ∂2u ∂x2w , ∂2u ∂x2w =8 ∆x2uO−uw−∆x 2 ∂u ∂xw(2.128) Sustituyendo en uE, uE=uw+3∆x 2 ∂u ∂xw +9∆x2 8 8 ∆x2uO−uw−∆x 2 ∂u ∂xw(2.129) Despejando ∂u ∂xw , ∂u ∂xw =9uO−uE−8uw 3∆x(2.130) Esta estimación del flujo difusivo (2.130) también será de gran utilidad en otras condiciones de contorno que se analicen más adelante en este estudio. Proporciona una forma general de abordar el flujo difusivo en las fronteras donde las velocidades son nulas, permitiendo adaptar eficazmente el tratamiento de los términos difusivos en diversas situaciones de contorno. JX w≈ − µw 9uO−uE−8uw 3∆x(2.131) Esto permite obtener una aproximación del valor del flujo JX w en función de uO , uE y uW , pero, dado que uW es nulo por la condición de la pared, JX w , queda como la expresión 2.132. El mismo razonamiento, sirve para JY w, quedando en función de vOyvE. JX w≈ − µw 9uO−uE 3∆x(2.132) JY w≈ − µw 9vO−vE 3∆x(2.133) Respecto al valor de la presión pw , como se ha indicando anteriormente, es desconocido pero se sabe que no es nulo, por lo que es necesario realizar una estimación. Mazumder [ 13 ] sugiere en este punto, emplear la teoría de la capa límite, la cual es propia de la mecánica de fluidos y se aplica en zonas cercanas a superficies donde los efectos viscosos son significativos, como es el caso que se 44 Capítulo 2. Formulación matemática está analizando. Por lo tanto, es apropiado considerar que dicha celda se puede encontrar dentro de la región conocida como capa límite. En esta teoría, se demuestra que la variación de la velocidad es mucho mayor en el eje perpendicular a la superficie que en el paralelo a ella. Es decir, según se ilustra en la figura (2.7), la componente vertical de la velocidad vexperimentará un mayor incremento que la componente horizontal u. Al relacionar este fenómeno con la ecuación de continuidad, se puede deducir que los gradientes de presión en dirección y son mucho mayores que los gradientes de presión en x . Esto significa que la presión en dirección xdesde la superficie varía muy poco dentro de la propia capa límite. Figura 2.7 Capa límite. ∂p ∂y≫∂p ∂x(2.134) En conclusión, una buena estimación podría ser considerar que la presión pw es igual a la presión en el centro de la celda pO . El error asociado a esta aproximación sería del orden de ∆x , cuyo valor suele ser pequeño debido al tamaño de la malla empleada. pw≈pO(2.135) Otro enfoque posible para estimar pw podría ser de naturaleza matemática, recurriendo nuevamente a una aproximación de Taylor de segundo orden basada en pO . Al desarrollar los términos y retener hasta el segundo término, es necesario aproximar el gradiente de presiones en O , el cual puede aproximarse mediante el gradiente en la cara e , que se conoce en función de las presiones pO y pE . Esto conduce a otra posible estimación, como se muestra en la ecuación 2.138 para el valor de pw . Al observar esta estimación, se ve claramente que la primera aproximación realizada, basada en la teoría de la capa límite, sería equivalente a considerar una aproximación de primer orden de Taylor. 2.4 Particularización de las ecuaciones para las condiciones de contorno 45 pw=pO−∆x 2 ∂p ∂xO (2.136) ∂p ∂xO ≈∂p ∂xe =pE−pO ∆x(2.137) pw=pO−∆x 2 pE−pO ∆x=3 2pO−1 2pE(2.138) Una vez realizadas las estimaciones para Jw y pw , simplemente se procede a sustituir estos valores en las ecuaciones de cantidad de movimiento y luego a reagrupar los términos para obtener los nuevos coeficientes de enlace AX i,AY iy las fuentes SX. AX O=|(ρu)e|+(ρu)e 2+µe ∆x+3µw ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (2.139) AX E=−|(ρu)e|−(ρu)e 2+µe ∆x+µw 3∆x∆y(2.140) AX W=0(2.141) Aprox. 1ºOrden SX O=pO−pE 2∆y(2.142) Aprox. 2ºOrden SX O= (pO−pE)∆y(2.143) AY O=|(ρu)e|+(ρu)e 2+µe ∆x+3µw ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (2.144) AY E=−|(ρu)e|−(ρu)e 2+µe ∆x+µw 3∆x∆y(2.145) AY W=0(2.146) Es evidente que los coeficientes de enlace afectados son únicamente los términos O , E y W . En el caso del coeficiente W , se anula debido a que no se incorpora flujo desde la pared. Sin embargo, en los coeficientes O y E , surge un nuevo término como resultado de la estimación del flujo difusivo en la pared. Además, la expresión de la fuente SX se altera, dependiendo de si se elige la aproximación de primer orden o de segundo. 46 Capítulo 2. Formulación matemática Interpolación de velocidades en las caras por el método PWIM El siguiente ajuste que se debe considerar, y que se ve influenciado por las condiciones de contorno, concierne a la interpolación PWIM para determinar el valor de las velocidades en las caras de la celda O . En este escenario específico donde el muro se ubica al Oeste, las velocidades en las caras Oeste yEste se ven afectadas. En el caso de ˆuw , por la condición de contorno, es nula. En cuanto a ˆue , es la más afectada, ya que, al revisar la expresión 2.78 deducida anteriormente con la interpolación de PWIM, esta depende del gradiente de presiones en la celda O , que se aproximaba como la diferencia ponderada de las presiones pW y pE . Sin embargo, ahora, la estimación realizada para pw es diferente. Por coherencia, se debe aplicar la misma aproximación anterior, ya sea de primer orden o de segundo orden. Por ende, la nueva expresión queda como: ˆuw|O=0(2.147) Aprox. 1ºOrden ˆue=ˆuO+ˆuE 2+1 2 pEE −pO 2AX OE ∆y+∆y 2AX OOpE−pO 2−"1 AX OE +1 AX OO#pE−pO 2∆y (2.148) Aprox. 2ºOrden ˆue=ˆuO+ˆuE 2+1 2 pEE −pO 2AX OE ∆y+∆y 2AX OO [pE−pO]−"1 AX OE +1 AX OO#pE−pO 2∆y (2.149) En este contexto, es crucial tener precaución en otros dos casos particulares que involucran celdas adyacentes a la pared, ya sea en el lado derecho o en la parte superior, dado que también se deben ajustar de manera especial las velocidades en las caras. Esto se debe a que, en estas celdas, como se muestra en las figura 2.8, el término afectado es el gradiente de presiones de la celda E o N , donde las presiones pee o pnn son desconocidas y se aproximan mediante una aproximación de Taylor de primer o segundo orden. Por lo tanto, siguiendo el mismo razonamiento anterior, se adaptan de la siguiente manera: Figura 2.8 Adaptación de velocidades ueyvnen condiciones cercanas a la pared. 2.4 Particularización de las ecuaciones para las condiciones de contorno 47 Aprox. 1ºOrden ˆue=ˆuO+ˆuE 2+∆y 2AX OEpE−pO 2+1 2 pE−pW 2AX OO ∆y−"1 AX OE +1 AX OO#pE−pO 2∆y (2.150) ˆvn=ˆvO+ˆvN 2+∆x 2AY ONpN−pO 2+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y (2.151) Aprox. 2ºOrden ˆue=ˆuO+ˆuE 2+∆y 2AX OE [pE−pO]+1 2 pE−pW 2AX OO ∆y−"1 AX OE +1 AX OO#pE−pO 2∆y(2.152) ˆvn=ˆvO+ˆvN 2+∆x 2AY ON [pN−pO]+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y(2.153) Corrección de presiones (Ec. de Continuidad) La corrección de presiones derivada de la ecuación de continuidad también experimenta ajustes debido a la presencia de la pared. Dado que el inicio del desarrollo de esta ecuación se centra en las velocidades de las caras, las cuales se expresan posteriormente como correcciones que involucran las presiones, es necesario analizar la velocidad en la cara oeste. Dado que esta velocidad es conocida y nula, su corrección u′ w debe ser igualmente nula. Por consiguiente, la expresión 2.154 se anula, eliminando así estos términos de los coeficientes de enlace Ap i . En cuanto a la fuente, el término ˆuw también debe ser anulado, ya que no hay flujo a través de la pared. Los ajustes en las expresiones quedan de la siguiente manera: u′ w="1 AX OW +1 AX OO#p′ W−p′ O 2∆y=0(2.154) Ap O=ρe∆y2 2"1 AX OE +1 AX OO#+ρn∆x2 2"1 AY ON +1 AY OO#+ρs∆x2 2"1 AY OS +1 AY OO# (2.155) Ap W=0(2.156) Sp O=−[(ρeˆue)∆y+ (ρnˆvn−ρsˆvs)∆x](2.157) 48 Capítulo 2. Formulación matemática Actualización de velocidades En el último paso, que implica actualizar las velocidades y presiones para la próxima iteración, el ajuste es sencillo, simplemente se debe imponer la condición de contorno de que la velocidad en la pared sea nula. Para el caso aquí presentado de la pared oeste, entonces: uk+1 wO=0(2.158) 2.4.2 Pared con deslizamiento Figura 2.9 Pared superior con deslizamiento. La siguiente condición de contorno que se desea analizar es un caso particular en el que la pared tiene deslizamiento. Esto significa que la pared se mueve de manera paralela al fluido, como se ilustra en la figura 2.9. Aunque no es una condición común, es importante estudiarla, ya que en los dominios posteriores planteados en este estudio, esta condición se presenta en la resolución de problemas típicos de cavitación. Lógicamente, si la pared se desliza paralelamente al fluido con una velocidad ulid , este movimien- to se transmitirá al fluido en contacto con la superficie debido a los efectos viscosos. Es importante tener en cuenta que la velocidad vertical vn sigue siendo cero, ya que de lo contrario implicaría un movimiento de la pared y, por consiguiente, un cambio en el volumen, lo cual no es el caso. Ecuación de cantidad de movimiento Al analizar los flujos JX n y JY n , se observa cómo efectivamente quedan los flujos difusivos debido a los efectos viscosos. Aprovechando la estimación del gradiente de velocidades realizada anteriormente, se observa que para el eje X aparece un nuevo término que incluye la velocidad de la pared, mientras que para el eje Y, los términos son los mismos que para una pared sin deslizamiento. JX n=ρnvnun−µn ∂u ∂yn =−µn ∂u ∂yn ≈ − µn 9uO−uS−8un 3∆y(2.159) JY n=ρwvnvn−µn ∂v ∂yn =−µn ∂v ∂yn ≈ − µn 9vO−vS−8vn 3∆y(2.160) 2.4 Particularización de las ecuaciones para las condiciones de contorno 49 JX n≈ − µn 9uO−uS−8ulid 3∆y(2.161) JY n≈ − µn 9vO−vS 3∆y(2.162) Por lo tanto, los cambios importantes para esta condición de contorno están relacionados principalmente con el término ulid , que se suma a la fuente de las presiones como una fuente más del sistema. AX OˆuO+AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS=SX O+µn 8∆x 3∆yulid (2.163) AX O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +3µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (2.164) AX N=0(2.165) AX S=−|(ρv)s|+(ρv)s 2+µs ∆y+µn 3∆y∆x(2.166) Interpolación de velocidades en las caras por el método PWIM En lo que respecta a la interpolación de PWIM, no se observan cambios significativos con respecto a las expresiones utilizadas para una pared sin deslizamiento. Simplemente es necesario aplicar la condición de contorno a ˆvn y tener la precaución, ya mencionada, para las celdas cercanas a la pared, donde es necesario aproximar la presión pnn utilizando las mismas expresiones que se utilizaron previamente (2.151 y 2.153). ˆvn=0(2.167) Corrección de presiones (Ec. de Continuidad) Para la corrección de las presiones, se mantiene la misma dinámica que en el caso de una pared sin deslizamiento. No se ven afectadas por el término ulid , por lo que su tratamiento es idéntico. Es decir, aplicando la misma lógica, el término v′ n es nulo, ya que no requiere corrección, lo que conlleva el ajuste de los términos AX O , AX N , AY O y AY N . Del mismo modo, debido a la ausencia de flujo desde la pared, se ajusta el término SP O. Actualización de velocidades En la actualización de las velocidades, se sigue un proceso similar. Simplemente se aplica la condición de contorno para la siguiente iteración, sin cambios significativos en comparación con el caso de una pared sin deslizamiento. 3 Resolución computacional de dominios fluidos U na vez cubiertos los fundamentos y el desarrollo de las ecuaciones que sustentan el algoritmo SIMPLE, se procederá a someterlo a prueba mediante una serie de casos en los que el dominio fluido presente una geometría simple. Siguiendo el enfoque propuesto en el desarrollo, el mallado debe ser estructurado, co-ubicado y ortogonal. Se ha optado por analizar tres casos clásicos: el flujo en un canal de sección constante, el problema de cavidad impulsada y el flujo en un canal con cambio de sección. Estos casos abarcan las diversas condiciones de contorno que se han desarrollado durante la particularización de las ecuaciones. Para todos los casos se estudiará el campo de presiones, el campo de velocidades, las líneas de corriente y los perfiles de velocidad en diferentes secciones según el número de Reynolds del dominio. Esta evaluación permitirá verificar la validez del método y del código implementado. 3.1 Flujo en canal de sección constante 3.1.1 Descripción del problema El problema de flujo en un canal de sección constante es un aspecto fundamental en la mecánica de fluidos, que se centra en el estudio del movimiento de líquidos o gases confinados dentro de un conducto o canal cuya geometría no varía a lo largo de su longitud. Su estudio es de gran importancia debido a su amplia gama de aplicaciones en diversas áreas de la ingeniería y la física. Comprende aspectos esenciales para el diseño y la operación eficiente de sistemas de transporte de fluidos, como tuberías en sistemas de suministro de agua, redes de distribución de gas o conductos de refrigeración. Comprender cómo fluyen los fluidos en estos canales es crucial para optimizar el diseño y el rendimiento de estos sistemas además de garantizar un transporte eficaz y seguro de los recursos. Para comenzar a aplicar el algoritmo SIMPLE y poder validar el modelo implementado, se ha optado por analizar el comportamiento de un fluido en un canal de sección constante como el que se muestra en la figura 3.1. Esto indica que tanto la forma como el área de la sección transversal, no cambian. El canal definido tiene una altura de H y una longitud de 5H , siendo la longitud característica H=0.01 m. Se decide implementar una configuración clásica del problema de mecánica de fluidos en la que se define una velocidad uniforme de entrada uin y una presión de salida pout que bien podría ser la presión atmosférica. Esto provocará el movimiento del fluido de izquierda a derecha, en cuyo avance, debido a los efectos viscosos es de esperar que los perfiles de velocidad vayan cambiando. 57 58 Capítulo 3. Resolución computacional de dominios fluidos Figura 3.1 Canal de sección constante. 3.1.2 Condiciones de contorno Revisando las condiciones de contorno del caso planteado, consta de dos paredes sin deslizamiento que son las que limitan la sección del canal. Al actuar como paredes, la velocidad del fluido que esté en contacto con la superficie deberá ser nula. El fluido entra por la cara izquierda del conducto con una velocidad uniforme uin , mientras que sale por la derecha con una presión de salida uniforme pout . En lo referente al mallado del dominio fluido, se ha implementado un mallado uniforme y coubicado. La longitud característica H se ha subdividido en 20 celdas cuadradas, cada una con una longitud de 0.5 mm. Esta división hace que el domino quede caracterizado con una matriz de 20x100 celdas. Como ya se indicó anteriormente en la particularización de las ecuaciones, hay que tener especial cuidado con las celdas que se encuentran cercanas a la salida y a la pared superior, puesto que se deben adaptar las aproximaciones realizadas para mantener la coherencia del planteamiento. Si se revisa el código en el anexo C, lo más sencillo es recorrer las filas y columnas en distintos bucles "For" diferenciando las zonas con condiciones de contorno y las celdas interiores. De forma detallada, se representa en la figura 3.2 las diferentes zonas descritas. Para realizar un análisis fluidodinámico más efectivo, se empleará el número de Reynolds como parámetro esencial. Esto permitirá una mejor comprensión y comparación de cómo el flujo se comporta en diferentes condiciones. Con este fin, se elegirán las propiedades del fluido de tal manera que el número de Reynolds pueda ser modificado fácilmente mediante la introducción del parámetro de velocidad de entrada uin . En este contexto, se definirán la densidad del fluido como ρ=1kg/m3 y su viscosidad como µ=10−5kg/(m·s) . Además, se fijará una longitud característica del problema, H=0.01 m. De esta manera, el único parámetro que variará será la velocidad de la entrada. Este último será ajustado en un rango que va desde 0.1 m/s hasta 0.4 m/s, lo que resultará en un intervalo de valores de Re comprendido entre 100 y 400. Re =ρuin H µ(3.1) 3.1 Flujo en canal de sección constante 59 Figura 3.2 Condiciones de contorno del canal de sección constante. •Zona 1: Entrada de velocidad uniforme •Zona 2: Pared superior sin deslizamiento •Zona 3: Salida de presión uniforme •Zona 4: Pared inferior sin deslizamiento •Zona 5: Celdas cercanas a pared superior •Zona 6: Celdas cercanas a la salida a presión uniforme. 3.1.3 Resultados y análisis En esta sección, se examinará la evolución del dominio fluido en relación con tres aspectos fundamentales: el campo de presiones, el campo de velocidades y los perfiles de velocidad, en función de la variación del número de Reynolds. Para poder validar el modelo implementado, es de esperar que a medida que aumente el Reynolds, los efectos viscosos sean menores y por tanto sus efectos sobre los perfiles de velocidad serán menos pronunciados. Para investigar este fenómeno, se realiza un análisis detallado de los perfiles de velocidad horizontal en distintas ubicaciones a lo largo del canal, partiendo desde su entrada. Este enfoque permitirá observar cómo la viscosidad del fluido influye en la distribución y la forma de los perfiles de velocidad a lo largo del recorrido del flujo, esperando obtener la clásica forma de campana cuando el flujo está completamente desarrollado. Otro aspecto crucial para validar el modelo, será la evolución de los residuos. Es fundamental verificar que el algoritmo utilizado para el análisis sea capaz de converger hacia una solución precisa y coherente. Para ello, se van a monitorear de cerca la evolución de los residuos a lo largo del proceso de simulación. La convergencia adecuada de estos residuos proporcionará una indicación confiable de que el algoritmo está progresando hacia una solución correcta y estable. Esta evaluación de la convergencia no solo asegurará la fiabilidad de los resultados obtenidos, sino que también permitirá comprender mejor la calidad y la precisión del modelo utilizado en este estudio. 60 Capítulo 3. Resolución computacional de dominios fluidos (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.3 Residuos y campo de presiones Re =100 . (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.4 Residuos y campo de presiones Re =200 . (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.5 Residuos y campo de presiones Re =300 . 3.1 Flujo en canal de sección constante 61 (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.6 Residuos y campo de presiones Re =400 . (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.7 Residuos y campo de presiones Re =450 . Tabla 3.1 Canal longitudinal constante: iteraciones y residuos en función de Re. Re Iteraciones Residuo u Residuo v Residuo p 100 122 6,14E-11 1,43E-11 8,39E-11 200 128 8,01E-11 1,84E-11 9,66E-11 300 152 1,98E-11 1,15E-11 9,59E-11 400 186 2,04E-11 1,68E-11 9,27E-11 450 204 2,10E-11 1,79E-11 8,78E-11 62 Capítulo 3. Resolución computacional de dominios fluidos Al analizar los residuos de todos los casos, se observa que descienden por debajo de la tolerancia establecida de 10−11 , lo que indica una convergencia satisfactoria en todos los resultados y en ningún caso se llega al máximo de iteraciones. Además, todos los residuos alcanzados, son del mismo orden. Inicialmente, los residuos experimentan un aumento hasta transcurrir unas cuantas iteraciones, momento en el cual los resultados comienzan a converger y, por ende, los residuos empiezan a disminuir. Como se refleja en los gráficos, la ecuación que presenta mayor dificultad para alcanzar la convergencia es la corrección de las presiones, que, como se recuerda, está relacionada con la ecuación de continuidad. Es probable que esto se deba a que es la ecuación que hace de nexo entre las dos ecuaciones de cantidad de movimiento. Durante las simulaciones, se evidenció que esta ecuación es un factor crítico para lograr la convergencia, ya que frecuentemente se detectó como la primera en no converger. Este contratiempo desencadenaba un efecto dominó en el resto de variables, generando resultados insatisfactorios. Se llegó a la conclusión de que es fundamental resolver con precisión la corrección de las presiones, pues su resolución inexacta afecta significativamente la estabilidad del proceso y a la validez de los resultados obtenidos en la simulación. A medida que el número de Reynolds aumenta, se observa cómo las curvas de los residuos muestran una mayor fluctuación, atribuible al aumento de la turbulencia del flujo y, por tanto, a mayores variaciones. Esto se refleja en los campos de presión, donde debido a esta turbulencia, los gradientes de presión entre diferentes zonas del dominio se incrementan, resultando en una mayor división del dominio. Numéricamente, se aprecia que, por ejemplo, en el caso de Re =100 , la diferencia de presión entre la entrada y la salida es de aproximadamente ∆p=0.01 Pa, mientras que para el último caso de Re =450 , asciende a ∆p=0.08 Pa. Es natural que esto afecte el número de iteraciones necesarias para alcanzar la convergencia, como se refleja en la tabla 3.1, pasando de 122 iteraciones hasta 204 para el último caso, lo que se traduce en un mayor tiempo de computación. A pesar de esto, en todos los casos, las líneas de corriente son prácticamente rectas, lo que indica un flujo laminar sin que se mezclen las capas de fluido, excepto en la entrada que son algo más curvadas hacia el centro. En las esquinas de la entrada, como era de esperar, se aprecia un fuerte gradiente de presiones debido a que el fluido debe frenarse por el contacto con las paredes debido a los esfuerzos viscosos. Para mantener el caudal, la presión debe aumentar para incrementar la velocidad del fluido de las capas adyacentes, lo que resulta en una aceleración del fluido en su parte central, adquiriendo la forma característica de campana, como se puede observar en las figuras 3.39, 3.41, 3.42, 3.11 y 3.12. Efectivamente, al revisar los campos de velocidades, se observa que los vectores de velocidad cercanos a la entrada poseen una componente vertical que dirige el fluido desde las paredes hacia el centro de la sección. A medida que el fluido avanza por el canal, la velocidad comienza a aumentar hasta que el flujo se desarrolla completamente. Para comprender mejor este fenómeno, se presentan los perfiles de velocidad horizontal tomados en diferentes secciones y se comparan con el perfil uniforme en la entrada y el perfil presentado en la salida. En todos los casos, se aprecia el desarrollo del flujo y el incremento progresivo de la velocidad desde las paredes hacia el centro. Inicialmente, la velocidad se incrementa cerca de las paredes debido al gradiente de presión en las esquinas de la entrada, que impulsa el fluido hacia el centro, mientras que la parte central ve aumentada su velocidad debido a los efectos difusivos de esas líneas de corriente, hasta que finalmente las líneas se unen en el centro e impulsan la línea central del flujo. 3.1 Flujo en canal de sección constante 63 (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.8 Campo de velocidades y perfiles, Re =100 . (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.9 Campo de velocidades y perfiles, Re =200 . (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.10 Campo de velocidades y perfiles, Re =300 . 64 Capítulo 3. Resolución computacional de dominios fluidos (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.11 Campo de velocidades y perfiles, Re =400 . (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.12 Campo de velocidades y perfiles, Re =450 . Si se comparan los perfiles de velocidad, se aprecia que a medida que la turbulencia en el flujo aumenta, las fuerzas inerciales predominan sobre las fuerzas viscosas. Este desequilibrio implica que el flujo requiere más tiempo para desarrollarse a lo largo del canal. Por lo tanto, resulta evidente que el flujo con un número de Reynolds de Re =100 se encuentra prácticamente desarrollado en la mitad del dominio, mientras que para Re =450 , el desarrollo del flujo se extiende a lo largo de casi toda la longitud del dominio hasta la salida. Este fenómeno se ilustra claramente en las figuras 3.13, 3.14 y 3.15, donde se comparan los perfiles de velocidad horizontal en diferentes secciones del canal ( x=0.0175 m, x=0.025 m y x=0.0375 m). Es notable cómo en los flujos que aún no están completamente desarrollados, la parte central del perfil de velocidad es casi uniforme, hasta que los efectos viscosos alteran la distribución, llevándola a una forma parabólica. Este fenómeno, también explica las fluctuaciones de los residuos al aumentar el Reynolds y el límite encontrado en Re =460 . Este límite en el Re encontrado durante las simulaciones se podría explicar porque con el dominio planteado, el flujo no termina de desarrollarse por completo, y por la tanto el problema no converge. Para poder solucionarlo, sería necesario ampliar la longitud del canal y mallar el dominio con celdas más pequeñas, aunque esto conllevaría a un aumento significativo del tiempo de simulación. 3.1 Flujo en canal de sección constante 65 En términos generales, las soluciones obtenidas en cuanto al campo de presiones, líneas de corriente y campo de velocidades muestran coherencia y concuerdan con los fenómenos esperados en la realidad. Por lo tanto, a primera vista, el algoritmo parece funcionar de acuerdo con las expectativas establecidas. Figura 3.13 Comparación de perfiles de velocidad horizontal, x=0.0175m. Figura 3.14 Comparación de perfiles de velocidad horizontal, x=0.025m. Figura 3.15 Comparación de perfiles de velocidad horizontal, x=0.0375m. 72 Capítulo 3. Resolución computacional de dominios fluidos (a) Campo de presiones (b) Campo de velocidades Figura 3.24 Campo de presiones y velocidades Re =100 . (a) Campo de presiones (b) Campo de velocidades Figura 3.25 Campo de presiones y velocidades Re =200 . (a) Campo de presiones (b) Campo de velocidades Figura 3.26 Campo de presiones y velocidades Re =500 . 3.2 Problema de cavidad impulsada 73 (a) Campo de presiones (b) Campo de velocidades Figura 3.27 Campo de presiones y velocidades Re =1000 . (a) Campo de presiones (b) Campo de velocidades Figura 3.28 Campo de presiones y velocidades Re =1500 . (a) Campo de presiones (b) Campo de velocidades Figura 3.29 Campo de presiones y velocidades Re =2000 . 74 Capítulo 3. Resolución computacional de dominios fluidos Revisando el campo de velocidades horizontales, muestran cómo las capas superiores se aceleran debido a los efectos difusivos. Para números de Reynolds bajos, estas capas parecen simétricas respecto a la sección central. Sin embargo, a medida que aumenta la velocidad de la tapa superior, el fluido del lado izquierdo no responde al cambio de velocidad hasta que está cerca de la pared superior, debido a las fuerzas inerciales dominantes sobre las viscosas. Del mismo modo, a mayor velocidad de la tapa, mayor es el efecto sobre las capas inferiores, que deben acelerar para mantener el caudal. Esto se refleja claramente en las figuras 3.24 - 3.29, donde la zona de velocidades negativas intensas, representada en color azul oscuro, se desplaza hacia abajo del dominio. Al representar los vectores de velocidad sobre el campo de velocidades horizontales, es posible identificar el centro de rotación de la recirculación. Debido al aumento de las velocidades en las capas superiores y su influencia en las inferiores, este centro de rotación se desplaza hacia abajo. Este fenómeno es más evidente al observar los perfiles de velocidad horizontal en una sección específica del dominio. En la figura 3.30, se aprecia cómo el punto de cruce de la curva con el valor nulo desciende con el aumento del número de Reynolds, comenzando en Re =100 por y=7.5 mm y llegando a descender para Re =2000 hasta y=5.5 mm, lo que indica un desplazamiento hacia abajo del centro de rotación. Además, en la parte superior de este mismo gráfico, se puede confirmar como, efectivamente, las capas de fluido se aceleran de forma exponencial pero solo a partir de y=9mm. Por debajo de esa altura, la velocidad aumenta, pero de forma lineal. En conclusión, las soluciones derivadas para el fenómeno de la cavidad impulsada reflejan una coherencia y consistencia notables, alineándose de manera congruente con las expectativas iniciales. No obstante, un paso crucial en la validación y contextualización de estos hallazgos radica en su comparación con los estudios previos existentes, lo cual será abordado en la siguiente sección de discusión de resultados. Este análisis comparativo no solo enriquecerá la comprensión de los resultados obtenidos, sino que también contribuirá a situar el presente estudio dentro del panorama académico actual. Figura 3.30 Comparación de perfiles de velocidad horizontal, x=0.05 m. 3.3 Flujo en canal con cambio de sección 75 3.3 Flujo en canal con cambio de sección 3.3.1 Descripción del problema Por último, se ha propuesto una variación del problema de canal de sección constante, en el que esta vez, se propone un flujo por un canal cuya sección se incrementa en forma de escalón, siendo un caso particularmente interesante dentro del estudio de la mecánica de fluidos. En este escenario, la sección transversal del canal experimenta un cambio abrupto en su área, creando un escalón o discontinuidad en la geometría del canal. Esto genera una serie de fenómenos hidrodinámicos que son importantes de comprender y analizar. En el caso de una sección de canal con área constante, el flujo tiende a ser más uniforme y predecible. Sin embargo, al introducir un escalón, se producen efectos como la separación del flujo, la formación de remolinos, cambios en la velocidad y la presión. Estos fenómenos pueden tener un impacto significativo en diversos aspectos, desde la eficiencia del transporte de fluidos hasta la seguridad de estructuras hidráulicas. En el diseño de perfiles aerodinámicos, como alas de aeronaves o cuerpos de cohetes, la comprensión de cómo el flujo reacciona ante cambios repentinos de geometría es esencial para optimizar la eficiencia y la estabilidad aerodinámica. Se concluye que, se trata de un problema un poco más complejo con una geometría relativamente sencilla y cuyo estudio es de utilidad para poner a prueba de nuevo el algoritmo SIMPLE y ver cómo se comporta en escenarios no uniformes. La geometría seleccionada consta de una entrada de altura H y una longitud de 2H . Tras esta entrada, la sección transversal se incrementa en forma de escalón, pasando a una altura total de 2H y una longitud de 8H hasta la salida del canal, siendo la longitud H=0.01m . Como se muestra en la figura 3.31, el flujo se considera que fluye desde la izquierda hasta la derecha, entrando con una velocidad uniforme uin y saliendo a una presión uniforme pout . Al igual que con el canal de sección constante, es de esperar que los perfiles de velocidad vayan cambiando debido a los efectos viscosos, pero en este caso, también habrá una zona de fuerte transición debido al cambio abrupto de la sección transversal. Continuando con la misma metodología, el análisis se hará en función del número de Reynolds. Figura 3.31 Flujo en canal con cambio de sección. 76 Capítulo 3. Resolución computacional de dominios fluidos 3.3.2 Condiciones de contorno La geometría planteada, consta de cuatro paredes sin deslizamiento que delimitan el canal por el que se mueve el flujo y donde el fluido en contacto tendrá velocidad nula. Nótese que las dos paredes adicionales, son las correspondientes al cambio de la sección transversal, siendo las zonas 4 y 6 iguales en cuanto a particularización de las ecuaciones. El fluido entra por la cara izquierda del conducto de altura H con una velocidad uniforme uin , mientras que sale por la derecha con una altura de 2H y una presión de salida uniforme pout . Al igual que en el canal de sección constante, se deben aplicar las mismas adaptaciones a las celdas que se encuentran cercanas a la salida y a la pared superior. En el dominio fluido se ha implementado una malla uniforme y co-ubicada, formada por una matriz de 40x200 celdas. Es decir, cada celda tiene una longitud de 0.5 mm. Si se revisa el código en el anexo C, es recomendable separar el dominio en dos zonas. Una primera zona cuya altura sea igual a la del escalón (for j=1: ny_in) y otra, desde dicha altura hasta la pared superior (for j=ny _in+1:ny). De esta forma, tenemos dos dominios diferenciados con dos longitudes diferentes y condiciones de contorno distintas. Lo recomendable es recorrer el dominio con bucles "For" diferenciando las zonas y sus condiciones de contorno. De forma detallada, se representa en la figura 3.32 las diferentes zonas descritas: Figura 3.32 Condiciones de contorno del canal con cambio de sección. •Zona 1: Entrada de velocidad uniforme •Zona 2: Pared superior sin deslizamiento •Zona 3: Salida de presión uniforme •Zona 4: Pared inferior sin deslizamiento •Zona 5: Pared izquierda sin deslizamiento •Zona 6: Pared inferior en la entrada sin deslizamiento •Zona 7: Celdas cercanas a pared superior •Zona 8: Celdas cercanas a la salida a presión uniforme. 3.3 Flujo en canal con cambio de sección 77 En referencia a las condiciones iniciales necesarias para iniciar el proceso iterativo, es crucial no comenzar desde cero. Iniciar con las variables en un estado nulo no resulta apropiado para este caso. Si se considera esta opción, durante la primera iteración, se experimentará un flujo de entrada con una velocidad uin que inducirá un desequilibrio significativo en la ecuación de continuidad. Este desequilibrio implicará correcciones de presión de gran magnitud, lo cual invalida las suposiciones de pequeñas correcciones efectuadas durante el desarrollo de las ecuaciones. Además, por la experiencia adquirida durante el estudio, se ha visto que esta ecuación es fundamental para la estabilidad del proceso iterativo. Por consiguiente, existe una alta probabilidad de que si las correcciones son demasiado grandes, la solución diverja y no se logre la convergencia esperada. Así pues, para evitar este gran gradiente adverso de presión al inicio de las iteraciones, es esencial imponer uin en todo el dominio. De esta manera, a medida que avancen las iteraciones, se corregirá gradualmente este efecto y se llegará a la convergencia de los resultados. Siguiendo el método empleado hasta ahora, se usará el número de Reynolds como parámetro. En este contexto, para este caso concreto, la longitud característica del problema es 2H=0.02m que es la altura total del conducto. De la misma forma, se define la densidad del fluido como ρ=1000 kg/m3 y su viscosidad como µ=2·10−5kg/(m·s) , de tal manera que el número de Reynolds pueda ser modificado fácilmente mediante la introducción del parámetro de velocidad de entrada uin . Esta velocidad será ajustada en un rango que va desde 0.1 m/s hasta 0.5 m/s, lo que resultará en un intervalo de valores de Re comprendido entre 100 y 500. Re =ρuin 2H µ(3.3) 3.3.3 Resultados y análisis Una vez más, los resultados que se analizarán abarcan el campo de presiones, el campo de velocidades y los perfiles de velocidad. En este último caso, se está interesado en analizar el efecto que tiene el aumento brusco de la sección. Se espera, que este cambio de sección genere una zona de recirculación del fluido que será tanto mayor, cuanto más turbulento sea el flujo, es decir, cuanto mayor sea el Reynolds. Además, debido a ese aumento en la sección, por conservación del caudal, las velocidades deberán disminuir y volver a desarrollarse hasta tener un perfil en forma de campana. Para verificar este fenómeno, se realiza un análisis detallado de los perfiles de velocidad horizontal en distintas ubicaciones a lo largo del canal, partiendo desde la zona del cambio de sección. Al igual que en los casos anteriores, también se monitorearán los residuos para confirmar la convergencia de la solución y ver las limitaciones que tiene el algoritmo. Las particularidades de este problema permitirán analizar cómo se comporta el algoritmo ante geometrías no uniformes. Tabla 3.3 Cambio de sección: iteraciones y residuos en función de Re. Re Iteraciones Residuo u Residuo v Residuo p 50 467 3,21E-11 1,06E-11 9,50E-11 100 346 3,69E-11 6,23E-12 9,66E-11 150 304 3,21E-11 1,37E-11 9,96E-11 200 1445 8,39E-11 1,80E-11 9,84E-11 300 1273 4,55E-11 1,91E-11 9,89E-11 78 Capítulo 3. Resolución computacional de dominios fluidos (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.33 Residuos y líneas de corriente Re =50 . (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.34 Residuos y líneas de corriente Re =100 . (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.35 Residuos y líneas de corriente Re =150 . 3.3 Flujo en canal con cambio de sección 79 (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.36 Residuos y líneas de corriente Re =200 . (a) Evolución de los residuos (b) Campo de presiones y líneas de corriente Figura 3.37 Residuos y líneas de corriente Re =300 . Este es un problema de mayor complejidad y, por ende, requiere un número considerablemente mayor de iteraciones para converger. A pesar de esto, se logra alcanzar valores dentro de la tolerancia establecida antes de llegar al límite de iteraciones globales. Es importante destacar que para los casos en los que Re =200 y Re =300 , fue necesario realizar ajustes en el código para lograr la convergencia. Estos ajustes implicaron extender la longitud del dominio hasta 20H y refinar el mallado, resultando en un sistema de ecuaciones más extenso, con una resolución de 80x800 nodos. Naturalmente, estos cambios incrementaron significativamente el tiempo de simulación de estos casos, pasando de unos pocos minutos a varias horas. En las gráficas de los residuos, se distinguen dos zonas principales: durante las primeras iteraciones, se realizan ajustes iniciales en las variables hasta alcanzar la convergencia, y estas curvas presentan una forma más irregular. Luego, sigue una fase donde los residuos comienzan a disminuir, pero muestran un patrón sinusoidal. Durante estas fluctuaciones, las correcciones de velocidad vertical y de presión parecen estar sincronizadas, alcanzando sus valores extremos en las mismas iteraciones, mientras que la corrección de velocidad horizontal muestra un comportamiento desfasado, con sus extremos ocurriendo en iteraciones opuestas. Estas fluctuaciones, junto con la tardanza de la 80 Capítulo 3. Resolución computacional de dominios fluidos ecuación de continuidad en converger, resultan en un número dispar de iteraciones entre los distintos casos, ya que dependen de la evolución de estos ciclos sinusoidales hacia la tolerancia establecida. Como Mazumder [ 13 ] señala, estas fluctuaciones son inherentes al proceso iterativo, ya que cada vez que se resuelve la cantidad de momento, se invalidan los resultados de la ecuación de continuidad y viceversa, cuando se corrige la presión, los resultados de la velocidad ya no son válidos y hay que resolverlos de nuevo. Al final, en última instancia, estos residuos serán tan pequeños, que apenas influirán en la solución final. Al examinar los campos de presión de los diferentes casos, se observa que en la región de entrada, antes del escalón, el comportamiento de la presión es similar al caso del canal de sección constante, donde la presión disminuye gradualmente hasta el escalón. Sin embargo, al llegar a la zona de cambio de sección, la presión vuelve a disminuir, aunque de manera no uniforme. Posteriormente, hay una zona de transición donde la presión se vuelve más uniforme, con un gradiente de presión bajo, hasta la salida del conducto. En el caso más laminar con Re =50 , la diferencia de presión entre la entrada y la salida es de ∆p=0.0092 Pa, mientras que para el flujo más turbulento con Re =300, la diferencia asciende a ∆p=0.1208 Pa. Revisando las líneas de corriente del dominio fluido, se aprecia que son predominantemente paralelas entre sí, indicando poca mezcla entre el fluido. Sin embargo, después del escalón, debido a las fuerzas inerciales, se forma una zona de recirculación de fluido que no se mezcla con el resto del flujo. Esta zona de recirculación es más extensa a medida que aumenta el número de Reynolds. Comparando las figuras 3.33 y 3.37 con Re =50 y Re =300 , respectivamente, se observa que en la primera, la zona llega hasta una longitud de x=0.035 m, mientras que para Re =300 , la zona se extiende hasta x=0.07 m. Dado que esta zona era considerablemente extensa en los últimos casos en comparación con la longitud original del dominio, fue necesario extender el canal para permitir el desarrollo del flujo y lograr la convergencia para Re =200 y Re =300 . Como el flujo se mueve de izquierda a derecha, la recirculación se produce en sentido de las agujas del reloj y aumenta en velocidad con el aumento de Re . En la zona intermedia, las líneas tienden a ser horizontales durante el desarrollo del flujo, y en la salida, parecen abrirse hacia las zonas de baja presión en las esquinas. Una vez más, al analizar los campos de velocidad obtenidos, se observa claramente cómo el flujo se desarrolla rápidamente en la entrada antes del escalón y cómo, después del escalón, se desvía hacia el nuevo centro de sección para luego comenzar a desarrollarse nuevamente. Naturalmente, debido a la conservación del caudal, la velocidad se reduce sustancialmente al aumentar la sección. Además, se aprecia claramente la zona de recirculación tras el escalón, donde la velocidad horizontal se vuelve negativa, y como su extensión es mayor a medida que aumenta Re. Para obtener una visión más detallada de la evolución de la velocidad en el canal, se han generado perfiles de velocidad en distintas secciones clave del dominio fluido: en la entrada, antes del escalón, junto al escalón, después del escalón, en una sección intermedia del dominio y en la salida. Estos perfiles revelan patrones significativos en el comportamiento del flujo. Inicialmente, la velocidad entra en el canal de manera uniforme desde la parte superior y se acelera rápidamente, siguiendo la curva de Poiseuille. Sin embargo, para el caso con Re =300 , no alcanza a completar su desarrollo antes del escalón, apreciándose una parte central del perfil, aún uniforme. Posterior al cambio de sección, la velocidad máxima disminuye y se desplaza hacia abajo, aunque este efecto es menos notable en casos con Re altos debido al dominio de las fuerzas inerciales. En esta etapa, en la parte inferior de la curva de velocidad, ya aparecen velocidades negativas, indicando la presencia de recirculación en esa región. 3.3 Flujo en canal con cambio de sección 81 (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.38 Campo de velocidades y perfiles, Re =50 . (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.39 Campo de velocidades y perfiles, Re =100 . (a) Campo de velocidades (b) Evolución de perfil de vel. horizontal Figura 3.40 Campo de velocidades y perfiles, Re =150 . 88 Capítulo 4. Discusión de los resultados En el análisis de estos resultados, se observan las mismas tendencias que se identificaron en el capítulo anterior. Sin embargo, hay una diferencia notable: ahora, con la disponibilidad de un canal de mayor longitud, se evidencia cómo la velocidad alcanza un máximo umax,300 =0.4486 m/s y se estabiliza en este valor. Por lo tanto, parece que finalmente se ha logrado desarrollar completamente el flujo. Esta conclusión se ve respaldada al comparar el perfil de velocidades en la salida con el perfil de Poiseuille, donde se observa una coincidencia efectiva, con un error relativo entre ellos de apenas el 0.31 % , y porque al representar el desarrollo de la velocidad central en la figura 4.6, se ve como tiende a una asíntota horizontal que marca la velocidad máxima. Al calcular nuevamente la longitud de entrada para Re =300 , como el 99 % de umax,300 =0.4486 m/s, se obtiene un valor de L=0.1284 m, que al adimensionalizarlo y compararlo con las curvas de Atkinson, Chen o Durst, muestra una convergencia mucho más cercana a dichas curvas, con errores significativamente reducidos en comparación con los obtenidos previamente. En resumen, los resultados alcanzados en las simulaciones son consistentes con los resultados analíticos y experimentales, demostrando una precisión satisfactoria para un estudio académico con objetivos didácticos centrados en el desarrollo de códigos propios para la resolución de dominios fluidos. Resalta la efectividad del enfoque metodológico empleado y su capacidad para generar resultados confiables y coherentes. Figura 4.4 Determinación de la longitud de entrada para Re =300. Tabla 4.3 Error relativo de la velocidad máxima para Re =300. Re Error relativo umax 300* 0,31% Tabla 4.4 Error relativo de la longitud de entrada para Re =300. Re Atkinson Chen Durst 300* 7,12% 3,13% 3,63% 4.1 Flujo en canal de sección constante 89 Figura 4.5 Comparación con perfil de Poiseuille para Re =300. Figura 4.6 Comparación de la longitud de entrada para Re =300. 90 Capítulo 4. Discusión de los resultados 4.2 Problema de cavidad impulsada Tal como se ha mencionado en capítulos anteriores, el problema de cavidad impulsada es usado como un problema clásico para la validación de códigos CFD por su geometría y condiciones de contorno simples, y lo que es más importante, ha sido fundamento de muchos y diversos estudios por lo que hay mucha literatura disponible para poder contrastar los resultados de este estudio. Los resultados obtenidos se van a comparar con los mostrados por Kundu et. al. en su libro Fluid Mechanics [ 10 ], donde resuelven este mismo problema aplicando el método MacCormack y lo comparan con los resultados alcanzados por Hou et al. con el método Lattice Boltzmann. El método MacCormack es un algoritmo numérico utilizado en la resolución de ecuaciones en derivadas parciales (EDP) especialmente en sistemas dependientes del tiempo y el espacio, lo que hace que sea un método muy empleado en la dinámica de fluidos computacional. Este método se basa en la división del dominio en una malla discreta, donde se calculan las soluciones en pasos de tiempo sucesivos. La idea clave del método MacCormack es dividir la resolución del problema en dos pasos: una predicción y una corrección. En el paso de predicción, consiste en la aproximación de las derivadas espaciales y temporales mediante la diferencia progresiva. Luego, en la etapa de corrección, se aproximan las derivadas espaciales mediante una diferencia regresiva. Combinando estas soluciones, logra obtener soluciones estables y precisas para una amplia gama de problemas. Como señala Carreño en su TFG [ 18 ], todo principiante en CFD comienza utilizando este tipo de métodos de diferencias finitas por la familiaridad con las ecuaciones diferenciales y después paulatinamente tienden a métodos basados en volúmenes finitos. Se trata de un método muy sencillo de programar para iniciarse en la mecánica de fluidos computacional. Por otro lado, se encuentra el método Lattice Boltzmann [ 1 ] que se trata de una técnica de simulación numérica. A diferencia de los métodos tradicionales basados en la resolución directa de ecuaciones de Navier-Stokes, Lattice Boltzmann simula el comportamiento del fluido a nivel microscópico, es decir, modela la distribución y comportamiento de partículas ficticias en una malla regular. En este método, las partículas se mueven en direcciones discretas, a través de los nodos de la malla dependiendo del conjunto de velocidades discretas. Cada partícula lleva un peso asociado, que representa su distribución de probabilidad. Durante cada paso de tiempo, las partículas interactúan localmente según reglas predefinidas, que conservan la masa y el momento. El algoritmo consta de dos pasos fundamentales: uno de colisión y otro de propagación. En el paso de colisión se redistribuyen las distribuciones de probabilidad de las partículas en cada nodo en función de sus interacciones. Posteriormente, en el paso de propagación, se mueven las distribuciones de probabilidad de acuerdo con las velocidades discretas, simulando el movimiento de las partículas. Algunas de sus ventajas, son su capacidad para simular flujos complejos en geometrías irregulares, su paralelización eficiente en sistemas computacionales modernos y su flexibilidad para modelar diferentes tipos de fluidos y fenómenos físicos. La evaluación de los resultados obtenidos en este estudio, a través de la aplicación del método SIMPLE, será llevada a cabo comparándolos con los enfoques de MacCormack y Lattice Boltzmann. Esta comparación se centrará en contrastar la ubicación del centro de rotación del fluido en la cavidad y en analizar los perfiles de velocidad en la sección intermedia para dos valores de Reynolds, Re =100 yRe =400. 4.2 Problema de cavidad impulsada 91 Al analizar minuciosamente la localización del centro de rotación, en general se observa que se encuentran muy próximos unos de otros. Esta proximidad es aún mayor para Re =100 en comparación con Re =400 , lo cual se atribuye al incremento de la turbulencia, y cómo es implementada por cada modelo. No obstante, el margen de error máximo se mantiene por debajo del 5 % , mostrando resultados positivos que se alinean con las expectativas y muestran una desviación mínima con respecto a los hallazgos reportados en la literatura. Sin embargo, es probable que estos resultados puedan mejorarse mediante el empleo de una malla más refinada, tal como se logró en el caso del flujo en canal de sección constante. Ubicación (X,Y) Lattice Boltzmann (0,6196 , 0,7373) MacComack (0,62 , 0,74) SIMPLE (0,6125 , 0,75) Figura 4.7 Comparación de la localización del centro de rotación para Re =100 . Ubicación (X,Y) Lattice Boltzmann (0,5608 , 0,6078) MacComack (0,57 , 0,61) SIMPLE (0,59 , 0,63) Figura 4.8 Comparación del campo de velocidades para Re =400 . Tabla 4.5 Errores relativos de la posición del centro de rotación. Re eLB eMC 100 1,33% 1,05% 400 4,436% 4,437% 92 Capítulo 4. Discusión de los resultados Inspirado por la figura 4.9 presentada por Kundu et al. [ 10 ] en su libro, donde compara los perfiles de velocidad adimensionalizados en la sección intermedia de la cavidad, se ha optado por usar esta figura de referencia para contrastar los resultados alcanzados con el método SIMPLE. Se han extrapolado los resultados y dividido en las figuras para Re =100 y para Re =400 , lo que permite evaluar con mayor detalle los valores logrados para los perfiles de velocidad. Figura 4.9 Comparación de resultados entre el método MacCormack y Hou et al [10]. En el caso de menor turbulencia, Re =100 , se aprecia que las curvas son bastante similares entre sí, reproduciendo casi con exactitud el perfil calculado por otros métodos. Reproduce con exactitud las condiciones de contorno impuestas tanto en la pared inferior donde la velocidad es nula, como la condición de arrastre de la tapa superior con una velocidad de Ulid . El perfil de velocidades trascurre casi paralelamente a la dibujada por Kundu. La mayor discrepancia parece encontrarse en el valor máximo de la velocidad negativa, aproximadamente en y≈0.5. Como era de esperar, por resultados previos, para el caso de Re =400 , las diferencias entre las curvas son mucho más notables, aunque el perfil general sigue una forma similar. Se siguen cumpliendo las condiciones de contorno impuestas por la pared inferior y el movimiento de la tapa superior. Comienzan desde la tapa inferior a tomar valores negativos rápidamente, para luego ir aumentando la velocidad de forma lineal, hasta cerca de la tapa superior. Ambas curvas muestran cómo el flujo cerca de la tapa superior, alrededor de y≈0.9 , cambia su comportamiento lineal para acelerarse hacia la velocidad nominal Ulid. El mayor error se localiza nuevamente en el valor máximo negativo de la velocidad que debería presentar el perfil. A pesar de estas diferencias, los resultados siguen siendo satisfactorios y se concluye que, a medida que aumenta la turbulencia, sería recomendable una mejor discretización del mallado para una captura más precisa. Una vez más, los resultados obtenidos a través del método SIMPLE siguen la tendencia esperada en comparación con los datos experimentales. Esto evidencia su consistencia y confiabilidad, destacando como un enfoque de implementación sencilla para cualquier académico y que logra resultados satisfactorios. 4.2 Problema de cavidad impulsada 93 Figura 4.10 Comparación de perfiles de velocidad para Re =100 . Figura 4.11 Comparación de perfiles de velocidad para Re =400 . 94 Capítulo 4. Discusión de los resultados 4.3 Flujo en canal con cambio de sección Por último, pero no menos importante, se procede a contrastar los resultados alcanzados para el problema de flujo en canal con cambio de sección. Al igual que el problema de cavidad impulsada, el estudio del flujo tras un cambio de sección, es un problema clásico de la mecánica de fluidos en los que existen numerosos estudios disponibles en la literatura, lo que permite usarse también como un problema para validar códigos CFD. De todos los estudios disponibles en la literatura, se ha seleccionado la publicación Numerical Investigation of Fluid Flow Through A 2D Backward Facing Step Channel[ 15 ] escrito por Hossain et al, donde compara los resultados experimentales medidos por Armaly et al, con los simulados mediante ANSYS Fluent. De esta forma, se podrá comparar los resultados tanto con datos experimentales como con datos CFD de un software comercial como es ANSYS Fluent. La geometría planteada por Hossain et al, es casi idéntica a la planteada en el problema de cambio de sección para este estudio. Se dispone de un canal con una sección de entrada de altura h , y tras recorrer una distancia 4h , el canal se ensancha hasta una altura de 2h hasta su salida. Para permitir el desarrollo completo del flujo tras el escalón, la longitud del canal hasta la salida se alarga hasta 20h . El planteamiento del problema es idéntico al de este trabajo, ir variando el valor de Reynolds y analizar el comportamiento del fluido en sus tres fases: laminar, transitorio y turbulento. El valor de Reynolds lo basa en el diámetro H=2h y como fluido de trabajo usa el agua, con lo que puede determinar el valor de la velocidad de entrada. Para resolver los casos, configura el programa con una solución basa en presiones y dentro de las opciones, selecciona el método SIMPLE de segundo orden upwind. Como método de convergencia usa la técnica estándar de discretización por mínimos cuadrados y un valor de convergencia de 10−5. Figura 4.12 Esquema de separación del flujo [15]. Al analizar las soluciones para este caso, tanto Hossain et al. como muchos otros autores, remarcan una particularidad de este caso, y es que a medida que aumenta la turbulencia en el flujo, se observa la generación de varias zonas de recirculación, donde el flujo se desprende de la pared y pasada una cierta distancia, se vuelve a reinyectar. Hossain et al. esquematiza este fenómeno en la figura 4.12. Inicialmente, para flujos laminares entre Re =100 a Re =300 , solo aparece una única 4.3 Flujo en canal con cambio de sección 95 zona de separación del flujo (I) que se encuentra junto al escalón. A medida que aumenta Re , en torno a Re =350 , se aprecia una segunda zona de separación del flujo (II) pero situada en la parte superior del canal. Si el flujo se vuelve más turbulento, es posible ver una tercera zona de separación (III), fenómeno que parece volverse recurrente cuanto más turbulento es el flujo de entrada, dando lugar a nuevas zonas de separación en la parte superior del canal e inferior de manera secuencial. A modo de validación de datos, Hossain et al. escoge comparar la longitud de reinyección de las dos primeras zonas de separación del flujo X1 , X4 y X5 . Contrasta estos datos con otros autores, entre ellos con Armaly et al. que proporciona datos medidos de forma experimental para flujos laminares y turbulentos, dando como resultado el gráfico mostrado en la figura 4.13. Figura 4.13 Validación de datos, Hossain et al. [15]. Volviendo a este estudio, es interesante usar esta misma metodología para validar los resultados obtenidos, ya que se compararán al mismo tiempo con datos experimentales y con datos numéricos calculados por el software comercial ANSYS Fluent que usa el mismo método que se ha tratado de implementar en Matlab. Sin embargo, la comparación solo se limitará al rango laminar por dos motivos: es el rango donde las hipótesis impuestas en el desarrollo siguen siendo validas y, porque, como comenta Hossain en el propio artículo y se corrobora en el propio gráfico, los datos a partir de Re =600 , empiezan a mostrar una gran desviación respecto a los datos experimentales. Llegados a este punto, para poder comparar resultados, se tuvieron que realizar más simulaciones con el código desarrollado, de esta manera se lograron obtener más casos con los que poder comparar. Se simularon tres casos más, hasta Re =500 puesto que, como se ha comentado anteriormente, Hossain explica que los datos a partir de Re =600 , mostraban desviaciones importantes. Los parámetros de simulación se adaptaron a cada caso según el Re para lograr un compromiso entre convergencia y tiempo de computación, seleccionando los siguiente parámetros: •Re =300:Tol. convergencia 10−11,n=40, sistema de ecuaciones de 80x800 nodos. •Re =400:Tol. convergencia 10−5,n=80, sistema de ecuaciones de 160x1600 nodos. •Re =500:Tol. convergencia 10−5,n=80, sistema de ecuaciones de 160x1600 nodos. A continuación se muestran los resultados obtenidos para estas simulaciones en las figuras 4.14, 4.15 y 4.16. 96 Capítulo 4. Discusión de los resultados (a) Evolución de los residuos (b) Evolución de perfil de vel. horizontal (c) Campo de presiones y líneas de corriente (d) Campo de velocidades Figura 4.14 Resultados para Re =350. 4.3 Flujo en canal con cambio de sección 97 (a) Evolución de los residuos (b) Evolución de perfil de vel. horizontal (c) Campo de presiones y líneas de corriente (d) Campo de velocidades Figura 4.15 Resultados para Re =400. 104 Capítulo 5. Conclusiones para obtener resultados confiables y consistentes en la simulación. La otra clave esencial para el funcionamiento del método SIMPLE aquí implementado, es la incorporación del método PWIM (Pressure-Weigthed Interpolation Method). Gracias a él, es posible emplear un mallado estructurado co-ubicado, ya que evita el problema de tablero de ajedrez de las oscilaciones de presión. Este problema de tablero de ajedrez es resuelto mediante PWIM porque tiene en cuenta para la interpolación de las velocidades en las caras, no solo las velocidades, si no además las presiones de los volúmenes adyacentes. Tal como se mostró en el desarrollo de las ecuaciones, el uso del mallado estructurado co-ubicado facilita enormemente el planteamiento y proyección de las ecuaciones en dichos volúmenes finitos, necesarias para la formación de los sistemas de ecuaciones que se deben resolver de forma iterativa. Se ha constatado que, como era evidente, la complejidad de la geometría que se desea simular afecta directamente la dificultad y el tiempo necesarios para aplicar las condiciones de contorno en los sistemas de ecuaciones. En estos casos, es necesario subdividir mucho más el dominio para recoger adecuadamente las condiciones de contorno. Además, se ha corroborado la relación directa entre la discretización del dominio y los tiempos de computación. Sin un modelo adecuado para flujos turbulentos, fue necesario aumentar la discretización del dominio, buscando un equilibrio entre la tolerancia para el criterio de convergencia, la discretización y los tiempos de simulación. Para tener una idea de los tiempos, en sistemas sencillos, los tiempos de simulación fueron de unos pocos minutos, mientras que en mallados de mayor complejidad, los tiempos de cálculo superaron las cuatro horas. Aunque se consideró la introducción de una ecuación de estado para relacionar la densidad y presión, finalmente se descartó esta opción, ya que requería incluir la ecuación de energía para calcular las temperaturas, puesto que los problemas con temperatura constante no son relevantes en aplicaciones prácticas. Por este motivo, se decidió dejarlo como una futura línea de estudio. El código implementado en este estudio basado en el método SIMPLE fue evaluado y validado para tres geometrías diferentes: canal de sección constante, cavidad impulsada y canal con cambio de sección, logrando buenos resultados en todos los casos, y que fueron contrastados con ejemplos presentes en la literatura. En el caso del canal de sección constante se estudió la evolución de los perfiles de velocidad horizontal a lo largo del canal para determinar la longitud de entrada en función del número de Reynolds, alcanzando resultados similares a los estimados por Atkinson, Chen y Durst. Respecto a la cavidad impulsada, se compararon la posición del centro de rotación dentro de la cavidad y los perfiles de velocidad horizontal en la mitad de la sección, logrando una precisión muy buena. Por último, para el canal con cambio de sección, se compararon las imágenes de las líneas de corriente frente a las obtenidas por simulación con ANSYS Fluent por Hossain. También se comparó la longitud de reinyección del flujo tras el escalón con los datos experimentales y CFD, donde quedaron claras las limitaciones del código implementado con el aumento de la turbulencia. Revisando el objetivo principal de este estudio, se considera que ha sido cumplido satisfactoriamente. El trabajo presentado constituye un manual académico adecuado para estudiantes o investigadores interesados en iniciarse en la dinámica de fluidos computacional (CFD). Resume los conceptos fundamentales necesarios para entender los desarrollos en este campo, ya que la literatura a menudo supone ser muy extensa y difusa. El estudio, explica detalladamente el método SIMPLE desde las ecuaciones de Navier-Stokes hasta el desarrollo del método iterativo y, aplica dicho método a ejemplos sencillos que permite familiarizarse con los sistemas de ecuaciones, condiciones de contorno y códigos implementados. Además, el Anexo B ofrece una guía práctica sobre las modificaciones necesarias en las ecuaciones según las condiciones de contorno más comunes, lo 105 que facilitará a los estudiantes e investigadores desarrollar sus propios códigos. En resumen, con este trabajo se ha logrado profundizar en la base de los programas CFD comerciales, para que dejen de ser empleados como una mera caja negra, y entender los procesos que siguen para resolver y simular la dinámica de los dominio fluidos. Con las hipótesis establecidas al inicio de este estudio, el rango de aplicación se limita a geometrías sencillas 2D, malladas con elementos rectangulares, para flujos laminares a bajo número de Reynolds y flujos incompresibles. Por este motivo, las futuras líneas de investigación deben ser orientadas a superar estos límites, priorizando según la necesidad de los dominios a resolver: • Extensión a dominios 3D: con la configuración actual y manteniendo un mallado ortogonal, bastaría con incorporar la tercera componente de velocidad w y la tercera ecuación de momento para un mallado hexaédrico. • Incorporación de un modelo de turbulencia: como por ejemplo el modelo k - ε , para aumentar el rango de Reynolds y poder simular flujos más turbulentos. Sería cuestión de incorporar el modelo de viscosidad cinemática modelada por k - ε y al final de cada iteración, tras la actualización de velocidades y presiones, actualizar también estos parámetros. • Simulación de flujos compresibles: bastaría con incluir la ecuación de energía y estado para modelar la variabilidad de la densidad. • Adaptación a mallados no estructurados: el objetivo sería poder usar mallados triangulares (2D) o tetraédricos (3D) para conseguir reproducir geometrías mucho más complejas. Esto conllevaría a volver a coger las ecuaciones generales y volver a proyectarlas según la orientación de las caras de cada triángulo, así como implementar una enumeración de los triángulos para posteriormente plantear los sistemas de ecuaciones que se deben resolver iterativamente. Estas líneas de investigación permitirán que el código implementado en este estudio, basado en el método SIMPLE, continúe evolucionando y ampliando sus capacidades hasta equipararse con las funcionalidades más avanzadas de los softwares comerciales de simulación de dinámica de fluidos computacional. Apéndice A Método Line-by-line El método line-by-line [ 2 ] es una técnica fundamental en el ámbito de los métodos computacionales, particularmente en la resolución de problemas que involucran ecuaciones diferenciales parciales (EDPs) y sistemas de ecuaciones lineales. Fue sugerido por Patankar en 1980. Este método está basado en la combinación de dos métodos: el algoritmo TDMA para problemas unidimensionales y el método iterarivo de Gauss-Seidel. Se caracteriza por descomponer el dominio del problema en una serie de líneas (o segmentos), donde se aplican algoritmos y técnicas numéricas para resolver el problema en cada una de ellas de manera secuencial. Estos segmentos se barren de forma horizontal o vertical en cada iteración, determinando el valor de todos los nodos de dicha columna o fila. Suponga que se encuentra en la iteración k , y se pretende resolver la fila j , donde las incógnitas son φk i,j(i=1,2...n) . El resto de variables en los nodos de la fila superior e inferior se consideran valores conocidos. En el caso de la fila inferior, los valores φk i,j−1(i=1,2...n) , acaban de ser calculados y actualizados en esta misma iteración k . En cuanto a los valores nodales de la fila superior φk−1 i,j+1(i=1,2...n) , fueron calculados en el barrido anterior, en la iteración k−1 . La relación entre todas estas variables nodales se puede expresar de forma genérica como la ecuación discretizada A.1. ai−1,jφk i−1,j+ai,jφk i,j+ai+1,jφk i+1,j=di,j−ai,j−1φk i,j−1+ai,j+1φk−1 i,j+1(A.1) Figura A.1 Mallado para resolución mediante el método line-by-line con barrido vertical. 107 108 Apéndice A. Método Line-by-line Ahora solo es cuestión de rescribir la ecuación para poder aplicar el método TDMA, es decir, se debe convertir en un sistema tridiagonal. Aplicando las relaciones mostradas a continuación, se observa como se puede construir una matriz A tridiagonal conformada por los factores que multiplican a las incógnitas, un vector d donde se encuentran las fuentes, y un vector x con las incógnitas de los valores nodales de la fila j. Ai,jφk i−1,j+Bi,jφk i,j+Ci,jφk i+1,j=Di,j(A.2) Ai,j=ai−1,ji=2,3,...n Bi,j=ai,ji=1,2,...n Ci,j=ai+1,ji=1,2,...n−1 Di,j=di,j−ai,j−1φk i,j−1+ai,j+1φk−1 i,j+1i=1,2,...n A x =d(A.3) A=         B1,jC1,j0 A2,jB2,jC2,j A3,jB3,j ... ......Cn−1,j 0An,jBn,j         ,x=       φk 1,j φk 2,j φk 3,j . . . φk n,j        ,d=       D1,j D2,j D3,j . . . Dn,j        Una vez resueltos y actualizados los valores nodales φk i,j(i=1,2...n) de la fila j , se repite este mismo proceso para las sucesivas filas hasta terminar el barrido. Tras esto, el método recomienda que se realice otro barrido, pero esta vez por columnas, consiguiendo así una mejor convergencia. Algunas ventajas que tiene este método son: • Eficiencia Computacional: Al resolver el problema en segmentos pequeños, el método lineby-line puede ser más eficiente computacionalmente que enfoques globales que requieren la discretización de todo el dominio simultáneamente. • Flexibilidad: El método line-by-line es altamente flexible y puede adaptarse a una amplia variedad de problemas, desde la simulación de fluidos hasta la propagación de ondas electromagnéticas. • Paralelización: La naturaleza secuencial del método line-by-line permite su fácil paralelización, lo que lo hace adecuado para implementaciones en sistemas de cómputo distribuido o con unidades de procesamiento gráfico (GPU). 109 Pero también, tiene algunas limitaciones y consideraciones: • Acoplamiento entre Segmentos: El acoplamiento entre segmentos puede introducir desafíos adicionales, especialmente en problemas no lineales o altamente no lineales, donde la convergencia puede ser más difícil de alcanzar. • Precisión: La precisión de la solución puede variar dependiendo de la densidad de la discretización y la naturaleza del problema. Es importante realizar análisis de sensibilidad y estudios de convergencia para garantizar la validez de los resultados obtenidos. Como se comentó en el capítulo 2, dentro de las funciones para la resolución de las ecuaciones de cantidad de movimiento y corrección de presiones, es necesario implementar la identificación de estas variables y coeficientes, de manera que se pueda formar la matriz tridiagonal A y el vector d . Una vez hecho esto, el sistema tridiagonal se resuelve con la función propia A\d implementada en Matlab. Apéndice B Condiciones de contorno El propósito de este Anexo B es proporcionar una guía detallada de las ecuaciones principales utilizadas en el algoritmo SIMPLE. En primer lugar, se presentan las ecuaciones generales necesarias para las celdas interiores. Luego, según la condición de contorno específica que se aplique, se indican los ajustes necesarios en las ecuaciones para implementar dicha condición en el código. Estos ajustes se resaltan en rojo para una fácil identificación. Para utilizar esta guía, primero se debe buscar la condición de contorno deseada y luego consultar las ecuaciones correspondientes que se desean resolver, ya sean: •Ecuación de cantidad de movimiento en X •Ecuación de cantidad de movimiento en Y •Interpolación PWIM para velocidades en las caras •Corrección de presiones (Ec. continuidad) •Actualización de velocidades Es importante tener en cuenta que, en el caso específico de las esquinas del dominio fluido, será necesario combinar adecuadamente los cambios provocados en las ecuaciones debido a la interacción de ambas condiciones de contorno. 111 112 Capítulo B. Condiciones de contorno B.1 Ecuaciones generales para celdas interiores Figura B.1 Celda interior de referencia. Ecuación de cantidad de movimiento en X AX OˆuO+AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS=SX O(B.1) AX O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.2) AX E=−|(ρu)e|−(ρu)e 2+µe ∆x∆y(B.3) AX W=−|(ρu)w|+(ρu)w 2+µw ∆x∆y(B.4) AX N=−|(ρv)n|−(ρv)n 2+µn ∆y∆x(B.5) AX S=−|(ρv)s|+(ρv)s 2+µs ∆y∆x(B.6) SX O=1 2(pW−pE)∆y(B.7) B.1 Ecuaciones generales para celdas interiores 113 Ecuación de cantidad de movimiento en Y AY OˆvO+AY WˆvW+AY EˆvE+AY NˆvN+AY SˆvS=SY O(B.8) AY O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.9) AY E=−|(ρu)e|−(ρu)e 2+µe ∆x∆y(B.10) AY W=−|(ρu)w|+(ρu)w 2+µw ∆x∆y(B.11) AY N=−|(ρv)n|−(ρv)n 2+µn ∆y∆x(B.12) AY S=−|(ρv)s|+(ρv)s 2+µs ∆y∆x(B.13) SY O=1 2(pS−pN)∆x(B.14) Interpolación PWIM para velocidades en las caras ˆue=ˆuO+ˆuE 2+1 2 pEE −pO 2AX OE ∆y+1 2 pE−pW 2AX OO ∆y−"1 AX OE +1 AX OO#pE−pO 2∆y(B.15) ˆvn=ˆvO+ˆvN 2+1 2 pNN −pO 2AY ON ∆y+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y(B.16) ˆuw|O=ˆue|W(B.17) ˆvs|O=ˆvn|S(B.18) 120 Capítulo B. Condiciones de contorno Aprox. 1ºOrden SY O=pO−pN 2∆x(B.69) Aprox. 2ºOrden SY O= (pO−pN)∆x(B.70) Interpolación PWIM para velocidades en las caras Aprox. 1ºOrden ˆvn=ˆvO+ˆvN 2+1 2 pNN −pO 2AY ON ∆y+∆x 2AY OOpN−pO 2−"1 AY ON +1 AY OO#pN−pO 2∆y (B.71) Aprox. 2ºOrden ˆvn=ˆvO+ˆvN 2+1 2 pNN −pO 2AY ON ∆y+∆x 2AY OO [pN−pO]−"1 AY ON +1 AY OO#pN−pO 2∆y(B.72) ˆvs|O=0(B.73) Corrección de presiones (Ec. continuidad) Ap O=ρe∆y2 2"1 AX OE +1 AX OO#+ρw∆y2 2"1 AX OW +1 AX OO#+ +ρn∆x2 2"1 AY ON +1 AY OO#(B.74) Ap S=0(B.75) Sp O=−[(ρeˆue−ρwˆuw)∆y+ (ρnˆvn)∆x](B.76) Actualización de velocidades vs|O=0(B.77) B.5 Pared superior sin deslizamiento 121 B.5 Pared superior sin deslizamiento Figura B.6 Pared superior sin deslizamiento. Ecuación de cantidad de movimiento en X AX O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +3µO ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.78) AX N=0(B.79) AX S=−|(ρv)s|+(ρv)s 2+µs ∆y+µO 3∆y∆x(B.80) Ecuación de cantidad de movimiento en Y AY O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +3µO ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.81) AY N=0(B.82) AY S=−|(ρv)s|+(ρv)s 2+µs ∆y+µO 3∆y∆x(B.83) 122 Capítulo B. Condiciones de contorno SY O=1 2(pS−pN)∆x(B.84) Aprox. 1ºOrden SY O=pS−pO 2∆x(B.85) Aprox. 2ºOrden SY O= (pS−pO)∆x(B.86) Interpolación PWIM para velocidades en las caras ˆvn=0(B.87) Para la fila j=m−1: Figura B.7 Celda cerca de la pared superior. Aprox. 1ºOrden ˆvn=ˆvO+ˆvN 2+∆x 2AY ONpN−pO 2+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y(B.88) Aprox. 2ºOrden ˆvn=ˆvO+ˆvN 2+∆x 2AY ON [pN−pO]+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y(B.89) B.5 Pared superior sin deslizamiento 123 Corrección de presiones (Ec. continuidad) Ap O=ρe∆y2 2"1 AX OE +1 AX OO#+ρw∆y2 2"1 AX OW +1 AX OO#+ +ρs∆x2 2"1 AY OS +1 AY OO#(B.90) Ap N=0(B.91) Sp O=−[(ρeˆue−ρwˆuw)∆y+ (−ρsˆvs)∆x](B.92) Actualización de presiones vk+1 n=0(B.93) 124 Capítulo B. Condiciones de contorno B.6 Pared superior con deslizamiento Figura B.8 Pared superior con deslizamiento. Ecuación de cantidad de movimiento en X AX OˆuO+AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS=SX O+µO 8∆x 3∆yulid (B.94) AX O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +3µO ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.95) AX N=0(B.96) AX S=−|(ρv)s|+(ρv)s 2+µs ∆y+µO 3∆y∆x(B.97) Ecuación de cantidad de movimiento en Y AY O=|(ρu)e|+(ρu)e 2+µe ∆x+|(ρu)w|−(ρu)w 2+µw ∆x∆y+ +3µO ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.98) AY N=0(B.99) B.6 Pared superior con deslizamiento 125 AY S=−|(ρv)s|+(ρv)s 2+µs ∆y+µO 3∆y∆x(B.100) SY O=1 2(pS−pN)∆x(B.101) Aprox. 1ºOrden SY O=pS−pO 2∆x(B.102) Aprox. 2ºOrden SY O= (pS−pO)∆x(B.103) Interpolación PWIM para velocidades en las caras ˆvn=0(B.104) Para la fila j=m−1: Figura B.9 Celda cerca de la pared superior con deslizamiento. Aprox. 1ºOrden ˆvn=ˆvO+ˆvN 2+∆x 2AY ONpN−pO 2+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y (B.105) Aprox. 2ºOrden ˆvn=ˆvO+ˆvN 2+∆x 2AY ON [pN−pO]+1 2 pN−pS 2AY OO ∆y−"1 AY ON +1 AY OO#pN−pO 2∆y(B.106) 126 Capítulo B. Condiciones de contorno Corrección de presiones (Ec. continuidad) Ap O=ρe∆y2 2"1 AX OE +1 AX OO#+ρw∆y2 2"1 AX OW +1 AX OO#+ +ρs∆x2 2"1 AY OS +1 AY OO#(B.107) Ap N=0(B.108) Sp O=−[(ρeˆue−ρwˆuw)∆y+ (−ρsˆvs)∆x](B.109) Actualización de presiones vk+1 n=0(B.110) B.7 Entrada a velocidad constante 127 B.7 Entrada a velocidad constante Figura B.10 Entrada del fluido a velocidad constante. Ecuación de cantidad de movimiento en X AX OˆuO+AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS=SX O+ρin uin uin +8µin 3∆xuin∆y(B.111) AX O=|(ρu)e|+(ρu)e 2+µe ∆x+3µin ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.112) AX E=−|(ρu)e|−(ρu)e 2+µe ∆x+µin 3∆x∆y(B.113) AX W=0(B.114) Aprox. 1ºOrden SX O=pO−pE 2∆y(B.115) Aprox. 2ºOrden SX O= (pO−pE)∆y(B.116) Ecuación de cantidad de movimiento en Y AY O=|(ρu)e|+(ρu)e 2+µe ∆x+3µin ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.117) AY E=−|(ρu)e|−(ρu)e 2+µe ∆x+µin 3∆x∆y(B.118) AY W=0(B.119) 128 Capítulo B. Condiciones de contorno Interpolación PWIM para velocidades en las caras Aprox. 1ºOrden ˆue=ˆuO+ˆuE 2+1 2 pEE −pO 2AX OE ∆y+∆y 2AX OOpE−pO 2−"1 AX OE +1 AX OO#pE−pO 2∆y (B.120) Aprox. 2ºOrden ˆue=ˆuO+ˆuE 2+1 2 pEE −pO 2AX OE ∆y+∆y 2AX OO [pE−pO]−"1 AX OE +1 AX OO#pE−pO 2∆y (B.121) ˆuw|O=uin (B.122) Corrección de presiones (Ec. continuidad) Ap O=ρe∆y2 2"1 AX OE +1 AX OO#+ +ρn∆x2 2"1 AY ON +1 AY OO#+ρs∆x2 2"1 AY OS +1 AY OO#(B.123) Ap W=0(B.124) Sp O=−[(ρeˆue−ρinuin)∆y+ (ρnˆvn−ρsˆvs)∆x](B.125) Actualización de velocidades uw|O=uin (B.126) B.8 Entrada a presión constante 129 B.8 Entrada a presión constante Figura B.11 Entrada del fluido a presión constante. Ecuación de cantidad de movimiento en X Aprox. 1ºOrden uin =uk O(B.127) Aprox. 2ºOrden uin =3 2uk O−1 2uk E(B.128) AX OˆuO+AX WˆuW+AX EˆuE+AX NˆuN+AX SˆuS=SX O+ρin uin uin +8µin 3∆xuin∆y(B.129) AX O=|(ρu)e|+(ρu)e 2+µe ∆x+3µin ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.130) AX E=−|(ρu)e|−(ρu)e 2+µe ∆x+µin 3∆x∆y(B.131) AX W=0(B.132) SX O=pin −pO+pE 2∆y(B.133) Ecuación de cantidad de movimiento en Y AY O=|(ρu)e|+(ρu)e 2+µe ∆x+3µin ∆x∆y+ +|(ρv)n|+(ρv)n 2+µn ∆y+|(ρv)s|−(ρv)s 2+µs ∆y∆x (B.134)