Problemas de convección-difusión térmica mediante FVM en geometrías variables con mallas estructuradas y no estructuradas
Abstract
En el ámbito de la simulación numérica, el Método de Volúmenes Finitos (FVM) se erige como una de las principales herramientas para abordar problemas en la mecánica de fluidos mediante programas de Dinámica de Fluidos Computacional (CFD). Este Trabajo de Fin de Grado se propone extender la formulación de este método, que Alejandro Cumplido desarrolló con notable éxito en geometrías rectangulares, hacia geometrías variables. El objetivo es establecer las bases de la formulación del FVM en mallas estructuradas mediante un cambio de variables, así como en mallas no estructuradas. El proyecto incluye la elaboración de códigos en el lenguaje de programación MATLAB para resolver problemas de convección-difusión térmica en la mecánica de fluidos. Estos problemas requieren a su vez calcular el campo de velocidades, abordado también en este trabajo mediante la resolución de la ecuación de continuidad para flujos incompresibles según la teoría ideal. En esencia, este proyecto persigue como objetivo exponer de manera sencilla el funcionamiento interno de los programas CFD, los cuales son a menudo utilizados como cajas negras por el usuario final
Full text
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 Grado Grado en Ingeniería Aeroespacial. Mención en Vehículos Aeroespaciales Problemas de convección-difusión térmica mediante FVM en geometrías variables con mallas estructuradas y no estructuradas. Autor: Pedro Vico Cantero 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 Grado Grado en Ingeniería Aeroespacial. Mención en Vehículos Aeroespaciales Problemas de convección-difusión térmica mediante FVM en geometrías variables con mallas estructuradas y no estructuradas. Autor: Pedro Vico Cantero 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 Grado: Problemas de convección-difusión térmica mediante FVM en geometrías variables con mallas estructuradas y no estructuradas. Autor: Pedro Vico Cantero 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 n primer lugar, deseo expresar mi más profundo agradecimiento al Dr. Miguel Pérez-Saborid Sánchez- Pastor por brindarme la oportunidad de trabajar en este proyecto. Un tema que, al principio, me generaba cierta inquietud y, al mismo tiempo, despertaba una gran curiosidad, se reveló con el tiempo como una fuente de gran fascinación. Agradezco también el haber puesto a mi disposición su valiosa experiencia y vasto conocimiento, así como su constante atención y dedicación. Estos elementos no solo han guiado el proyecto de manera ejemplar, sino que han contribuido significativamente al desarrollo de mi conocimiento y a mi crecimiento personal. En segundo lugar, extiendo mi gratitud a todos mis familiares y amigos, cuyos constantes ánimos y apoyo incondicional han sido el cimiento sobre el cual se ha construido este proceso. Sin su ayuda, no habría sido posible llegar hasta donde me encuentro hoy. En especial, quiero agradecer a mis hermanas, por ser un ejemplo a seguir y una fuente de inspiración constante, y a mis padres, por su paciencia y por permitir que mis sueños estén cada vez más cerca de convertirse en realidad. Por último, pero no menos importante, quiero expresar mi más sincero reconocimiento a todos los profesores que me han acompañado a lo largo de mi educación, quienes con su dedicación y sabiduría han moldeado mi camino. Agradezco también a los investigadores, cuyo incansable afán de conocimiento y descubrimiento sigue abriendo nuevas fronteras del saber. Y, en general, a todas aquellas personas curiosas que se dedican a desvelar los misterios que tanto nos asombran y que, con su esfuerzo y pasión, contribuyen al progreso de la humanidad. Su labor es una fuente indiscutible de inspiración para nuestro avance colectivo. Pedro Vico Cantero Escuela Técnica Superior de Ingeniería Sevilla, 2024 I
Resumen E n el ámbito de la simulación numérica, el Método de Volúmenes Finitos (FVM) se erige como una de las principales herramientas para abordar problemas en la mecánica de fluidos mediante programas de Dinámica de Fluidos Computacional (CFD). Este Trabajo de Fin de Grado se propone extender la formulación de este método, que Alejandro Cumplido desarrolló con notable éxito en geometrías rectangulares, hacia geometrías variables. El objetivo es establecer las bases de la formulación del FVM en mallas estructuradas mediante un cambio de variables, así como en mallas no estructuradas. El proyecto incluye la elaboración de códigos en el lenguaje de programación MATLAB para resolver problemas de convección-difusión térmica en la mecánica de fluidos. Estos problemas requieren a su vez calcular el campo de velocidades, abordado también en este trabajo mediante la resolución de la ecuación de continuidad para flujos incompresibles según la teoría ideal. En esencia, este proyecto persigue como objetivo exponer de manera sencilla el funcionamiento interno de los programas CFD, los cuales son a menudo utilizados como cajas negras por el usuario final. III
Índice de Figuras 1.1 Mallado estructurado de un perfil aerodinámico (imagen tomada de [1]) 4 1.2 Mallado no estructurado de un perfil aerodinámico (imagen tomada de [1]) 5 2.1 Mallado de una geometría 12 2.2 Plano computacional 15 2.3 Flujo y vector normal en la frontera 17 2.4 Mallado de la geometría del problema estudiado 18 2.5 Posición de los coeficientes en la matriz dispersa 20 2.6 Representación del problema descrito 21 2.7 Mallado de una geometría 22 2.8 Comparación de soluciones 22 2.9 Comparación de soluciones curvas de nivel 23 2.10 Caudal a través de una línea 24 2.11 Conducto bidimensional 25 2.12 Campo de velocidades a través de un conducto exponencial 28 2.13 Líneas de corriente a través de un conducto exponencial 29 2.14 Función de corriente en un conducto exponencial 29 2.15 Comparación de soluciones para v=10−6m/s31 2.16 Comparación de soluciones para v=10−4m/s31 2.17 Comparación de soluciones para v=10−2m/s32 2.18 Comparación de velocidades para v=10−4m/s32 3.1 Mallado estructurado por bloques para un motor 35 3.2 Celda típica de un mallado no estructurado 36 3.3 Mallado no estructurado de un perfil 36 3.4 Indexación global de una malla no estructurada 40 3.5 55 3.6 Comparación de soluciones 55 3.7 Comparación de soluciones curvas de nivel 55 3.8 Mallado generado por Fluent 56 3.9 Comparación de soluciones para dos mallados no estructurados 58 3.10 Mallado no estructurado generado por Fluent 62 3.11 Campo de velocidades en una malla no estructurada 62 3.12 Función de corriente 62 3.13 Comparación de soluciones para ve=10−6m/s64 3.14 Comparación de soluciones para ve=10−4m/s65 3.15 Comparación de soluciones para ve=10−2m/s65 A.1 Entorno de trabajo de ANSYS 69 A.2 70 A.3 71 XI
XII Índice de Figuras A.4 Creación de una superficie a través de bordes 72 A.5 72 A.6 73 A.7 73 A.8 74 A.9 74 A.10 75 A.11 76
Índice de Tablas 3.1 Métodos del objeto triangulation de MATLAB 56 XIII
Índice Abreviado Resumen III Abstract V Notación VII Índice de Figuras XI Índice de Tablas XIII Índice Abreviado XV 1 Introducción 1 1.1 Fundamentos del CFD 1 1.2 Funcionamiento de los códigos CFD 3 1.3 Historia del CFD 6 1.4 Ecuaciones de la mecánica de fluidos 7 2 Mallado Estructurado en coordenadas curvilíneas 11 2.1 FVM en coordenadas curvilíneas 11 2.2 Problema de difusión pura, estacionario y bidimensional 14 2.3 Problema de convección-difusión térmica 23 3 FVM para problemas de convección-difusión térmica. Mallado no estructurado 35 3.1 Introducción 35 3.2 Problema difusivo, estacionario, bidimensional 37 3.3 Problema de convección-difusión térmica estacionario y bidimensional 58 4 Conclusiones y trabajos futuros 67 Apéndice A Solución problema difusivo con Fluent 69 Apéndice B Códigos de MATLAB 77 B.1 Mallado estructurado 77 B.2 Mallado no estructurado 89 Bibliografía 119 XV
Índice Resumen III Abstract V Notación VII Índice de Figuras XI Índice de Tablas XIII Índice Abreviado XV 1 Introducción 1 1.1 Fundamentos del CFD 1 1.1.1 Método de las Diferencias Finitas 2 1.1.2 Método de los Volúmenes Finitos 2 1.1.3 Método de los Elementos Finitos 3 1.2 Funcionamiento de los códigos CFD 3 1.2.1 Preprocesador 3 1.2.2 Solver 5 1.2.3 Post-procesador 6 1.3 Historia del CFD 6 1.4 Ecuaciones de la mecánica de fluidos 7 2 Mallado Estructurado en coordenadas curvilíneas 11 2.1 FVM en coordenadas curvilíneas 11 2.2 Problema de difusión pura, estacionario y bidimensional 14 2.2.1 Mallado del dominio 14 2.2.2 Discretización de las ecuaciones 15 2.2.3 Aplicación práctica: código de MATLAB 17 2.3 Problema de convección-difusión térmica 23 2.3.1 Función de corriente para flujo estacionario e incompresible 24 2.3.2 Aproximación del término convectivo 29 2.3.3 Aplicación práctica 31 3 FVM para problemas de convección-difusión térmica. Mallado no estructurado 35 3.1 Introducción 35 3.2 Problema difusivo, estacionario, bidimensional 37 3.2.1 Discretización de las ecuaciones. Flujos normales y tangenciales 37 3.2.2 Estructuración del Dominio: Gestión y Almacenamiento de Datos Geométricos 39 3.2.3 Estructuración de la ecuación en forma matricial 49 3.2.4 Aplicación práctica. Problema de difusión 54 3.2.5 Aplicación práctica 2D. Malla no estructurada triangular 56 3.3 Problema de convección-difusión térmica estacionario y bidimensional 58 4 Conclusiones y trabajos futuros 67 XVII
XVIII Índice Apéndice A Solución problema difusivo con Fluent 69 Apéndice B Códigos de MATLAB 77 B.1 Mallado estructurado 77 B.1.1 Script configuración problema difusivo y convectivo-difusivo 77 B.1.2 Función que resuelve el campo potencial 79 B.1.3 Problema convectivo-difusivo 84 B.2 Mallado no estructurado 89 B.2.1 Problema difusivo con mallado cuadrangular 89 B.2.2 Solución función de corriente 99 B.2.3 Problema convectivo-difusivo con una malla externa 109 Bibliografía 119
1 Introducción L a Dinámica de Fluidos Computacional (CFD, por sus siglas en inglés) se ha consolidado como una herramienta indispensable en diversas áreas de la ingeniería y la investigación científica. La capacidad de simular y analizar fenómenos complejos de flujo de fluidos con alta precisión ha revolucionado industrias como la automotriz, la aeroespacial, la ambiental y muchas otras. La razón de ser de los CFD radica en su habilidad para proporcionar soluciones detalladas y cuantitativas a problemas que serían extremadamente difíciles o costosos de abordar experimentalmente. El CFD se basa en la resolución numérica de las ecuaciones de Navier-Stokes, que describen el comportamiento de los fluidos en movimiento. Este campo ha avanzado significativamente gracias al desarrollo de métodos numéricos robustos y al incremento en la capacidad de cómputo de los ordenadores modernos. En este capítulo se expondrán los fundamentos que sustentan los códigos de Dinámica de Fluidos Computacional (CFD). Se ofrecerá una breve reseña histórica de su desarrollo y evolución, destacando los hitos más significativos en el campo. Finalmente, se presentarán las ecuaciones de la mecánica de fluidos que constituyen el núcleo teórico de las simulaciones CFD, proporcionando un marco completo para comprender cómo estas ecuaciones se aplican en la modelización y análisis de fenómenos de flujo. 1.1 Fundamentos del CFD En el análisis de modelización matemática de fenómenos físicos, existen dos categorías principales de modelos matemáticos utilizados: los empíricos y los basados en física fundamental. Los modelos empíricos se desarrollan a partir de datos experimentales y experiencia práctica, generando ecuaciones que describen relaciones causales específicas. Por ejemplo, la obtención de una correlación para calcular coeficientes de transferencia de calor como función de las propiedades termofísicas y la velocidad del flujo. Aunque son fáciles de entender y emplear, su rango de aplicabilidad siempre está limitado por las condiciones específicas del experimento realizado. Por el contrario, los modelos basados en física fundamental, como las leyes de conservación, tienen un amplio alcance de aplicación. Aunque más complejos de resolver, ofrecen una descripción más universal de los fenómenos físicos. La Dinámica de Fluidos Computacional (CFD) se basa en la aplicación de métodos numéricos para resolver las ecuaciones fundamentales basadas en principios físico-matemáticos que gobiernan el comportamiento de los fluidos, principalmente las ecuaciones de Navier-Stokes. Estas ecuaciones diferenciales en derivadas parciales (PDE por sus siglas en inglés) describen cómo cambian la velocidad, la presión, la temperatura y otras propiedades de un fluido en el espacio y el tiempo. Debido a su naturaleza compleja y no lineal, la solución analítica directa de las ecuaciones de Navier-Stokes es generalmente imposible, especialmente para problemas prácticos que implican geometrías complicadas. De este modo, los códigos buscan resolver numéricamente estas ecuaciones, permitiendo que los resultados se expresen mediante gráficos mucho más representativos. Estos gráficos son más accesibles y comprensibles que un conjunto de ecuaciones, cuyo entendimiento está limitado a aquellos con un conocimiento avanzado y una capacidad interpretativa elevada en la naturaleza físico-matemática del problema. 1
8Capítulo 1. Introducción y cada una de las partículas además de conocer sus condiciones iniciales de posición y velocidad resultaría prohibitivamente costoso. La alta densidad de partículas permite realizar la aproximación del fluido como medio continuo, lo que implica considerar que incluso para un elemento infinitamente pequeño (en el sentido del cálculo infinitesimal) existe un número lo suficientemente grande de partículas para estimar un promedio de velocidad, presión, temperatura, densidad, etc. Las ecuaciones que gobiernan el movimiento de los fluidos se derivan de que el comportamiento dinámico de estos se basa en leyes de conservación: •Conservación de la masa: La masa del fluido se conserva, conforme al principio de continuidad. • Conservación de la cantidad de movimiento: Esta ecuación se deriva de la Segunda Ley de Newton, la cual establece que la tasa de cambio de la cantidad de movimiento de un fluido es igual a la suma de las fuerzas actuantes sobre él. • Conservación de la energía: La ecuación de conservación de la energía establece que la diferencia entre el flujo de energía saliente y el flujo entrante, más cualquier generación interna de energía, es igual a cero. Esta ecuación se basa en la Primera Ley de la Termodinámica. El conjunto de ecuaciones que se derivan de las leyes de conservación describe el comportamiento de los fluidos en diversas condiciones, siendo fundamentales en la ingeniería y las ciencias aplicadas. Estas ecuaciones, conocidas como las ecuaciones de Navier-Stokes, permiten predecir el movimiento de los fluidos y su interacción con los entornos circundantes, abarcando fenómenos complejos como el flujo de fluidos viscosos. Su validez se extiende a cualquier fluido newtoniano compresible, siempre y cuando se lo considere como continuo. Ecuación de Continuidad (Conservación de la Masa): ∂ρ ∂t+∇·(ρv) = 0(1.1) Ecuación de Cantidad de Movimiento (Conservación de la Cantidad de Movimiento): ∂ρv ∂t+v·∇(ρv) = −∇p+∇·(µ∇v)+ρfm(1.2) Ecuación de Energía (Conservación de la Energía): ρ∂e ∂t+ρv·∇e=−p∇·v+τ′:∇v+Qr+Qq+∇·(k∇T)(1.3) Para cerrar el sistema de ecuaciones, se utilizan ecuaciones de estado que relacionan las variables termodinámicas, como la presión, la densidad y la temperatura. p=p(ρ,T),e=e(ρ,T)(1.4) Gas ideal: p=ρRgT,e=CvT(1.5) En las ecuaciones de Navier-Stokes, t representa el tiempo, fm es el vector de fuerzas másicas por unidad de volumen, τ′:∇v es el tensor de tensiones viscosas, Qr y Qq corresponden a los flujos de radiación y reacción química, respectivamente. Rg=R/Mm=Cp−Cv=Cv(γ−1) incluye R como la constante universal de los
1.4 Ecuaciones de la mecánica de fluidos 9 gases, Mm como la masa molecular del gas, Cp y Cv como las capacidades caloríficas a presión constante y volumen constante, respectivamente, y γcomo el coeficiente de dilatación adiabática. Estas ecuaciones, que conforman un conjunto de siete ecuaciones con siete incógnitas ( ρ , v , T y P ), permiten modelar y predecir el comportamiento de los fluidos en una amplia variedad de aplicaciones, desde el flujo de aire alrededor de un avión hasta el movimiento de los océanos, aplicando adecuadas condiciones de contorno en las superficies que delimitan el dominio del fluido. La similitud que presentan estas ecuaciones entre sí, permite expresar de forma genérica en una única ecuación, como la evolución de la propiedad φ . De este modo, se llega a la expresión conocida como ecuación de transporte para la propiedad φ: ∂φ ∂t+∇·(ρφv) = ∇·(Γ∇φ)+Sφ(1.6) Esta ecuación consta, por un lado, de un término de variación temporal y otro convectivo, que trasporta la propiedad φ según la velocidad v , y al otro lado, del término difusivo y de generación de calor, respectivamente. Se trata de una ecuación diferencial, la cual supone el punto de partida para los algoritmos computacionales del Finite Volume Method (FVM). El paso característico que da nombre a este método consiste en integrar dicha ecuación en un volumen de control apropiado: ZVC ∂(ρφ) ∂tdV +ZVC ∇·(ρφv)dV =ZVC ∇·(Γ∇φ)dV +ZVC SφdV (1.7) Las integrales de volumen de la componente convectiva a la izquierda de la ecuación y de la componente difusiva a la derecha pueden reescribirse como integrales alrededor de una superficie haciendo uso del Teorema de la Divergencia de Gauss: ∂ ∂tZVC ρφdV+ZA n·(ρφv)dA =ZA n·(Γ∇φ)dA +ZVC SφdV (1.8) Esta nueva expresión representa el balance de flujo de la propiedad φ para un volumen de control. El primer término describe la tasa de cambio de φ en el volumen de control. El producto escalar n·(ρφv) expresa la componente del flujo de φ que sale en dirección normal a la superficie que delimita el volumen de control. De este modo, esta expresión representa la tasa de disminución neta de φ debido a la convección. El flujo difusivo, en su esencia, se manifiesta como positivo cuando se desplaza en la dirección opuesta al gradiente de φ . Así, un valor positivo en la expresión n·(Γ∇φ) = −n·(−Γ∇φ) refleja flujo de calor en la dirección de la normal interior. Puede, por tanto, interpretarse esta formulación como la componente del flujo difusivo positivo en dirección normal interior a la superficie que encierra al volumen de control. El último término se refiere al aumento de φ debido a las fuentes que generan esta propiedad en el interior del fluido. La integración de la ecuación de transporte 1.6 representa finalmente una ecuación de balance de la propiedad φ para una región finita del dominio (a nivel macroscópico), que no es más que la conservación de esta propiedad en un volumen de control. En esencia, el Método de Volúmenes Finitos (FVM) consiste en resolver la ecuación de balance integral en cada uno de los volúmenes de control en los que se ha discretizado una geometría, considerando unas condiciones de contorno específicas en su frontera. Estos volúmenes de control, o celdas, se componen generalmente de líneas rectas en el caso bidimensional y planos en el espacio tridimensional. En cada celda se formula la ecuación de balance integral utilizando los valores de la propiedad φ en las celdas colindantes, lo que requiere una nomenclatura representativa para esquematizar claramente su formulación. Este nombramiento es especialmente crucial en mallas estructuradas, donde todas las celdas vecinas se involucran en la ecuación de la celda en la que se resuelve la ecuación de balance integral, etiquetada como O. Esto se detallará en profundidad en la siguiente sección, correspondiente a la formulación del FVM para mallas estructuradas.
10 Capítulo 1. Introducción En definitiva, los programas de CFD permiten modelar el comportamiento de fluidos bajo diversas condiciones, optimizar diseños y predecir el rendimiento de sistemas antes de su implementación. Sin embargo, a pesar de su amplio uso y los beneficios que ofrecen, los programas de CFD suelen ser considerados como una caja negra para muchos usuarios. Esto se debe a que, en la mayoría de los casos, los algoritmos y métodos numéricos subyacentes son desconocidos para los usuarios finales. Incluso en el caso de programas de código abierto, la complejidad y la falta de documentación clara pueden hacer que entender su funcionamiento interno sea una tarea ardua y desalentadora. Este trabajo tiene como objetivo ofrecer un entendimiento profundo de cómo funcionan estos programas desde adentro, con un enfoque particular en el estudio de los mallados estructurados y no estructurados. Los mallados son un componente crucial en cualquier simulación de CFD, ya que definen cómo se discretiza el dominio de cálculo y, por ende, influyen directamente en la precisión y eficiencia de las soluciones numéricas. Es imprescindible destacar el papel fundamental que desempeñan las obras de Alejandro Cumplido [ 4 ], Mazumder [ 8 ] y Versteeg [ 3 ] como fuentes de referencia e influencia en el desarrollo de este trabajo. Sin embargo, la aportación singular de este TFG reside en los programas de MATLAB desarrollados por el autor, los cuales implementan las ideas formuladas en las obras mencionadas y aquí descritas. En el Capítulo 2, se estudia el FVM utilizando un mallado estructurado en geometrías variables para abordar problemas de convección-difusión térmica y el cálculo de velocidades en problemas con fenómenos de convección térmica. Se examina la aplicación de las ecuaciones transformadas a un nuevo espacio de coordenadas, destacando la aproximación de las derivadas de la incógnita y los esquemas de interpolación en las caras de las celdas. El Capítulo 3 amplía la formulación del FVM a celdas de forma arbitraria, lo que implica un estudio detallado de su conectividad. Además, en este capítulo se resuelve la ecuación de continuidad para determinar el campo de velocidades, conocimiento esencial para abordar problemas de convección. Finalmente, en el Capítulo 4, se lleva a cabo un análisis exhaustivo del trabajo desarrollado, sugiriendo posibles vías para la extensión de este proyecto. Los códigos desarrollados se presentan en su totalidad en el Apéndice B. Además, ha sido crucial validar el modelo; los ejemplos expuestos en cada sección se han comparado con las soluciones proporcionadas por el programa comercial Fluent®. El procedimiento seguido para obtener dichas soluciones se detalla en el Apéndice A.
2 Mallado Estructurado en coordenadas curvilíneas C omo se ha mencionado previamente, la solución numérica de la ecuación de transporte a través del Método de los Volúmenes Finitos (FVM) requiere la discretización de la geometría en volúmenes de control cerrados, cada uno con una normal exterior. En el centro de estos volúmenes de control (nodos), se almacena la información de la incógnita φ . Además, es necesario incorporar las condiciones de contorno. Estas condiciones pueden consistir en la especificación de la temperatura en la frontera o en la imposición del flujo de calor a través de la misma. Es conveniente utilizar el FVM para resolver problemas cuya solución analítica es inabordable o extremadamente compleja. Abordar estos problemas implica adaptarse a la morfología específica del sistema en cuestión, que puede ser no uniforme, a diferencia del cuadrado o rectángulo regular estudiado en [ 4 ]. Para geometrías que no son rectangulares, la creación de una malla estructurada cartesiana se vuelve altamente complicada, ya que las celdas dejan de ser rectángulos y la aproximación de derivadas en ambas direcciones resulta imprecisa y sumamente compleja. En la década de los 80, los investigadores llegaron a la conclusión de que la única forma precisa de llevar a cabo esta discretización era transformar las coordenadas cartesianas en coordenadas curvilíneas que se adaptaran a la geometría del problema [ 8 ]. En este nuevo espacio, la representación puede simplificarse como un cuadrado, donde las fronteras de este cuadrado coinciden exactamente con las fronteras en el dominio cartesiano. 2.1 FVM en coordenadas curvilíneas Resolver las ecuaciones para un mallado uniforme, donde todas las celdas tienen las mismas dimensiones y geometría, simplifica el problema. Sin embargo, al estudiar una morfología nueva que incluye celdas no uniformes, el cálculo de áreas, flujos normales y derivadas aproximadas de la incógnita se complica considerablemente. Por esta razón, se recurre a un cambio de variables apropiado mediante el cual la geometría se transforma en una más simple, típicamente un cuadrado (en 2D). En este nuevo espacio, todas las celdas tienen el mismo tamaño y las ecuaciones diferenciales se expresan en términos de las coordenadas transformadas. Las coordenadas cartesianas (x,y) representan las variables físicas habituales, conocido como el plano físico, mientras que su transformación a las nuevas variables (ξ,η) se denomina el plano computacional. En la figura Figura 2.1 se muestra un ejemplo en el que se ilustra lo anteriormente comentado. Es importante destacar que en sistemas no uniformes, las celdas tienden a agruparse en áreas de mayor complejidad geométrica, donde ocurren variaciones significativas en las propiedades, o en zonas de particular interés, entre otros factores. En este contexto, se considerará el caso bidimensional, donde la transformación directa de (x,y) a (ξ,η) en la mayoría de los casos se puede expresar como una función explícita. Sin embargo, en la transformación inversa, no es común contar con expresiones explícitas. El desarrollo que se presenta a continuación se basa en el enfoque propuesto por Mazumder [ 8 ], el cual se describe con mayor detalle en las secciones 2.8 y 6.6 de su trabajo. 11
12 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 (a) Plano físico. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b) Plano computacional. Figura 2.1 Mallado de una geometría. Esta transformación permite expresar las variables (x,y)como función de (ξ,η): x=x(ξ,η) y=y(ξ,η)(2.1) Como es sabido, la ecuación de transporte es una ecuación diferencial en derivadas parciales. En ella, es necesario expresar mediante la regla de la cadena las derivadas de la transformación inversa, ∂ξi ∂xj , las cuales no suelen estar disponibles mediante una expresión matemática: ∂φ ∂x=∂φ ∂ξ ·∂ξ ∂x+∂φ ∂η ·∂η ∂x ∂φ ∂y=∂φ ∂ξ ·∂ξ ∂y+∂φ ∂η ·∂η ∂y (2.2) De este modo, se presentan una serie de definiciones matemáticas relacionadas con la transformación de variables. Estas definiciones serán de gran utilidad para llevar la ecuación de transporte al espacio computacional. En particular, se introduce la matriz jacobiana [J] , su determinante J , la matriz de cofactores [β], y la matriz [B], que resulta del producto tensorial de [β]consigo misma. [J] = ∂x ∂ξ ∂x ∂η ∂y ∂ξ ∂y ∂η ,J=∂x ∂ξ ·∂y ∂η −∂x ∂η ·∂y ∂ξ (2.3) [β] = ∂y ∂η −∂y ∂ξ −∂x ∂η ∂x ∂ξ ,Bjk =βikβi j =β11β11 +β21β21 β12β11 +β22β21 β11β12 +β21β22 β12β12 +β22β22 (2.4)
2.1 FVM en coordenadas curvilíneas 13 Para mayor claridad, se ha optado por representar estas definiciones en formato matricial. Sin embargo, en lo sucesivo, se utilizará la notación tensorial debido a su claridad en la expresión de los resultados matemáticos, de los cuales es importante destacar la derivada ∂φ ∂xi =1 J ∂ ∂ξjβi jφ=1 J ∂φ ∂ξj βi j,(2.5) donde, debido al Teorema de Schwarz, la suma ∂βi j/∂ξj es igual a cero. Recuperando la ecuación de transporte de φ, se pueden aplicar los resultados anteriores para llevar a cabo la transformación: ∇·(ρu u uφ−Γ∇φ) = ∂ ∂xiρuiφ−Γ∂φ ∂xi=Sφ(2.6) Siguiendo la metodología empleada por Mazumder, puede reescribirse esta expresión de la siguiente manera: 1 J ∂ ∂ξjρUjφ−Γ JBk j ∂φ ∂ξk=Sφ(2.7) Donde, en esta ecuación, Uj=uiβi j representa la componente de velocidad normal a la superficie ξj constante. Dado que se está haciendo uso del Método de los Volúmenes Finitos, la ecuación se integra en un volumen de control rectangular, de lados ∆ξ y ∆η , resultando en una expresión donde las variables están evaluadas en las caras de la celda: ρU1φ−Γ JB11 ∂φ ∂ξ1 +B21 ∂φ ∂ξ2e ∆ξ2−ρU1φ−Γ JB11 ∂φ ∂ξ1 +B21 ∂φ ∂ξ2w ∆ξ2+ ρU2φ−Γ JB12 ∂φ ∂ξ1 +B22 ∂φ ∂ξ2n ∆ξ1−ρU2φ−Γ JB12 ∂φ ∂ξ1 +B22 ∂φ ∂ξ2s ∆ξ1= =SOJO∆ξ1∆ξ2 (2.8) donde, en esta ecuación, el subíndice en minúscula hace referencia a las caras de la celda y en subíndice O al nodo de la celda. En las secciones subsiguientes, se procederá a discretizar esta ecuación para abordar individualmente los casos de difusión pura y convección-difusión. En dichas secciones, se detallará cómo aproximar las derivadas y la variable en las caras, ya que, como se mostrará, este proceso no es trivial, especialmente para el caso de convección-difusión. Antes de proceder con la discretización, es esencial definir la geometría que será objeto de estudio en este trabajo. Esta geometría quedará delimitada por cuatro caras que serán nombradas en función de los puntos cardinales. Tres de estas caras son líneas rectas: la cara oeste, que tiene una altura de H1 ; la cara sur, que posee una longitud de L ; y la cara este, que tiene una altura de H2 . Además, está la cara norte, que conecta las caras oeste y este, y se define mediante una función arbitraria ys(x) . Un ejemplo de lo comentado se muestra en la Figura 2.1-(a). Para discretizar el dominio, simplemente se toman M+1 vértices equiespaciados en la dirección ξ y N+1 vértices equiespaciados en la dirección η , ambos tomando valores comprendidos entre 0 y 1, donde M y N son el número de celdas en la dirección ξ y η respectivamente. El cambio de variables que se plantea está basado en el comentado en [9]:
14 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas x=ξL(2.9a) y=ηys(x(ξ)),(2.9b) cuya matriz jacobiana correspondiente a este cambio de variables se expresa como: [J] = L0 ηys′(ξL)L ys(ξL) (2.10) Con esta formulación establecida, en la siguiente sección se explorará el problema estacionario de difusión pura en un contexto bidimensional. 2.2 Problema de difusión pura, estacionario y bidimensional En el caso estacionario de difusión térmica con generación de calor, al considerar la expresión integrada de la ecuación de transporte (Ecuación 2.8) y al omitir los términos convectivos, se obtiene: −Γ JB11 ∂φ ∂ξ +B21 ∂φ ∂η e ∆η+Γ JB11 ∂φ ∂ξ +B21 ∂φ ∂η w ∆η− −Γ JB12 ∂φ ∂ξ +B22 ∂φ ∂η n ∆ξ+Γ JB12 ∂φ ∂ξ +B22 ∂φ ∂η s ∆ξ= =SOJO∆ξ∆η (2.11) El proceso que procede, de manera similar a lo que Alejandro Cumplido realiza en [ 4 ], consiste en buscar una expresión matricial que permita representar la ecuación 2.11 como un sistema de ecuaciones lineales Aφ=bpara todas las celdas. De esta manera, podrá resolverse utilizando el comando barra invertida ‘\’ de MATLAB. En esencia, A es la matriz de tamaño MN ×MN que contiene los coeficientes que multiplican a cada magnitud φk , φ es el vector incógnita de tamaño MN ×1 , y b es el vector del mismo tamaño que φ que contiene la generación de calor y las condiciones de contorno en las celdas de la frontera. Por tanto, se discretizará la ecuación comentada aproximando las derivadas de φ por su valor en los nodos y se incluirán las condiciones de contorno. Luego, se organizarán estos coeficientes de manera que puedan ser almacenados en la matriz A. 2.2.1 Mallado del dominio El dominio en cuestión tiene la configuración de un cuadrado, con lados de una unidad de longitud, y se encuentra dividido en M×N celdas. Para facilitar la referencia a estas celdas, se ha adoptado la convención de utilizar índices, tal como se describe en la fuente de referencia [4]. Por un lado se tienen los índices (i,j) , donde i=1,2,...,M identifica las columnas de izquierda a derecha, y j=1,2,...,N , las filas desde abajo hacia arriba. Por otro lado, dadas las dimensiones del sistema matricial, será necesario utilizar un índice k para numerar cada una de las celdas, desde la esquina inferior izquierda hasta la esquina superior derecha, ya sea recorriendo las columnas o las filas. En este trabajo, se ha seleccionado el recorrido de las celdas por filas sin que exista una razón específica que lo justifique. La relación que existe entre estos índices se expresa como:
2.2 Problema de difusión pura, estacionario y bidimensional 15 k= (j−1)M+i(2.12) Durante el proceso de discretización, es esencial nombrar adecuadamente a todos los nodos vecinos involucrados en la ecuación. Como se ilustra en la Figura 2.2, estos nodos se identifican con los puntos cardinales en mayúscula, coincidiendo con aquellos que rodean a la celda k . Del mismo modo, las caras de la celda se etiquetan tabmién según los puntos cardinales, pero utilizando subíndices en minúscula. ξ η OE W N SSW SE NW NE ∆ξ ∆η 1 1 (a) O s e n w (b) Figura 2.2 Plano computacional. Por último, cabe mencionar que los incrementos ∆ξ y ∆η son constantes en todo el dominio, con valores de 1/M y 1/N respectivamente, ya que la información sobre la concentración de celdas en una localización específica se almacena en las funciones x=x(ξ,η)ey=y(ξ,η). 2.2.2 Discretización de las ecuaciones El siguiente paso consiste en desarrollar la ecuación de transporte 2.11 y particularizarla para cada volumen de control. Para aproximar las derivadas con respecto a ξ en las caras este y oeste, se pueden obtener empleando dos expresiones del polinomio de Taylor tal como indica [ 8 ] en la sección 6.6, y del mismo modo para las derivadas en la dirección η en las caras norte y sur. Esto se conoce como el método de las diferencias centradas, que para mallas uniformes, estas aproximaciones son de segundo orden. ∂φ ∂ξ e ≈φE−φO ∆ξ,∂φ ∂ξ w ≈φO−φW ∆ξ, ∂φ ∂η n ≈φN−φO ∆η,∂φ ∂η s ≈φO−φS ∆η (2.13) Las derivadas de φ con respecto a η en las caras este y oeste, y con respecto a ξ en las caras norte y sur, expresiones que en problemas de geometrías cuadrangulares no aparecían, representan la no ortogonalidad de las líneas ξ y η . Para aproximar su valor, se estima primero el valor de los vértices de la cara según los valores de φ en las celdas que contienen el vértice, y posteriormente se aplica el método de las diferencias centradas como se ha hecho para las anteriores derivadas: ∂φ ∂η e ≈ φN+φNE +φO+φE 4−φS+φSE +φO+φE 4 ∆η=φN+φNE −φS−φSE 4∆η(2.14)
16 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas Del mismo modo para las caras restantes: ∂φ ∂η w ≈φN+φNW −φS−φSW 4∆η, ∂φ ∂ξ s ≈φE+φSE −φW−φSW 4∆ξ, ∂φ ∂ξ n ≈φE+φNE −φW−φNW 4∆ξ. (2.15) Sin embargo, estas aproximaciones son únicamente válidas en las celdas interiores, ya que en las caras de la frontera no se dispondrán de los valores más allá de la celda Ohacia afuera. En este caso, la expresión resultante dependerá de si la condición de contorno impuesta es de temperatura o de flujo impuesto, que son las dos que se tratarán en este trabajo. Condición de Dirichlet: temperatura impuesta Se trata del caso en el que la incógnita φ es conocida en la frontera. Sin embargo, al aplicar el Método de los Volúmenes Finitos, la implementación de esta condición no es directa puesto que el centro de las celdas no coincide con la frontera. Por lo tanto, es necesario obtener una expresión que aproxime la derivada en esta zona con precisión de segundo orden. Para ello, se recurre nuevamente a dos expresiones del polinomio de Taylor. Se ejemplificará para el caso de la cara este, pero el proceso es completamente análogo para el resto de caras: I: φO=φECC +∂φ ∂ξ e−∆ξ 2+∂2φ ∂ξ2e−∆ξ 221 2+O(∆ξ3) II: φW=φECC +∂φ ∂ξ e−3∆ξ 2+∂2φ ∂ξ2e−3∆ξ 221 2+O(∆ξ3) (2.16) Operando con estas dos expresiones, se cancelan los términos de segundo orden haciendo 9(I)−(II) , de lo que se deduce: ∂φ ∂ξ e =8φECC +φW−9φO 3∆ξ(2.17) Del mismo modo, para el resto de caras se obtienen expresiones similares: ∂φ ∂η n =8φNCC +φS−9φO 3∆η ∂φ ∂ξ w =9φO−φE−8φWCC 3∆ξ ∂φ ∂η s =9φO−φN−8φSCC 3∆η (2.18) Condición de Neumann: flujo impuesto Este tipo de condición de contorno impone el valor de la derivada de la incógnita en la dirección normal exterior a la frontera:
2.2 Problema de difusión pura, estacionario y bidimensional 17 qCCf=n n nf·(−Γ∇φ)(2.19) n n nf f f q=qFCC O Figura 2.3 Flujo y vector normal en la frontera. Para desarrollar esta expresión, es necesario traducir la ecuación al ámbito computacional. Para ello, se acude a la expresión 2.5: n n nf·(−Γ∇φ) = ni −Γβi j J ∂φ ∂ξj!=−nx Γβ11 Jf ∂φ ∂ξ −Γβ12 Jf ∂φ ∂η !−ny Γβ21 Jf ∂φ ∂ξ +Γβ22 Jf ∂φ ∂η !(2.20) Esto quiere decir que la expresión 2.20 equivale a un valor conocido y puede ser sustituido en la ecuación 2.11. Sin embargo, esta suma de derivadas no aparece directamente en la ecuación de balance energético en coordenadas computacionales 2.11, pero sí puede encontrarse al desarrollar los términos Bi j reagrupando los términos: Γ JB12 ∂φ ∂ξ +B22 ∂φ ∂η n ∆ξ=Γ J(β12β11 +β22β21)∂φ ∂ξ +(β12β12 +β22β22)∂φ ∂η n ∆ξ= =Γ Jnβ21 ∂φ ∂ξ n +β22 ∂φ ∂η nβ22 +Γ Jnβ11 ∂φ ∂ξ n +β12 ∂φ ∂η nβ12∆ξ= =Γ Jnβ11 ∂φ ∂ξ n +β12 ∂φ ∂η n | {z } −qxn (β12 −nx) + Γ Jnβ21 ∂φ ∂ξ n +β22 ∂φ ∂η n | {z } −qyn β22 −ny +Γ Jn nxβ11 ∂φ ∂ξ n +β12 ∂φ ∂η n+Γ Jn nyβ21 ∂φ ∂ξ n +β22 ∂φ ∂η n | {z } −qNCC ∆ξ (2.21) Se ha elegido denominar a las expresiones identificadas como qxn y qyn debido a que, teniendo en cuenta la expresión 2.5, representan la componente xy la componente ydel flujo −Γ∇φen la cara norte. 2.2.3 Aplicación práctica: código de MATLAB Hasta este punto, se ha logrado expresar la ecuación 2.11 en función de la incógnita en los nodos de todas las celdas. Por lo tanto, es ahora factible implementar la solución del sistema de ecuaciones que plantea el
24 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas estimación más precisa. A partir de aquí, comienza un proceso iterativo donde la presión, velocidad y densidad se actualizan hasta que las variaciones son suficientemente pequeñas, alcanzando la convergencia. Para comenzar este proceso son necesarias estimaciones iniciales de los campos de velocidad y presión. Con el paso del tiempo, han surgido variaciones de este método que, a pesar de aumentar el número de cálculos, reducen el coste computacional y el tiempo de ejecución. Un ejemplo es el algoritmo SIMPLER (SIMPLE Revised), desarrollado por el mismo Patankar ([ 11 ]) en 1980, donde se deduce una nueva ecuación discretizada para la corrección de la presión obteniendo mejores aproximaciones por cada iteración. Un método alternativo se obtiene al omitir algunos términos menos relevantes en las ecuaciones de la velocidad de SIMPLE, nombrado como SIMPLEC (SIMPLE Consistent), desarrollado por Van Doormal y Raithby en 1984. Por último, dos años más tarde surge un nuevo método que podría entenderse como una extensión de SIMPLE: PISO (Pressure Implicit with Splitting Operators), que se basa en reforzar el proceso correctivo de las variables para mejorar cada iteración. Estos algoritmos resuelven el campo de velocidades para un sistema de ecuaciones cuando existe un acoplamiento entre la presión y velocidad. Sin embargo, una amplia parte de los problemas fluidodinámicos pueden ser simplificados bajo una serie de justificadas hipótesis. Tal es el caso de los flujos incompresibles. En efecto, como se verá a continuación, bajo dichas condiciones, el campo de velocidades es solenoidal, es decir, su divergencia es nula, por lo que en la ecuación de la energía 1.3 la presión desaparece, permitiendo calcular las velocidades de manera independiente. 2.3.1 Función de corriente para flujo estacionario e incompresible El movimiento estacionario de flujos incompresibles ( ρ= cte) puede ser descrito directamente mediante la ecuación de continuidad ∂u ∂x+∂v ∂y=0,(2.22) en la que dada su particularidad, el teorema de la igualdad de las derivadas cruzadas permite definir una función de corriente φde la que deriven las componentes de velocidad: u=∂φ ∂yv=−∂φ ∂x(2.23) Esta definición permite reducir el número de incógnitas de 2 para el caso bidimensional, y 3 para el tridimensional, a una sola para ambos casos. Además, tal y como se demuestra en [ 12 ], la función de corriente se mantiene constante a lo largo de las líneas de corriente y el flujo volumétrico o caudal a través de una línea que une dos puntos A y B en el caso bidimensional es igual a φA−φB . Esta propiedad favorece el establecimiento de la condición de contorno a través de un conducto como el de la figura, donde en la pared inferior se establece que φB=0y en la pared superior φA=q. A B n Figura 2.10 Caudal a través de una línea. Esto puede adaptarse sin mayor complicación a la geometría que fue definida al principio de esta sección. A través de la cara oeste entra el fluido a una determinada temperatura y velocidad, y por la cara este saldría el mismo caudal con condición de temperatura o flujo impuesto. Por otro lado, las caras norte y sur actúan como
2.3 Problema de convección-difusión térmica 25 las paredes del conducto. Para el caso concreto de flujos incompresibles, como se ha visto, la temperatura no interviene en este problema y pueden ser obviadas dichas condiciones. Acorde a la anterior descripción, se ilustra el problema del cálculo de velocidades: L H Figura 2.11 Conducto bidimensional. Dadas las cualidades del problema, es posible servirse del código de difusión para resolver este problema, ya que las dos únicas diferencias que los separan son que Γ=1 , y que en toda la frontera el valor de la función de corriente φes impuesto. La ecuación a resolver en el plano computacional resultaría: −B11 ∂φ ∂ξ +B21 ∂φ ∂η e ∆η+B11 ∂φ ∂ξ +B21 ∂φ ∂η w ∆η− −B12 ∂φ ∂ξ +B22 ∂φ ∂η n ∆ξ+B12 ∂φ ∂ξ +B22 ∂φ ∂η s ∆ξ=0 (2.24) En cuanto a las condiciones de contorno en las caras oeste y este, por donde entra y sale el fluido respectivamente, el valor de la función de corriente en un punto ( x0,y0 ) puede calcularse a través de la fórmula del flujo que atraviesa una línea desde la cara sur hasta el punto ( x0,y0 ). Para este particular caso, la línea será el segmento perteneciente a la cara oeste que comienza en la esquina suroeste y es de altura y. Q(y) = Zy 0 v·nds =Zy 0 u(y)dy =Zy 0 ∂φ ∂ydy =Zφ(y) φ0 dφ=φ(y)−0(2.25) Q(y) = Zy 0 v·nds =Zy 0 u(s)ds =uy (2.26) Donde se han supuesto, por simplicidad, velocidades horizontales u constantes en la dirección y en la entrada y la salida, sin que esta condición sea estrictamente necesaria. Por tanto, igualando las dos expresiones anteriores, se obtiene el valor de la función de corriente para cualquier punto de la frontera: φW(y) = uey(2.27) φE(y) = usy(2.28) φN=ueHe(2.29) φS=0(2.30)
26 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas Dado que la formulación se realiza en términos de las variables (ξ,η) , es preciso expresar este resultado en dicho dominio: φW(η) = ueηHe(2.31) φE(η) = usηHs(2.32) φN=ueHe(2.33) φS=0(2.34) Es importante tener en cuenta que solo una de las dos velocidades es proporcionada por los datos del problema, ya sea ue o us . Sin embargo, como no existe pérdida de masa, el flujo a la entrada y a la salida debe ser el mismo, lo que implica que: us=ueHe Hs (2.35) Al igual que en el problema de difusión, la incógnita φ será calculada en los nodos de las celdas. Sin embargo, dado que los términos convectivos de la ecuación 2.8 son evaluados en las caras, las velocidades también se calcularán en tales lugares. Así, una vez obtenida la incógnita, las velocidades pueden deducirse gracias a la expresión 2.5: u(x,y) = 1 J∂φ ∂ξ β21 +∂φ ∂η β22(2.36) v(x,y) = −1 J∂φ ∂ξ β11 +∂φ ∂η β12(2.37) El código será esencialmente similar al utilizado para la difusión, por lo que se puede partir de ese mismo. También se beneficiará del script de configuración, que se nombrará como “ config_conv.m ” para diferenciarlo del anterior y se empleará en el programa principal para calcular la distribución de temperaturas. La única modificación es la adición de dos parámetros: la velocidad de entrada ue y la capacidad calorífica a volumen constante del fluido Cv, factor presente en el término convectivo. El programa que se desarrolla para calcular la solución puede almacenarse como función, devolviendo como valores de salida las velocidades en las caras. Esta decisión se motiva por la eficiencia en el uso del espacio en memoria, ya que tras finalizar la ejecución de la función, todas sus variables locales son eliminadas. Por otro lado, las secciones que se diferencian del programa de difusión son tres. La primera es la sección de las condiciones de contorno de la función de corriente. Estas condiciones son constantes en las caras norte y sur, pero en las caras este y oeste son funciones lineales que van desde φS a φN . Además, es necesario proporcionar sus derivadas. 1% Datos del flujo 2 3H1 = H; 4H2 = ys(L); 5v2 = v1*H1/H2; 6 7deltaxi = 1/M; 8deltaeta = 1/N; 9A=sparse(n_c,n_c); 10 b=sparse(1,n_c); 11 phi_scc = 0; 12 phi_ncc = v1*H1; 13 phi_ecc = @(eta) v2*H2*eta; 14 dphi_e_deta = @(eta) v2*H2;
2.3 Problema de convección-difusión térmica 27 15 phi_wcc = @(eta) v1*H1*eta; 16 dphi_w_deta = @(eta) v1*H1; En segundo lugar, dentro del bucle, es preciso añadir a las expresiones de las derivadas de φ en las caras este y oeste de la celda, la expresión de la derivada de ∂φ/∂η en las caras este y oeste del dominio. Al resto de derivadas se le eliminarían todas las expresiones relacionadas con flujo de calor impuesto y la variable binaria bphiimp_f puede ser suprimida. 1dphideta_e = ... + dphi_e_deta(etae(k))*istepe(i,j)*b_v; 2dphideta_w = ... + dphi_w_deta(etaw(k))*istepw(i,j)*b_v; La tercera y última modificación del código consiste en calcular las velocidades en las caras de las celdas a partir de la solución. Nuevamente, se recorrerán las celdas una vez φ sea conocida. En este bucle, el algoritmo desarollado en la Sección 2.2.3 para discriminar los índices que representaban celdas más allá del dominio es reutilizado para crear nueve variables que tomen el valor de φ en la celda y en sus celdas vecinas ( φO , φE , φNE , ...). Por ende, pueden calcularse las derivadas para hacer el cálculo de las velocidades según las expresiones 2.37. Por último, pueden calcularse las velocidades en el nodo de la celda como la media de velocidades de sus cuatro caras. 1vx_f = zeros(M*N,4); 2vy_f = zeros(M*N,4); 3for j=1:N 4for i=1:M 5k = (j-1)*M+i; 6index = [k-M-1 k-M k-M+1 k-1 k k+1 k+M-1 k+M k+M+1]; 7starts_index = index(index > 0); 8end_index = index(index <= M*N & index > 0); 9eq_start = find(index == starts_index(1)); 10 eq_end = find(index == end_index(end)); 11 a=zeros(size(index)); 12 a(1,eq_start:eq_end) = phi1D(index(eq_start:eq_end))’; 13 14 phi_SW = a(1,1); phi_S = a(1,2); phi_SE = a(1,3); 15 phi_W = a(1,4); phi_O = a(1,5); phi_E = a(1,6); 16 phi_NW = a(1,7); phi_N = a(1,8); phi_NE = a(1,9); 17 18 dphidxi_e = (phi_E - phi_O)/deltaxi*stepe(i,j) + (phi_W + 8*phi_ecc(etae(k)) - 9*phi_O)/(3*deltaxi)*istepe(i,j); 19 dphideta_e = (phi_NE + phi_N - phi_SE - phi_S)/(4*deltaeta)*stepe(i,j)*stepn(i,j )*steps(i,j) ... 20 + (phi_ncc-(phi_S + phi_SE + phi_O + phi_E)/4)/deltaeta*istepn(i,j)*stepe(i,j ) ... 21 + ((phi_N + phi_NE + phi_O + phi_E)/4 - phi_scc)/(deltaeta)*isteps(i,j)*stepe (i,j) ... 22 + dphi_e_deta(etae(k))*istepe(i,j); 23 dphidxi_w = (phi_O - phi_W)/deltaxi*stepw(i,j) + (9*phi_O - phi_E - 8*phi_wcc( etaw(k)))/(3*deltaxi)*istepw(i,j); 24 dphideta_w = (phi_NW + phi_N - phi_SW - phi_S)*steps(i,j)/(4*deltaeta)*stepw(i,j )*stepn(i,j)*steps(i,j) ... 25 + (phi_ncc-(phi_S + phi_SW + phi_O + phi_W)/4)/deltaeta*istepn(i,j)*stepw(i,j ) ... 26 + ((phi_N + phi_NW + phi_O + phi_W)/4 - phi_scc)/(deltaeta)*isteps(i,j)*stepw (i,j) ... 27 + dphi_w_deta(etaw(k))*istepw(i,j); 28 dphidxi_s = (phi_E + phi_SE - phi_SW - phi_W)/(4*deltaeta)*steps(i,j)*stepe(i,j) *stepw(i,j) ...
28 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas 29 + (phi_ecc(etas(k))-(phi_SW + phi_W + phi_O + phi_S)/4)/(deltaeta)*istepe(i,j )*steps(i,j) ... 30 + ((phi_SE + phi_E + phi_O + phi_S)/4 - phi_wcc(etas(k)))/(deltaeta)*istepw(i ,j)*steps(i,j); 31 dphideta_s = (phi_O - phi_S)/deltaxi*steps(i,j) + (9*phi_O - phi_N - 8*phi_scc) /(3*deltaeta)*isteps(i,j); 32 dphidxi_n = (phi_NE + phi_E - phi_NW - phi_W)/(4*deltaeta)*stepn(i,j)*stepe(i,j) *stepw(i,j) ... 33 + (phi_ecc(etan(k))-(phi_NW + phi_W + phi_O + phi_N)/4)/(deltaeta)*istepe(i,j )*stepn(i,j) ... 34 + ((phi_NE + phi_E + phi_O + phi_N)/4 - phi_wcc(etan(k)))/(deltaeta)*istepw(i ,j)*stepn(i,j); 35 dphideta_n = (phi_N - phi_O)/deltaxi*stepn(i,j) + (phi_S + 8*phi_ncc - 9*phi_O) /(3*deltaeta)*istepn(i,j); 36 37 vy_f(k,1) = - 1/detJ(xis(k),etas(k))*(beta11(xis(k),etas(k))*dphidxi_s + beta12( xis(k),etas(k))*dphideta_s); 38 vy_f(k,2) = - 1/detJ(xie(k),etae(k))*(beta11(xie(k),etae(k))*dphidxi_e + beta12( xie(k),etae(k))*dphideta_e); 39 vy_f(k,3) = - 1/detJ(xin(k),etan(k))*(beta11(xin(k),etan(k))*dphidxi_n + beta12( xin(k),etan(k))*dphideta_n); 40 vy_f(k,4) = - 1/detJ(xiw(k),etaw(k))*(beta11(xiw(k),etaw(k))*dphidxi_w + beta12( xiw(k),etaw(k))*dphideta_w); 41 vx_f(k,1) = 1/detJ(xis(k),etas(k))*(beta21(xis(k),etas(k))*dphidxi_s + beta22( xis(k),etas(k))*dphideta_s); 42 vx_f(k,2) = 1/detJ(xie(k),etae(k))*(beta21(xie(k),etae(k))*dphidxi_e + beta22( xie(k),etae(k))*dphideta_e); 43 vx_f(k,3) = 1/detJ(xin(k),etan(k))*(beta21(xin(k),etan(k))*dphidxi_n + beta22( xin(k),etan(k))*dphideta_n); 44 vx_f(k,4) = 1/detJ(xiw(k),etaw(k))*(beta21(xiw(k),etaw(k))*dphidxi_w + beta22( xiw(k),etaw(k))*dphideta_w); 45 end 46 end 47 vx_c=sum(vx_f(:,[1:4])’)/4; 48 vy_c=sum(vy_f(:,[1:4])’)/4; De los datos obtenidos de esta función, se pueden representar la función de corriente, las líneas de corriente, y las velocidades resultantes. Figura 2.12 Campo de velocidades a través de un conducto exponencial.
2.3 Problema de convección-difusión térmica 29 Figura 2.13 Líneas de corriente a través de un conducto exponencial. Figura 2.14 Función de corriente en un conducto exponencial. 2.3.2 Aproximación del término convectivo Una vez que se ha obtenido la velocidad del flujo en cada una de las caras de la celda, es posible incluir dichos valores en la ecuación 2.8. Sin embargo, la aproximación de la incógnita en las caras de la celda requiere especial atención. En efecto, como se analizó anteriormente por Alejandro Cumplido en [ 4 ] en la Sección 2.2.2 de su trabajo, para garantizar la convergencia en problemas convectivo-difusivos, es necesario que el esquema de discretización empleado cumpla con tres propiedades: conservatividad, delimitación y transportividad. La conservatividad establece que el flujo de φ entrante a través de una cara del volumen de control debe ser idéntico al flujo entrante en el volumen de control adyacente a través de la cara fronteriza entre ambos. La delimitación implica el cumplimiento, por un lado, del criterio de Scarborough, que especifica que el elemento de la diagonal principal debe ser mayor que el resto de su fila, y, por otro lado, que todos los coeficientes del sistema tengan el mismo signo. Finalmente, la transportividad contempla la direccionalidad del flujo y la relación existente entre convección y difusión. Esta última condición se sirve del número de Peclet (Pe):
30 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas Pe =ρu Γ/∆x(2.38) donde ∆xes la longitud característica, en este caso, el ancho de la celda. En el trabajo mencionado, se demostró que el método de las diferencias centradas no cumple con la condición de delimitación, ya que para números de Peclet mayores que 2, no se satisface la condición de Scarborough. Además, este método no contempla la importancia relativa del flujo convectivo frente al difusivo. Por esta razón, este método puede no proporcionar resultados fiables y se opta por aplicar directamente el esquema upwind de primer orden, o FOU (First Order Upwind). En dicho método, el valor de la incógnita en una cara de la celda se establece en base al valor del nodo que se sitúa aguas arriba. Matemáticamente se expresa de la siguiente manera: φw=φW,φe=φPsi uw>0,ue>0 φw=φO,φe=φEsi uw<0,ue<0 φs=φS,φn=φOsi vs>0,vn>0 φs=φO,φn=φOsi vs<0,vn<0 Esto puede llevarse a la práctica mediante una función de MATLAB de la siguiente manera: 1function [phi_r] = upwind(v, phi_1, phi_2) 2phi_r = phi_1*(sign(v) + 1)/2 + phi_2*(1 - sign(v))/2; 3end Donde phi_1 y phi_2 son las incógnitas de los nodos de las celdas adyacentes que avanzan en el sentido positivo del eje ξ o η . Por ejemplo, para aproximar la cara oeste, el argumento phi_1 sería φW y phi_2 sería φO . Por otro lado el argumento v es la componente de velocidad en x para las caras este y oeste, e y para las caras norte y sur. Aplicar las condiciones de contorno para el esquema de upwind de primer orden no presenta mayor complicación. Si la cara tratada se encuentra en la frontera, la incógnita φ se sustituirá por su valor de condición de contorno cuando la temperatura sea impuesta. En caso contrario, se utilizará φO de la propia celda si se impone un flujo de calor. Para completar la formulación del término convectivo, es necesario calcular las velocidades contravariantes 2.37 y agregarlas a la ecuación: 1for j=1:N 2for i=1:M 3U1e = vx_f(k,2)*beta11(xie(k),etae(k)) + vy_f(k,2)*beta21(xie(k),etae(k)); 4U1w = vx_f(k,4)*beta11(xiw(k),etaw(k)) + vy_f(k,4)*beta21(xiw(k),etaw(k)); 5U2n = vx_f(k,3)*beta12(xin(k),etan(k)) + vy_f(k,3)*beta22(xin(k),etan(k)); 6U2s = vx_f(k,1)*beta12(xis(k),etas(k)) + vy_f(k,1)*beta22(xis(k),etas(k)); 7... 8phi_re = upwind(vx_f(k,2), phi_O, phi_E)*stepe(i,j) + phi_O*istepe(i,j)*bfimp_e + phi_ecc*istepe(i,j)*b_v*bphiimp_e; 9phi_rw = upwind(vx_f(k,4), phi_W, phi_O)*stepw(i,j) + phi_O*istepw(i,j)*bfimp_w + phi_wcc*istepw(i,j)*b_v*bphiimp_w; 10 phi_rn = upwind(vy_f(k,3), phi_O, phi_N)*stepn(i,j) + phi_O*istepn(i,j)*bfimp_n + phi_ncc*istepn(i,j)*b_v*bphiimp_n; 11 phi_rs = upwind(vy_f(k,1), phi_S, phi_O)*steps(i,j) + phi_O*isteps(i,j)*bfimp_s + phi_scc*isteps(i,j)*b_v*bphiimp_s; 12 13 eq = (rho*U1e*Cv*phi_re + qEcc*beta11(xie(k),etae(k)) + qE2*beta21(xie(k),etae(k )))*deltaeta ...
2.3 Problema de convección-difusión térmica 31 14 - (rho*U1w*Cv*phi_rw + qWcc*beta11(xiw(k),etaw(k)) + qW2*beta21(xiw(k),etaw(k )))*deltaeta ... 15 + (rho*U2n*Cv*phi_rn + qNcc*beta22(xin(k),etan(k)) + qN2 + qncc*b_v*istepn(i, j)*bfimp_n)*deltaxi ... 16 - (rho*U2s*Cv*phi_rs + qScc*beta22(xis(k),etas(k)) + qS2*beta12(xis(k),etas(k )))*deltaxi; 17 end 18 end 2.3.3 Aplicación práctica Para poner a prueba estos códigos, considérese la geometría de la figura 2.11, con dimensiones H=2 m y L=2 m, y una función superior dada por y=1.9775e−x/L+0.0225 . La velocidad del flujo entrante por la cara entrante tomará valores en el intervalo [10−6,10−2] m/s. El fluido estudiado es aire, con las siguientes propiedades: Cv=718 J/(kgK) , ρ=1.225 kg/m3 , Γ=0.0242 ( W/(m2K) ). En cuanto a las condiciones de temperatura, la cara oeste estará sometida a 30◦ C, las caras norte y este estarán aisladas térmicamente, y la cara sur recibirá un flujo de calor entrante de 7 W/m2. (a) Solución de MATLAB para M=100 yN=100.(b) Solución de Fluent para el problema convectivo. Figura 2.15 Comparación de soluciones para v=10−6m/s. (a) Solución de MATLAB para M=100 yN=100.(b) Solución de Fluent para el problema convectivo. Figura 2.16 Comparación de soluciones para v=10−4m/s. Varios comentarios pueden hacerse al observar los resultados obtenidos por ambos programas. En primer
32 Capítulo 2. Mallado Estructurado en coordenadas curvilíneas (a) Solución de MATLAB para M=100 yN=100.(b) Solución de Fluent para el problema convectivo. Figura 2.17 Comparación de soluciones para v=10−2m/s. lugar, a medida que aumenta la velocidad, los términos convectivos cobran mayor importancia, incrementando el transporte de la magnitud φ desde la cara oeste hacia una mayor extensión del dominio. Por lo tanto, se observa que una parte cada vez mayor de la región adopta un tono azulado (temperatura cercana a 30◦C ). En segundo lugar, el gradiente de temperaturas es más dilatado en la zona cercana a la pared en la figura generada por Fluent. Esto puede deberse a que este último considera la capa límite de esta zona, donde las velocidades varían desde 0 hasta la velocidad prevista por la teoría potencial, reduciendo en consecuencia la importancia relativa de la convección frente a la difusión en dicha zona. Además, el esquema empleado para resolver el campo de velocidades es el SIMPLE. Por lo tanto, estas razones podrían ser, en parte, el origen de las disimilitudes entre los valores de temperatura. (a) Solución de MATLAB para M=100 yN=100.(b) Solución de Fluent del campo de velocidades. Figura 2.18 Comparación de velocidades para v=10−4m/s. Otro motivo por el que las soluciones pueden diferir es la ecuación de la energía resuelta por Fluent. Según su guía teórica ([13]), esta ecuación se formula de la siguiente manera: ∂ ∂t(ρE)+∇·(v(ρE+p)) = ∇· keff∇T−∑ j hjJj+(τeff ·v)!+Sh(2.39) donde keff es la conductividad efectiva ( k+kt , siendo kt la conductividad térmica turbulenta), que en este caso es igual a la conductividad térmica del material ya que se ha supuesto flujo laminar, y Jj es el flujo difusivo de las especies j . Dado que no se trata de una mezcla, el flujo es laminar sin viscosidad, no hay
2.3 Problema de convección-difusión térmica 33 generación de calor y el régimen es estacionario, la ecuación resultante debería ser: ∇·(v(ρE+p)) = ∇·(k∇T),(2.40) en la cual aparecen varios términos que no han sido considerados en el modelo desarrollado en este trabajo. Estos términos a un gradiente de presiones no nulo. En efecto, para un flujo estacionario incompresible no viscoso y en ausencia de fuerzas másicas, el gradiente de presiones es: ∇p=−ρv·∇v(2.41) A pesar de estas diferencias, la distribución de la temperatura obtenida con los códigos de MATLAB es casi indistinguible para la mayoría de los casos, y proporciona un punto de partida sólido para futuros desarrollos.
40 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado – Número total de vértices: nc – Coordenadas de todos los nodos: (xc,yc) – Coordenadas de todos los centros de las caras: (xf,yf) – Coordenadas de todos los vértices: (xv,yv) •Información a nivel local – Número total de caras de cada celda: nf(ic) – Número total de vértices de cada celda: nv(ic) – Número total de vértices de cada cara: nvf(if) •Conectividad – Conectividad Celda-A-Cara: Lc f (ic,k) – Conectividad Cara-A-Celda: L f c(if,k) – Conectividad Celda-A-Vértice: Lcv(ic,k) – Conectividad Cara-A-Vértice: L f v(ic,k) Para abordar este mallado se podría comenzar ubicando primero los vértices mediante el cambio de variable 2.9, que podrían seleccionarse equiespaciados (o no) en los ejes (ξ,η) , ambos tomando valores comprendidos entre 0 y 1. La siguiente figura ilustra la indexación de los elementos de la malla en términos globales. Figura 3.4 Indexación global de una malla no estructurada. En este caso bidimensional, donde todas las celdas son cuadrangulares, la información local es invariable: nf(∀ic) = 4 , nv(∀ic) = 4 y nvf(∀if) = 2 . Una vez determinadas las posiciones de todos los vértices, la ubicación y numeración de los demás elementos pueden derivarse de estos. Para identificar elementos específicos, como el baricentro de las celdas, es esencial conocer previamente ciertas características del
3.2 Problema difusivo, estacionario, bidimensional 41 mallado, como los vértices de cada celda. Este procedimiento determina el orden en que se genera y almacena la información necesaria. A continuación, se presentan los códigos de MATLAB que calculan los datos de la malla en su totalidad. Estas líneas de código están basadas en los desarrollados por Alejandro Cumplido en [4]. 1%% Mallado 2 3H=2; 4L=2; 5Nx=5; % Número de vértices en la dirección x 6Ny=5; % Número de vértices en la dirección y 7M=Nx-1; % Número de celdas en x 8N=Ny-1; % Número de celdas en y 9 10 x_ini=linspace(0,L,Nx); 11 eta=linspace(0,1,Ny); 12 13 %% Definición de la geometría 14 % Límite superior y = s2(x) 15 16 % Elipse 17 a = L*1.2; b = H; 18 S2 = @(x) b*sqrt(1-x.^2/a^2); 19 s2=S2(x_ini); 20 21 xvert = zeros(Nx,Ny); 22 yvert = zeros(Nx,Ny); 23 for i=1:Nx 24 xvert(i,1:Ny) = x_ini(i); 25 yvert(i,1:Ny) = eta*s2(i); 26 end 27 28 %% Información a nivel Global 29 % Número total de nodos 30 n_c = M*N; 31 32 % Número total de vértices 33 34 n_v = Nx*Ny; 35 36 % Coordenadas de los vértices 37 38 iv = 0; 39 xv=zeros(1,n_v); 40 yv=zeros(1,n_v); 41 for j=1:Ny 42 for i=1:Nx 43 iv = iv + 1; 44 xv(iv) = xvert(i,j); 45 yv(iv) = yvert(i,j); 46 end 47 end 48 49 % Coordenadas de los nodos 50 51 ic = 0; 52 xc=zeros(1,n_c);
42 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado 53 yc=zeros(1,n_c); 54 for j=1:N 55 for i=1:M 56 ic = ic + 1; 57 iv = ic + (j-1); % Vértice asociado a la esquina inferior izquierda de la celda 58 xc(ic) = (xv(iv) + xv(iv+1))/2; 59 yc(ic) = (yv(iv) + yv(iv+Nx) + yv(iv+1) + yv(iv+Nx+1))/4; 60 end 61 end Posteriormente, se procede al procesamiento de la conectividad entre las celdas y las caras, así como entre las caras y las celdas, los vértices y las celdas, y finalmente entre las caras y los vértices. Es importante destacar que el cálculo del centro de las caras no se presenta a continuación, ya que implica utilizar la conectividad entre las caras y los vértices, obteniéndose como la media aritmética de los vértices de cada cara. Los siguientes fragmentos de código implementan las conexiones mencionadas anteriormente. 1%% Conectividad 2% Conectividad Celda-A-Cara 3 4ic = 0; 5Lcf = zeros(n_c,4); 6for j=1:N 7for i=1:M 8ic = ic + 1; 9Lcf(ic,1) = ic; 10 Lcf(ic,2) = M*Ny + N*i + j; 11 Lcf(ic,3) = ic + M; 12 Lcf(ic,4) = M*Ny + N*(i-1) + j; 13 end 14 end 15 16 n_f = Lcf(ic,2); 17 18 % Conectividad Cara-A-Celda 19 Lfc = zeros(n_f,2); 20 for iface=1:n_f 21 c1 = 0; c2 = 0; 22 for ic=1:n_c 23 for k=1:4 24 if Lcf(ic,k) == iface 25 if c1 == 0 26 c1 = ic; 27 else 28 c2 = ic; 29 end 30 end 31 end 32 end 33 Lfc(iface,1) = c1; 34 Lfc(iface,2) = c2; 35 if Lfc(iface,1) == 0 36 Lfc(iface,1) = Lfc(iface,2); 37 elseif Lfc(iface,2) == 0 38 Lfc(iface,2) = Lfc(iface,1); 39 end 40 end
3.2 Problema difusivo, estacionario, bidimensional 43 41 42 43 % Conectividad Celda-A-Vértice 44 ic = 0; 45 Lcv = zeros(n_c,4); 46 for j=1:N 47 for i=1:M 48 ic = ic + 1; 49 Lcv(ic,1) = ic + (j-1); 50 Lcv(ic,2) = ic + j; % Es uno más que el primer vértice de la celda 51 Lcv(ic,3) = Nx*j + (i+1); 52 Lcv(ic,4) = Nx*j + i; 53 end 54 end 55 56 % Conectividad Cara-A-Vértice 57 Lfv = zeros(n_f,2); 58 for iface=1:M*N 59 ic = Lfc(iface,2); 60 Lfv(iface,1) = Lcv(ic,1); 61 Lfv(iface,2) = Lcv(ic,2); 62 end 63 for iface=(M*N+1):M*Ny 64 ic = Lfc(iface,2); 65 Lfv(iface,1) = Lcv(ic,4); 66 Lfv(iface,2) = Lcv(ic,3); 67 end 68 for iface=(M*Ny+1):(M*Ny+N*M) 69 ic = Lfc(iface,2); 70 Lfv(iface,1) = Lcv(ic,1); 71 Lfv(iface,2) = Lcv(ic,4); 72 end 73 for iface=(M*Ny+N*M+1):n_f 74 ic = Lfc(iface,2); 75 Lfv(iface,1) = Lcv(ic,2); 76 Lfv(iface,2) = Lcv(ic,3); 77 end 78 79 % Centros de las caras 80 xf=zeros(n_f,1); 81 yf=zeros(n_f,1); 82 for iface=1:n_f 83 xf(iface) = (xv(Lfv(iface,1))+xv(Lfv(iface,2)))/2; 84 yf(iface) = (yv(Lfv(iface,1))+yv(Lfv(iface,2)))/2; 85 end Con respecto a estas últimas líneas, cabe realizar una aclaración, existen caras que únicamente pertenecen a una celda, que es el caso de las fronterizas, por lo que en la matriz de conectividad L f c(ic,k) se le ha asignado la misma celda para k=1yk=2. En relación a la información del mallado, resta obtener los vectores normales nf, tangentes tf,δ, el área de la cara Af y el volumen de la celda. Es de vital importancia asegurarse de que los vectores normales son exteriores a la celda para cumplir con la ecuación que se ha obtenido. El área de la cara puede calcularse mediante el Teorema de Pitágoras: Af=q(x2−x1)2+(y2−y1)2(3.15)
44 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado donde los subíndices indican cada uno de los vértices de la cara. Posteriormente, el vector tangente unitario puede derivarse de este valor: tf=1 Afx2−x1 y2−y1(3.16) Las componentes del vector unitario pueden deducirse del sistema de ecuaciones siguiente: n·t=0 n×t=k⇒nf=ty,f −tx,f(3.17) Una manera de comprobar que el vector resultante es exterior consiste en verificar que el producto escalar n·lc f es positivo, y en caso de no serlo, bastaría con cambiar de signo ambos vectores normal y tangente. A continuación se muestra el código que aborda el procesamiento de estos términos anteriormente mencionados. 1%% CALCULO DE AREAS Y NORMALES Y TANGENTES UNITARIAS 2 3lfcx = zeros(n_c,4); 4lfcy = zeros(n_c,4); 5area = zeros(1,n_f); 6tx=zeros(n_c,4); 7ty=zeros(n_c,4); 8nx=zeros(n_c,4); 9ny=zeros(n_c,4); 10 for ic=1:n_c 11 for k=1:4 12 iface = Lcf(ic,k); 13 lfcx(ic,k) = xf(iface) - xc(ic); 14 lfcy(ic,k) = yf(iface) - yc(ic); 15 v1 = Lfv(iface,1); 16 v2 = Lfv(iface,2); 17 dx1 = xv(v2) - xv(v1); 18 dy1 = yv(v2) - yv(v1); 19 area(iface) = sqrt(dx1^2 + dy1^2); 20 tx(ic,k) = dx1/area(iface); 21 ty(ic,k) = dy1/area(iface); 22 nx(ic,k) = ty(ic,k); 23 ny(ic,k) = -tx(ic,k); 24 scalar_nl = nx(ic,k)*lfcx(ic,k) + ny(ic,k)*lfcy(ic,k); 25 if scalar_nl < 0 26 tx(ic,k) = -tx(ic,k); 27 ty(ic,k) = -ty(ic,k); 28 nx(ic,k) = -nx(ic,k); 29 ny(ic,k) = -ny(ic,k); 30 end 31 end 32 end Se recuerda que, de acuerdo con la Figura 3.2 el vector δ se obtiene al realizar el producto escalar n·l , que representa la distancia en dirección normal entre los centroides de las celdas. 1 2%% Cálculo de DELTA 3
3.2 Problema difusivo, estacionario, bidimensional 45 4delta = zeros(1,n_f); 5for iface=1:n_f 6c1 = Lfc(iface,1); 7c2 = Lfc(iface,2); 8for ic=1:n_c 9for k=1:4 10 if Lcf(ic,k) == iface 11 nx1 = nx(ic,k); 12 ny1 = ny(ic,k); 13 end 14 end 15 end 16 17 if c1 == c2 % Cara frontera 18 dx1 = xf(iface) - xc(c1); 19 dy1 = yf(iface) - yc(c1); 20 delta(iface) = abs(dx1*nx1 + dy1*ny1); 21 else 22 dx1 = xc(c2) - xc(c1); 23 dy1 = yc(c2) - yc(c1); 24 delta(iface) = abs(dx1*nx1 + dy1*ny1); 25 end 26 end Por otro lado, el término que multiplica a la generación de calor es el volumen de la celda, este valor se obtiene de la mitad del módulo del producto vectorial de los vectores tangentes de dos de las caras de la celda: ∆V=1 2|t1×t2|=1 2|tx1ty2−ty1tx2|(3.18) Estos vectores deben tener una longitud igual a las caras de la celda. El siguiente código computa el volumen de la celda: 1%% Cálculo del VOLUMEN 2vol(:) = zeros(1,n_c); 3for ic=1:n_c 4v1 = Lcv(ic,1); 5v2 = Lcv(ic,2); 6v3 = Lcv(ic,3); 7v4 = Lcv(ic,4); 8dx1 = xv(v2) - xv(v1); 9dy1 = yv(v2) - yv(v1); 10 vector1 = [dx1; dy1]; 11 dx2 = xv(v4) - xv(v1); 12 dy2 = yv(v4) - yv(v1); 13 vector2 = [dx2; dy2]; 14 semiVol1 = 0.5*abs(dx1*dy2 - dy1*dx2); % semiArea = 1/2*|t1 x t2| 15 16 dx3 = xv(v4) - xv(v3); 17 dy3 = yv(v4) - yv(v3); 18 vector3 = [dx3; dy3]; 19 dx4 = xv(v2) - xv(v3); 20 dy4 = yv(v2) - yv(v3); 21 vector4 = [dx4; dy4];
46 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado 22 semiVol2 = 0.5*abs(dx3*dy4 - dy3*dx4); % semiArea = 1/2*|t1 x t2| 23 24 vol(ic) = semiVol1 + semiVol2; 25 end En la ecuación 3.4 aparece la conductividad térmica en la cara Γf . Normalmente, será necesario expresar alguna propiedad en la cara a partir de su valor en los nodos. Lo más usual, es interpolar linealmente su valor a través de pesos inversamente proporcionales a la distancia: Γf=wfΓ1+(1−wf)Γ2donde wf=1/d1 1/d1+1/d2 ,(3.19) de modo que el valor de Γf se aproximará más a aquella cara cuya distancia sea mas próxima. En el caso extremo que d2⇒∞ , Γf→Γ1 , y lo mismo para el caso opuesto. También es común, aunque menos frecuentre, encontrar interpolaciones que derivan de leyes de la física como una interpolación cuadrática. En los códigos se ha optado por considerar la interpolación lineal siguiento el desarrollo de [ 4 ]. Esta interpolación puede obtenerse del siguiente códidgo: 1wfun = zeros(1,n_f); 2gamma = 0.6*ones(1,n_c); 3gammaf = zeros(1,n_f); 4for iface=1:n_f 5c1 = Lfc(iface); 6c2 = Lfc(iface); 7dx1 = xc(c1) - xf(iface); 8dy1 = yc(c1) - yf(iface); 9d1=sqrt(dx1^2 + dy1^2); 10 dx2 = xc(c2) - xf(iface); 11 dy2 = yc(c2) - yf(iface); 12 d2=sqrt(dx2^2 + dy2^2); 13 wfun(iface) = d2/(d1+d2); 14 gammaf(iface) = wfun(iface)*gamma(c1) + (1-wfun(iface))*gamma(c2); 15 end Asimismo y como puede observarse en la misma ecuación, también es necesario calcular el valor de φ en los vértices, esto puede realizarse siguiendo el mismo procedimiento, con la salvedad de que el número de celdas que intervienen en la interpolación puede ser cualquiera, puesto que un vértice puede estar influenciado por múltiples celdas, de este modo, el valor puede aproximarse como: φv=φi/di ∑N i=11/di = N ∑ i=1 wv,iφidonde wv,i=1/di ∑N i=11/di (3.20) El procedimiento en MATLAB que calcula esta propiedad en los vértices es: 1%% Funciones interpolación Nodo-A-Vértice 2 3dtotal = zeros(n_v,4); 4di=zeros(n_c,4); 5for ic=1:n_c 6for iv=1:4 7v1 = Lcv(ic,iv); 8di(ic,iv) = sqrt((xv(v1) - xc(ic))^2 + (yv(v1) - yc(ic))^2); 9dtotal(v1) = dtotal(v1) + 1/di(ic,iv); 10 end
3.2 Problema difusivo, estacionario, bidimensional 47 11 end 12 cwfun = zeros(n_c,4); 13 for ic=1:n_c 14 for iv=1:4 15 v1 = Lcv(ic,iv); 16 cwfun(ic,iv) = (1/di(ic,iv))/(dtotal(v1)); 17 end 18 end Para terminar con el planteamiento del mallado es necesario tener en cuenta las expresiones en las caras de la frontera, ya que la expresión 3.4 deberá ser adaptada a las condiciones impuestas en dichas caras. Estas expresiones ponen de manifiesto la importancia de poder localizar las caras y vértices de la frontera fácilmente. Para ello, se particularizará primero la ecuación de conservación de la propiedad φ para las dos condiciones de contorno estudiadas en este trabajo. Condición de Dirichlet: temperatura impuesta El valor de φ en la frontera es conocido y el sumatorio se divide entre las caras interiores y la cara exterior. El término (∇φ)f cc ·nf cc no es conocido y debe ser reescrito en términos de φP y φfcc . En caso de que la condición de contorno la temperatura sea constante a lo largo de toda esa frontera, el término tangencial quedaría anulado puesto que φ0 a(fcc)−φ0 b(fcc) sería igual a 0. Además, no existe celda vecina a través de la cara fronteriza, por lo que el vector lf cc debe ser redefinido, en la dirección desde el nodo de la celda P hasta el centro de la cara fcc. Nf,P ∑ f=1 f=f cc Γf"φnb(f)−φP δf − φ0 a(f)−φ0 b(f) δfAf!tf·lf#Af+ Γf cc "φf cc −φP δf cc − φ0 a(f cc)−φ0 b(f cc) δf ccAf cc !tf cc ·lf cc#Af cc+¯ Sφ∆V=0 (3.21) Condición de Neumann: flujo impuesto Ante esta situación, es conocido el flujo de φ en dirección normal (exterior) a la cara. Matemáticamente, se expresa como −(Γ∇φ)fcc ·nfcc =qCC f . De este modo, la expresión resultante queda de la siguiente manera: Nf,P ∑ f=1 f=f cc Γf"φnb(f)−φP δf − φ0 a(f)−φ0 b(f) δfAf!tf·lf#Af+qCC fAf cc +¯ Sφ∆V=0(3.22) El tratamiento de las caras fronterizas requiere una formulación específica, lo cual hace necesario que el programa identifique si las caras bajo análisis pertenecen a la frontera. Hay varias razones por las cuales es preferible almacenar las condiciones de contorno de manera separada. En primer lugar, a medida que la malla se refina, la proporción de celdas fronterizas respecto al total disminuye gradualmente. En segundo lugar, debido a la estructura interna del software, es más eficiente especificar un punto de inicio y una longitud en un arreglo (array) en lugar de definir bloques específicos dentro del mismo. Por último, contar con un vector dedicado para las caras fronterizas ayuda a evitar el uso excesivo de sentencias condicionales ( if ), lo cual puede afectar negativamente la eficiencia del programa. Por estos motivos se formulan las siguientes variables que facilitan la formulación de la ecuación de difusión en las caras de la frontera:
48 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado •Número total de caras en la frontera: nf b •Conectividad Cara-A-CaraFrontera:L f _b f (if) •Conectividad CaraFrontera-A-Cara:Lb f _f(if b) •Vector de valores binarios cuando la cara está en la frontera: b f ace(if) Además, para reconocer si en la cara de la frontera se impone la condición de Dirichlet o de Neumann, se definirán de manera similar los vectores bphiimp_bf , phi_bf y q_bf , donde la primera variable es un vector binario, cuyo papel es el mismo que el de bphiimp_f en el caso de mallas estructuradas de la anterior sección, es decir, bphiimp_f(ifaceb) = 1 si en la cara de la frontera ifaceb se impone temperatura y 0 en si se impone flujo. Por otro lado, las otras dos variables indican qué temperatura o qué flujo se impone las caras si se dan las correspondientes condiciones. Cabe destacar que estos vectores almacenan valores para todas las caras, pero su imposición se rige por la variable bphiimp_bf . A continuación se especifica el código que procesa toda esta información: 1%% Condiciones de contorno 2phi_scc = 40 + 273.15; bphiimp_s = 0; 3phi_ecc = 15 + 273.15; bphiimp_e = 0; 4phi_ncc = 40 + 273.15; bphiimp_n = 0; 5phi_wcc = 30 + 273.15; bphiimp_w = 1; 6 7% variable bfimp_f 1 si se impone el flujo en la cara ’f’, 0 si se impone 8% temperatura 9qwcc = 0; bfimp_w = not(bphiimp_w); 10 qncc = 15; bfimp_n = not(bphiimp_n); 11 qecc = 0; bfimp_e = not(bphiimp_e); 12 qscc = 0; bfimp_s = not(bphiimp_s); 13 14 %% Indexación de Caras Frontera 15 n_f_b = 2*(M+N); % Total en Contorno 16 bface = zeros(1,n_f); % 1 SI cara en Contorno, 0 SI cara Interior 17 f_to_bf = zeros(1,n_f); 18 bf_to_f = zeros(1,n_f_b); 19 bvert = zeros(1,n_v); 20 iface_b = 0; 21 bphiimp_bf = zeros(1,n_f_b); 22 phi_bf = zeros(1,n_f_b); 23 q_bf = zeros(1,n_f_b); 24 % Cara Sur 25 for i=1:M 26 ic = i; 27 iface = ic; 28 iface_b = iface_b + 1; 29 bface(iface) = 1; 30 f_to_bf(iface) = iface_b; 31 bf_to_f(iface_b) = iface; 32 bvert(Lfc(iface,1)) = 1; 33 phi_bf(iface_b) = phi_scc; 34 q_bf(iface_b) = qscc; 35 bphiimp_bf(iface_b) = bphiimp_s; 36 end 37 % Cara Este 38 for j=1:N 39 ic = i + M*(j-1); 40 iface = Lcf(ic,2);
3.2 Problema difusivo, estacionario, bidimensional 49 41 iface_b = iface_b + 1; 42 bface(iface) = 1; 43 f_to_bf(iface) = iface_b; 44 bf_to_f(iface_b) = iface; 45 bvert(Lfv(iface,1)) = 1; 46 phi_bf(iface_b) = phi_ecc; 47 q_bf(iface_b) = qecc; 48 bphiimp_bf(iface_b) = bphiimp_e; 49 end 50 % Cara Norte 51 for i=M:-1:1 52 ic = M*(j-1) + i; 53 iface = Lcf(ic,3); 54 iface_b = iface_b + 1; 55 bface(iface) = 1; 56 f_to_bf(iface) = iface_b; 57 bf_to_f(iface_b) = iface; 58 bvert(Lfv(iface,2)) = 1; 59 phi_bf(iface_b) = phi_ncc; 60 q_bf(iface_b) = qncc; 61 bphiimp_bf(iface_b) = bphiimp_n; 62 end 63 % Cara Oeste 64 for j=N:-1:1 65 ic = M*(j-1) + 1; 66 iface = Lcf(ic,4); 67 iface_b = iface_b + 1; 68 bface(iface) = 1; 69 f_to_bf(iface) = iface_b; 70 bf_to_f(iface_b) = iface; 71 bvert(Lfv(iface,2)) = 1; 72 phi_bf(iface_b) = phi_wcc; 73 q_bf(iface_b) = qwcc; 74 bphiimp_bf(iface_b) = bphiimp_w; 75 end Hasta este punto, se ha generado toda la información inherente al mallado, lo que permite abordar el proceso en el que se formula una ecuación matricial Aφ=b para resolver el problema de difusión, cuestión que se tratará en la siguiente subsección. 3.2.3 Estructuración de la ecuación en forma matricial A la hora de abordar la solucón de la ecuación 3.4 resulta conveniente reestructurarla, primero en términos independientes y no dependientes, y segundo, separar en elementos diagonales y no diagonales, llegando a la siguiente expresión: Nf,P ∑ f=1 ΓfAf δf φp+ Nf,P ∑ f=1 −Γf(P)Af(P) δf(P)!φf(P)=¯ Sφ∆V− Nf,P ∑ f=1 φ0 a(f)−φ0 b(f) δf!tf·lf=bP,(3.23) que es válida para las celdas interiores, y representaría representaría la fila ic(P) de la matriz A . En ella, se han mantenido en un lado de la ecuación los términos que multiplican a la incógnita de modo que el coeficiente del elemento diagonal principal sea positivo, y al otro lado los valores conocidos. Para facilitar la lectura e implementación de la misma, se nombran los coeficientes que multiplican a la incógnita de cada cara y el término de skewness como:
56 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado La comparación entre las soluciones obtenidas con MATLAB y Fluent muestra una coherencia notable en la distribución de las temperaturas a lo largo del dominio. Esta consistencia respalda la fiabilidad del algoritmo implementado en el programa propio. Es importante destacar que el uso de una malla no estructurada de 900 celdas en el programa de CFD demuestra una buena eficiencia computacional, logrando resultados precisos con una resolución relativamente alta. Esto es particularmente relevante en aplicaciones prácticas donde la optimización de recursos computacionales es crucial. 3.2.5 Aplicación práctica 2D. Malla no estructurada triangular Después de analizar y discutir los resultados obtenidos con la malla no estructurada rectangular, se procede al estudio de una malla no estructurada triangular, comúnmente empleada en programas CFD debido a su capacidad para adaptarse a geometrías complejas y su flexibilidad en la discretización del dominio. Para la generación de esta malla triangular, se ha utilizado el software Fluent, que proporciona herramientas avanzadas para la creación y optimización de mallas no estructuradas. Así, para la geometría analizada en difusión, se ha obtenido la siguiente malla: 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figura 3.8 Mallado generado por Fluent. Antes de presentar los resultados obtenidos, se explicará cómo pueden implementarse estas en los códigos de MATLAB. Primero, la malla se exporta como un archivo .STL (por sus siglas en inglés de Standard Tessellation Language). Estos archivos almacenan información sobre la geometría de una superficie en términos de triángulos, describiendo las coordenadas de los vértices y las normales de las caras triangulares. MATLAB dispone de una función para leer estos archivos y almacenarlos en el workspace como un objeto triangulation. Este objeto proporciona los métodos necesarios para manejar toda la conectividad necesaria para el mallado. La siguiente tabla describe algunos de los métodos más relevantes del objeto triangulation en MATLAB que serán empleados en el código: Tabla 3.1 Métodos del objeto triangulation de MATLAB. Método Descripción triangulation.Points Devuelve las coordenadas de los vértices de la malla. triangulation.ConnectivityList Devuelve la lista de conectividad que define cada triángulo mediante índices de vértices. triangulation.incenter Calcula los incentros de los triángulos. triangulation.edges Proporciona los vértices de las caras. triangulation.edgeAttachments Devuelve los triángulos conectados a una arista específica.
3.2 Problema difusivo, estacionario, bidimensional 57 A continuación, se muestran los códigos que obtienen la información del mallado a partir de las funciones que proporciona MATLAB, las cuales se encuentran en la Tabla 3.1. Es importante destacar que el resto del código permanece inalterado con respecto al programa de la malla no estructurada formada por cuadriláteros. 1%% ------------------------ Mallado ------------------------ %% 2 3TR = stlread("malla_elipse.stl"); 4 5% Número total de nodos y Coordenadas de los incentros de las celdas 6incenter = TR.incenter; 7xc_inc = incenter(:,1); 8yc_inc = incenter(:,2); 9n_c = length(incenter); 10 11 % Número total de vértices y Coordenadas de los vértices 12 vertices = TR.Points; 13 xv = vertices(:,1); 14 yv = vertices(:,2); 15 n_v = length(vertices); 16 L=max(xv); 17 idH = find(xv == 0); H = max(yv(idH)); 18 19 % Número total de caras y Centros de las caras 20 % Edges (vert1 --> vert2) 21 edges = TR.edges; 22 v1face = edges(:,1); 23 v2face = edges(:,2); 24 n_f = length(edges); 25 % Faces coordinates 26 xface = [xv(v1face(1:n_f)) xv(v2face(1:n_f))]; 27 yface = [yv(v1face(1:n_f)) yv(v2face(1:n_f))]; 28 29 xf = (xface(:,1) + xface(:,2))/2; 30 yf = (yface(:,1) + yface(:,2))/2; 31 32 %% Conectividades 33 % Conectividad Celda-A-Vértice y Coordenadas de los nodos 34 Lcv = TR.ConnectivityList; 35 xc=zeros(n_c,1); 36 yc=zeros(n_c,1); 37 for ic=1:n_c 38 xc(ic) = sum(xv(Lcv(ic,:)))/3; 39 yc(ic) = sum(yv(Lcv(ic,:)))/3; 40 end 41 42 % Conectividad Cara-A-Vértice 43 Lfv = [v1face v2face]; 44 45 % Conectividad Cara-A-Celda 46 for iface=1:n_f 47 face_con = edgeAttachments(TR,v1face(iface),v2face(iface)); 48 Lfc(iface,1:2) = face_con{:}; 49 end 50 51 % Conectividad Celda-A-Cara 52 for ic=1:n_c 53 k=1;
58 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado 54 for iface=1:n_f 55 if ismember(ic, Lfc(iface,:)) 56 Lcf(ic,k) = iface; 57 k = k+1; 58 end 59 end 60 end Para el mismo problema expuesto en la sección anterior, donde la geometría está conformada por tres caras rectas y una cuarta que sigue la forma de una elipse, se comparan los resultados obtenidos utilizando diferentes configuraciones de mallado no estructurado: rectangular y triangular. (a) Mallado no estructurado rectangular. (b) Mallado no estructurado triangular. Figura 3.9 Comparación de soluciones para dos mallados no estructurados. Como puede observarse, la solución obtenida por el mallado no estructurado triangular presenta unos valores muy similares con respecto a la malla rectangular. Las mallas triangulares son más representativas en los programas de CFD por el buen carácter adaptativo a geometrías irregulares. Regularmente, las mallas que se generan en estos programas combinan ambos tipos de celdas para conseguir la mayor calidad posible en la malla. 3.3 Problema de convección-difusión térmica estacionario y bidimensional Para completar la formulación del FVM en mallas no estructuradas, esta sección se centra en el análisis de problemas convectivo-difusivos. Como se vió en capítulos anteriores, la ecuación que describe la transferencia de calor en procesos estacionarios convectivo-difusivos con generación de calor es la ecuación de la energía, y obviando los términos transitorios y de viscosidad, se obtiene: ∇·(ρCvvφ) = ∇·(Γ∇φ)+Sφ(3.28) Procediendo con el paso fundamental del FVM, la integración del término convectivo para un volumen de control de normal exterior n y su posterior aplicación del Teorema de la Divergencia de Gauss llevan al siguiente resultado: ZVC ∇·(ρCvvφ)dV =ZS (ρCvvφ)·ndA = Nf,P ∑ f=1 (ρv·nA)fφf= Nf,P ∑ f=1 (ρvnA)fφf= Nf,P ∑ f=1 Gfφf(3.29)
3.3 Problema de convección-difusión térmica estacionario y bidimensional 59 Donde se ha supuesto que las celdas se componen de segmentos rectos. En esta ecuación se ha definido el gasto másico Gf= (ρvnA)f , que indica la masa de fluido que atraviesa la cara de celda en la dirección normal exterior. Esta ecuación pone de manifiesto la necesidad de calcular la velocidad del fluido en dirección normal exterior a la cara f cuando pertenece a una celda dada. Al igual que en el capítulo anterior, el campo de velocidades puede determinarse resolviendo la ecuación de continuidad. Según la Teoría Potencial para flujos incompresibles, esto se logra a partir de la ecuación de Laplace para una función de corriente φ , aplicando las condiciones de contorno adecuadas: ∇2φ=0,(3.30) donde las velocidades se expresan como función de las derivadas espaciales de esta función: u=∂φ ∂yv=−∂φ ∂x(3.31) Las condiciones de contorno que se aplican a la geometría para la función de corriente φ en cada una de las caras son: φW(y) = uey(3.32) φE(y) = usy(3.33) φN=ueHe(3.34) φS=0(3.35) La dificultad de calcular las velocidades estriba en que las derivadas 3.31 no pueden calcularse directamente mediante expansiones en series de Taylor, ya que las celdas no se encuentran alineadas con los ejes x,y . Sin embargo, un procedimiento similar al desarrollo del gradiente de φ al realizado para el mallado no estructurado puede emplearse para obtener las componentes de la velocidad en los ejes cartesianos. Puede comenzarse de nuevo teniendo en cuenta que cualquier vector puede expresarse en cualquier sistema de coordenadas como se mostró en la ecuación 3.5: ∇φ=∂φ ∂xi+∂φ ∂yj=∂φ ∂nn+∂φ ∂tt(3.5) A continuación, considerando que las velocidades 3.31 pueden expresarse como el producto vectorial ∇φ×k , dicho producto puede realizarse en el sistema de coordenadas formado por los vectores normal y tangente a la cara y el vector k: v=∇φ×k= i j k φxφy0 0 0 1 = n t k φnφt0 0 0 1 =∂φ ∂tn−∂φ ∂nt(3.36) Teniendo en cuenta nuevamente que (∇φ)·n=∂φ ∂n y (∇φ)·t=∂φ ∂t , y expresando los vectores normal y tangente en sus componentes cartesianas, estas pueden sumarse, obteniendo: v=∇φ·t(nxi+nyj)−∇φ·n(txi+tyj) = (nx∇φ·t−tx∇φ·n)i+(ny∇φ·t−ty∇φ·n)j,(3.37) y puesto que previamente se han discretizado los productos escalares ∇φ·n y ∇φ·t , resultados expresados
60 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado en 3.13 y 3.11 respectivamente, puede implementarse en el código sin mayor dificultad. Para ello, únicamente habrá que añadir al programa el siguiente bloque de código en la sección posterior al bucle iterativo donde se calcula la solución: 1%% Cálculo de la velocidad 2 3phiv = zeros(1,n_v); 4for ic=1:n_c 5for k=1:3 6iv = Lcv(ic,k); 7if bvert(iv) == 0 8phiv(iv) = phiv(iv) + phi(ic)*cwfun(ic,k); 9else 10 phiv(iv) = phi_bv(iv); 11 end 12 end 13 end 14 15 gphi_t = zeros(n_f,1); 16 gphi_n = zeros(n_f,1); 17 vx_f = zeros(n_f,1); 18 vy_f = zeros(n_f,1); 19 signo = zeros(n_c,3); 20 vn=zeros(n_c,3); 21 22 for ic=1:n_c 23 for k=1:3 24 iface = Lcf(ic,k); 25 va = Lfv(iface,1); 26 vb = Lfv(iface,2); 27 ic12 = Lfc(iface,:); 28 ic2 = ic12(ic12~=ic); 29 signo(ic,k) = ((xv(va)-xv(vb))*tx(ic,k) + (yv(va)-yv(vb))*ty(ic,k))/(area(iface) ); 30 if not(bface(iface)) 31 tdotl = signo(ic,k)*((xv(va)-xv(vb))*(xc(ic2)-xc(ic)) + (yv(va)-yv(vb))*(yc( ic2)-yc(ic))); 32 gphi_t(iface) = signo(ic,k)*(phiv(va)-phiv(vb))/area(iface); 33 gphi_n(iface) = (phi(ic2)-phi(ic))/delta(iface) - (signo(ic,k)*(phiv(va)-phiv (vb)))/area(iface)*tdotl/(area(iface)*delta(iface)); 34 else 35 tdotl = signo(ic,k)*((xv(va)-xv(vb))*(xf(iface)-xc(ic)) + (yv(va)-yv(vb))*(yf (iface)-yc(ic))); 36 gphi_t(iface) = signo(ic,k)*(phiv(va)-phiv(vb))/area(iface); 37 gphi_n(iface) = (phi_bf(f_to_bf(iface))-phi(ic))/delta(iface) - (signo(ic,k) *(phiv(va)-phiv(vb)))*tdotl/(area(iface)^2*delta(iface)); 38 end 39 vx_f(iface) = gphi_t(iface)*nx(ic,k) - gphi_n(iface)*tx(ic,k); 40 vy_f(iface) = gphi_t(iface)*ny(ic,k) - gphi_n(iface)*ty(ic,k); 41 vn(ic,k) = gphi_t(iface); 42 end 43 end 44 45 vx_c = zeros(n_c,1); 46 vy_c = zeros(n_c,1); 47 for ic=1:n_c 48 d = 1./sqrt((xc(ic) - xf(Lcf(ic,:))).^2 + (yc(ic) - yf(Lcf(ic,:))).^2);
3.3 Problema de convección-difusión térmica estacionario y bidimensional 61 49 vx_c(ic) = sum(vx_f(Lcf(ic,:)).*d./sum(d)); 50 vy_c(ic) = sum(vy_f(Lcf(ic,:)).*d./sum(d)); 51 end En este bloque, primero se actualizan los valores de φ en los vértices acorde a la solución final, a continuación se entra en un bucle que calcula las componentes de velocidad de cada cara y la velocidad vn=v·n en base a los productos escalares ∇φ·n y ∇φ·t . Por último, aproxima el valor de la velocidad en los nodos mediante la interpolación basada en los pesos inversamente proporcionales a la distancia. Para comprobar que las velocidades obtenidas son las correctas, puede implementarse el siguiente código para visualizar el campo escalar φy el campo de velocidades v. 1Mv = 8; 2Nv = 6; 3xrepr = zeros(Mv*Nv,1); 4yrepr = zeros(Mv*Nv,1); 5ic_list = zeros(Mv*Nv,1); 6for j=1:Nv 7for i=1:Mv 8k = (j-1)*Mv+i; 9xrepr(k) = L/(Mv+1)/2 + (i-1)*(L/(Mv+1)); 10 yrepr(k) = S2(xrepr(k))/(Nv)/2 + (j-1)*S2(xrepr(k))/Nv; 11 end 12 end 13 for j=1:Nv 14 for i=1:Mv 15 k = (j-1)*Mv+i; 16 [~, ic_list(k)] = min(sqrt((xc-xrepr(k)).^2 + (yc-yrepr(k)).^2)); 17 end 18 end 19 figure(5) 20 hold on 21 axis([-0.1*L L*1.1 -0.1*H H*1.1]); 22 axis(’equal’) 23 tckness = 1.5; 24 linecolor = ’#CE9FFF’; 25 plot([0 L], [0 0], ’Color’,linecolor,’LineWidth’,tckness); 26 plot([0 0],[0 H1], ’Color’,linecolor,’LineWidth’,tckness); 27 plot([0:0.01:L], S2([0:0.01:L]), ’Color’,linecolor,’LineWidth’,tckness); 28 plot([L L], [0 H2], ’Color’,linecolor,’LineWidth’,tckness); 29 title(’Campo de velocidades’,’Interpreter’,’latex’) 30 quiver(xrepr, yrepr, vx_c(ic_list), vy_c(ic_list), ’Color’,’#388ADC’,’linewidth’,1.2); Resolviendo este problema con la formulación no estructurada del FVM en la geometría de cara norte exponencial sobre la malla de la Figura 3.10, y para una velocidad de entrada ve=0,01 m/s , se obtiene el campo de velocidades y la función de corriente de las Figuras 3.11 y 3.12. Nuevamente puede observarse la concordancia entre los resultados obtenidos con mallas estructuradas y los obtenidos con mallas no estructuradas. De igual modo que con el programa de flujo potencial en mallas estructuradas, este puede funcionar también como una función que proporciona al programa principal la velocidad normal a la cara vn=v·nf . El siguiente paso en la formulación del problema convectivo-difusivo consiste en computar el valor del gasto a través de la cara Gf. Esto puede llevarse a cabo en un par de líneas de código: 1vn = Flujo_Potencial_UNST_ext_mesh; 2Gf = rho*Cv*vn.*area(Lcf);
62 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado Figura 3.10 Mallado no estructurado generado por Fluent. Figura 3.11 Campo de velocidades en una malla no estructurada. Figura 3.12 Función de corriente. Es pertinente recordar que la variable vn tiene dimensiones nc×3 , donde nc representa el número de celdas y 3 el número de caras por celda. Las constantes ρ y Cv están definidas en el script de configuración. Una vez calculado el gasto másico, solo resta expresar el valor de φf en función de los valores en los nodos. La
3.3 Problema de convección-difusión térmica estacionario y bidimensional 63 forma óptima de lograrlo ha sido objeto de una investigación extensa, y varios métodos han sido analizados en el trabajo [ 4 ], al cual se remite al lector para una discusión más detallada. Por consiguiente, se presentarán únicamente los métodos de diferencias centradas y el esquema upwind, dado que son los más utilizados en el ámbito de CFD. En el método de diferencias centradas, φf puede expresarse en función de los valores de φ en las celdas adyacentes a la cara f , utilizando la interpolación basada en los pesos wf inversamente proporcionales a la distancia, tal como se presentó en la Sección 3.2.2: Nf,P ∑ f=1 Gfφf= Nf,P ∑ f=1 Gf[wfφP+(1−wf)φnb(f)](3.38) Puede demostrarse mediante expansiones en series de Taylor que el error cometido en esta interpolación es del orden del cuadrado de la distancia de la cara al centro de la celda, al menos cuando esta distancia es la misma a ambas celdas. Esta característica hace a este método atractivo desde un punto de vista matemático. Sin embargo, como se mencionó en la Sección 3.2.2, en convección, la cantidad φ se transporta en la dirección del flujo del fluido, por lo que el valor de φf será más similar al de la celda situada aguas arriba. Como es sabido, el esquema upwind de primer orden aproxima el valor φf como aquel de la celda adyacente que se encuentra aguas arriba de la cara. Esta idea puede exponerse matemáticamente para celdas de una malla no estructurada mediante la siguiente expresión: φf=φOsi Gf>0 φnb si Gf<0(3.39) La manera de aplicar esta condición de manera matemática para evitar el costoso uso de consecutivas comprobaciones IF, puede expresarse como: Nf,P ∑ f=1 Gfφf= Nf,P ∑ f=1Gf+|Gf| 2φO+Gf−|Gf| 2φnb(f)(3.40) La ventaja de este otro método frente al de las diferencias centradas, es que representa más fielmente la naturaleza física de la convección, sin embargo, matemáticamente el error decae con la distancia en lugar de con la distancia al cuadrado como lo hacía con las diferencias centradas. Además, este esquema puede llevar asociado una falsa difusión, que aparece en problemas puramente convectivos cuando la dirección del flujo no se alinea con los ejes coordenados. A pesar de ello, este es un modelo perfectamente válido para los casos aquí tratados. Para terminar, añadiendo el término convectivo a la ecuación 3.23, se formula la ecuación general que modela los problemas convectivo-difusivos estacionarios sin viscosidad en mallas no estructuradas: Nf,P ∑ f=1Gf+|Gf| 2φO+Gf−|Gf| 2φnb(f)+ Nf,P ∑ f=1 ΓfAf δf φp+ Nf,P ∑ f=1 −Γf(P)Af(P) δf(P)!φf(P)=¯ Sφ∆V− Nf(P) ∑ f=1 φ0 a(f)−φ0 b(f) δf!tf·lf=bP (3.41) La ecuación puede formularse en formato matricial a través de los coeficientes Ap,P , elemento diagonal, Ak,P , el elemento correspondiente a la celda que comparte con P la cara k de la celda P y Askew,P , que dado que se ha optado por el método explícito, este valor es conocido y puede almacenarse en el vector de términos independientes. Así, los coeficientes de la matriz son:
64 Capítulo 3. FVM para problemas de convección-difusión térmica. Mallado no estructurado Ap= Nf,P ∑ f=1 ΓfAf δf +Gf+|Gf| 2!(3.42) Af(P)=−Γf(P)Af(P) δf(P) −Gf+|Gf| 2(3.43) Sskew,P=− Nf(P) ∑ f=1 φ0 a(f)−φ0 b(f) δf!tf·lf(3.44) Es notable observar que el vector de términos independientes permanece inalterado. Esto se debe a que el fenómeno convectivo es proporcional a φ y no al gradiente ∇φ . Por consiguiente, la discretización no necesita utilizar el valor de φ en los nodos, por lo que su implementación en el código es bastante directa, y únicamente afecta a la definición de la matriz A: 1for ic=1:n_c 2for k=1:3 3iface = Lcf(ic,k); 4if bface(iface) == 0 % Cara interior 5icvec = Lfc(iface,:); 6A(ic,icvec(icvec ~= ic)) = -gammaf(iface)*area(iface)/delta(iface) + 0.5*(Gf( ic,k) - abs(Gf(ic,k))); 7A(ic,ic) = A(ic,ic) + 0.5*(Gf(ic,k) + abs(Gf(ic,k))); 8elseif bfimp_f(iface) == 1 % Flujo impuesto 9A(ic,ic) = A(ic,ic) + Gf(ic,k); 10 elseif bphiimp_bf(f_to_bf(iface)) == 1 % Temperatura impuesta 11 b(ic) = b(ic) - Gf(ic,k)*phi_bf(f_to_bf(iface)); 12 end 13 end 14 A(ic,ic) = A(ic,ic) + sum(gammaf(Lcf(ic,:)).*area(Lcf(ic,:))./delta(Lcf(ic,:)).*not( bfimp_f(Lcf(ic,:)))); 15 end Una vez incorporadas estas líneas, se está en condiciones de resolver el problema de convección-difusión. Para ello, téngase en cuenta el problema definido en la Sección 2.3.3, donde el fluido entra por la cara oeste a una temperatura de 30 ◦ C, en la cara sur hay un flujo de calor entrante de 7 W/m2 , y el resto se encuentran aisladas térmicamente. La malla sobre la que se resuelve el problema es la mostrada en la Figura 3.10. (a) Solución de MATLAB. (b) Solución de Fluent para el problema convectivo. Figura 3.13 Comparación de soluciones para ve=10−6m/s.
3.3 Problema de convección-difusión térmica estacionario y bidimensional 65 (a) Solución de MATLAB. (b) Solución de Fluent para el problema convectivo. Figura 3.14 Comparación de soluciones para ve=10−4m/s. (a) Solución de MATLAB. (b) Solución de Fluent para el problema convectivo. Figura 3.15 Comparación de soluciones para ve=10−2m/s. Llegados a este punto, se puede concluir con certeza que los códigos desarrollados para la formulación del FVM en problemas de convección-difusión térmica en geometrías variables han logrado proporcionar resultados satisfactorios. Esta conclusión se valida por la coherencia observada entre los resultados obtenidos con mallas estructuradas, no estructuradas, y los resultados proporcionados por el programa Fluent. El estudio de la generación de mallas no estructuradas no se aborda en este trabajo, ya que, como se ha mencionado previamente, su investigación puede requerir un esfuerzo tan extenso como el de la formulación del FVM. Para aquellos lectores interesados en los procedimientos empleados en la práctica para la generación de mallas tanto estructuradas como no estructuradas, se recomienda consultar los artículos de J. F. Thompson [6], N. P. Weatherhill [14], la obra de P. J. Frey [7], o el informe técnico de D. J. Mavriplis [15]. En el capítulo último, se realiza una síntesis de todo el trabajo llevado a cabo en esta memoria, además de presentar una lista de posibles investigaciones futuras que enriquecerían significativamente el contenido tratado en este proyecto. Finalmente, se han añadido dos apéndices que complementan la memoria. El primero detalla cómo se han abordado los problemas de convección-difusión mediante la ecuación de la energía en el software Fluent, describiendo la definición del dominio, el proceso de mallado y la obtención de las soluciones. El segundo apéndice incluye los códigos completos que generan las soluciones presentadas a lo largo de la memoria, evitando así redundancias en los capítulos principales.
72 Capítulo A. Solución problema difusivo con Fluent Figura A.4 Creación de una superficie a través de bordes. geométricos (como tetraedros o hexaedros en 3D) que permiten la discretización de las ecuaciones. Para iniciar el mallado en Fluent, se accede mediante clic derecho en la pestaña de Mesh y seleccionando Edit, donde se ajustan parámetros como la densidad y el tipo de elementos según las necesidades específicas del problema. Una vez configurado, al activar la función Generate o Mesh > Generate Mesh, Fluent genera automáticamente un mallado utilizando elementos tetraédricos para geometrías tridimensionales y elementos triangulares para geometrías bidimensionales. Este proceso se adapta dinámicamente según la complejidad del modelo, asegurando una transición suave entre las superficies y el volumen, y aplicando criterios de calidad para mantener la precisión de la simulación. Figura A.5 Adicionalmente, Fluent permite la creación de mallados estructurados mediante Face Meshing desde el menú Mesh. Aquí, se puede seleccionar entre varias alternativas como la generación de celdas cuadriláteras para una superficie. Para definir el número exacto de elementos a lo largo de cada eje, se utiliza la función de sizing dentro del mismo menú, donde se especifican las caras a afectar y se elige el tamaño deseado de los elementos o el número de divisiones a aplicar a lo largo de dichas caras. Es importante señalara que para que la rutina sizing funcione correctamente hay que indicar que las líneas que se han seleccionado son las
73 interiores. Al seleccionar una de estas líneas, en la esquina inferior izquierda de la ventana de la geometría aparecerán dos rectángulos indicando el interior y el exterior como los que aparecen en la figura A.5. Figura A.6 Como última etapa del mallado, es necesario indicar las caras que pertenecen a la frontera para posteriormente especificar las condiciones de contorno. Esto se hace a través de la opción seleccionar aristas en la barra de herramientas superior y con el botón derecho sobre la cara previamente seleccionada. Así, aparecerá un menú en el que se encuentra la opción de Create Named Selection que tras pulsar permite indicarle un nombre. Nuevamente, es necesario indicar que se trata de la línea de contorno de la superficie interior. Figura A.7 Para finalizar, es necesario comunicar la información de la malla a la solución para que pueda procesarla. Esta acción se realiza mediante la opción Update, disponible en el menú de Mesh.
74 Capítulo A. Solución problema difusivo con Fluent Solución del problema. Materiales y condiciones de contorno Para abordar el problema de difusión en ANSYS Fluent, se empleará un enfoque que integra varios modelos físicos y configuraciones esenciales para una simulación precisa. Primero, se utilizará el modelo de flujo laminar para capturar el comportamiento fluido de manera estable y predecible cuando se estudie el problema de convección-difusión. Este modelo se complementará con la inclusión de la ecuación de la energía, permitiendo analizar la transferencia térmica dentro del sistema. Estos modelos se encuentran en la columna izquierda Outline View dentro de Models. Figura A.8 El material seleccionado para este estudio será el agua, caracterizado por sus propiedades termodinámicas específicas que influirán en la distribución de la temperatura y la difusión de calor. Por defecto, este material no se encuentra cargado automáticamente, pero sí está disponible en la base de datos de Fluent y se puede añadir con el botón copy. Posteriormente, hay que indicar en las Cell Zone Conditions que el fluido del que se trata es el agua. Figura A.9
75 En cuanto a las condiciones de contorno, se configurarán como paredes para representar las superficies sólidas del dominio. Estas condiciones serán fundamentales para simular cómo el agua interactúa térmicamente con las estructuras circundantes. Figura A.10 El siguiente paso es realizar la inicialización en Fluent, que implica establecer condiciones iniciales precisas para variables físicas como velocidad, presión y temperatura antes de iniciar la simulación. Este paso es crucial para asegurar la convergencia del solver numérico y obtener resultados confiables. Para este propósito, se ha elegido una inicialización estándar que incluye la opción de calcular desde una cara específica. En problemas de convección, es beneficioso seleccionar la cara de entrada, ya que permite que Fluent determine automáticamente las condiciones iniciales basadas en las condiciones de entrada definidas en esa cara. El último paso consiste en ejecutar los cálculos en Fluent y examinar los resultados. La temperatura y otras variables pueden visualizarse utilizando los contornos (contours), lo cual permite analizar la distribución espacial y los gradientes dentro del dominio de la simulación. Además, Fluent ofrece una gráfica que muestra el error de cada iteración, facilitando la evaluación del proceso de convergencia del solver numérico y asegurando la precisión de los resultados obtenidos.
76 Capítulo A. Solución problema difusivo con Fluent Figura A.11
Apéndice B Códigos de MATLAB En este apéndice se incluyen los códigos de MATLAB más representativos desarrollados por el autor. Estos códigos han sido fundamentales para la implementación y validación de los métodos y algoritmos descritos a lo largo del presente trabajo. Se espera que estos códigos puedan servir como referencia y guía para futuros trabajos en el área. B.1 Mallado estructurado B.1.1 Script configuración problema difusivo y convectivo-difusivo 1% Parámetros del problema FVM. 2 3%% Mallado 4M=100; 5N=100; 6H=2; 7L=2; 8 9%% Datos del problema 10 % Gamma = 0.0242; % Conductividad Térmica (W/(m*K)) -> CTE Agua = 0.606 W/(m*K); Aire = 0.02572 W/(m*K); 11 Gamma = 0.6; 12 Cv = 720; % J/(K*kg) Agua = 4182 J/(K*kg); Aire = 718 J/(K*kg); 13 % Cv = 4182; 14 rho = 1.225; % kg/m^3 Agua = 998.2 kg/m^3; Aire = 1.225 kg/m^3; 15 % rho = 998.2; 16 v1 = 0.000; % m/s 17 18 nu = 1.5e-5; % Viscosidad cinemática del aire = 1.5e-5; agua = 1e-6; 19 20 %% Condiciones de contorno 21 % Condiciones de contorno de temperatura 22 % variable bphiimp_x 1 si se impone temperatura en la cara ’x’, 0 si se impone 23 % flujo 24 phi_scc = 40 + 273.15; bphiimp_s = 0; 25 phi_ecc = 15 + 273.15; bphiimp_e = 0; 26 phi_ncc = 40 + 273.15; bphiimp_n = 0; 27 phi_wcc = 30 + 273.15; bphiimp_w = 1; 28 29 % variable bfimp_x 1 si se impone el flujo en la cara ’x’, 0 si se impone 77
78 Capítulo B. Códigos de MATLAB 30 % temperatura 31 qwcc = 0; bfimp_w = not(bphiimp_w); 32 qncc = -15; bfimp_n = not(bphiimp_n); 33 qecc = 0; bfimp_e = not(bphiimp_e); 34 qscc = 0; bfimp_s = not(bphiimp_s); 35 36 Nx = M+1; % Número de vértices en x 37 Ny = N+1; % Número de vértices en y 38 39 40 %% Definición de la forma 41 % Función superior 42 43 % Línea horizontal 44 % ys = @(x) H; 45 % yp = @(x) 0; 46 47 % Exponencial 48 % sB=0.75; a=(sB-H)/(exp(-1)-1); 49 % ys = @(x) a*exp(-x/L) + H - a; 50 % yp = @(x) -a*exp(-x/L)/L; 51 52 % Elipse 53 a = L*1.2; b = H; 54 ys = @(x) b*sqrt(1-x.^2/a^2); 55 yp = @(x) -b*x/a^2/sqrt(1-x.^2/a^2); 56 57 % Circunferencia (Pendiente infinita --> da error) 58 % x0 = 1; y0 = 2; R=1; 59 % ys = @(x) y0 - sqrt(R^2 - (x - x0).^2)/4; 60 % yp = @(x) -(x - x0)./sqrt(R^2 - (x - x0).^2)/4; 61 62 % Parábola 63 % x0 = 1; y0 = 1.7; a=0.3; 64 % ys = @(x) a*(x-x0).^2 + y0; 65 % yp = @(x) 2*a*(x-x0); 66 67 normal_nx = @(xi) -yp(xi*L)/sqrt(1 + yp(xi*L)^2); 68 normal_ny = @(xi) 1/sqrt(1 + yp(xi*L)^2); 69 70 %% Jacobiano 71 J = @(xi,eta) [L 0; eta*yp(xi*L)*L ys(xi*L)]; 72 detJ = @(xi,eta) L*ys(xi*L); 73 74 beta11 = @(xi,eta) ys(xi*L); 75 beta12 = @(xi,eta) -eta.*yp(xi*L)*L; 76 beta21 = @(xi,eta) 0; 77 beta22 = @(xi,eta) L; 78 79 B11 = @(xi,eta) beta11(xi,eta)^2 + beta21(xi,eta)^2; 80 B12 = @(xi,eta) beta12(xi,eta)*beta11(xi,eta)+beta22(xi,eta)*beta21(xi,eta); 81 B21 = @(xi,eta) beta11(xi,eta)*beta12(xi,eta)+beta21(xi,eta)*beta22(xi,eta); 82 B22 = @(xi,eta) beta12(xi,eta)^2+beta22(xi,eta)^2;
B.1 Mallado estructurado 79 B.1.2 Función que resuelve el campo potencial 1function [vx_c, vy_c, vx_f, vy_f] = Flujo_Potencial 2%% Convección pura en 2D en geometría variable dada por ys(x) 3% close all; clear; clc 4 5% Ejecuta el script de parámetros 6config_conv; 7 8%% Distribución xi y eta. 9xi=zeros(M,N); eta = zeros(M,N); 10 yvert = zeros(M,N); 11 for j=1:Ny 12 xi(1:Nx, j) = linspace(0,1,Nx); 13 end 14 for i=1:Nx 15 eta(i,1:Ny) = linspace(0,1,Ny); 16 yvert(i,1:Ny) = linspace(0,ys(xi(i)*L),Ny); 17 end 18 xvert = xi*L; 19 20 %% FVM 21 % COORDENADAS DE LOS VERTICES, CARAS Y NODOS 22 n_c = M*N; % Número total de caras 23 24 xi1 = reshape(xi(1:M,1:N), 1, n_c); eta1 = reshape(eta(1:M,1:N), 1, n_c); 25 xi2 = reshape(xi(1:M,2:N+1), 1, n_c); eta2 = reshape(eta(1:M,2:N+1), 1, n_c); 26 xi3 = reshape(xi(2:M+1,2:N+1), 1, n_c); eta3 = reshape(eta(2:M+1,2:N+1), 1, n_c); 27 xi4 = reshape(xi(2:M+1,1:N), 1, n_c); eta4 = reshape(eta(2:M+1,1:N), 1, n_c); 28 % Nodos 29 xinodo = (xi1+xi4)/2; 30 etanodo = (eta1 + eta2)/2; 31 % Caras 32 xie = xi4; etae = etanodo; 33 xin = xinodo; etan = eta2; 34 xiw = xi1; etaw = etanodo; 35 xis = xinodo; etas = eta1; 36 37 % Vectores canónicos globales 38 39 phi_SW = zeros(1,2*M+4); phi_SW(1,1) = 1; 40 phi_S = zeros(1,2*M+4); phi_S(1,2) = 1; 41 phi_SE = zeros(1,2*M+4); phi_SE(1,3) = 1; 42 phi_W = zeros(1,2*M+4); phi_W(1,M+1) = 1; 43 phi_O = zeros(1,2*M+4); phi_O(1,M+2) = 1; 44 phi_E = zeros(1,2*M+4); phi_E(1,M+3) = 1; 45 phi_NW = zeros(1,2*M+4); phi_NW(1,2*M+1) = 1; 46 phi_N = zeros(1,2*M+4); phi_N(1,2*M+2) = 1; 47 phi_NE = zeros(1,2*M+4); phi_NE(1,2*M+3) = 1; 48 b_v = [zeros(1,2*M+3) -1]; 49 50 istepe = zeros(M,N); istepe(M,:) = 1; stepe = not(istepe); 51 isteps = zeros(M,N); isteps(:,1) = 1; steps = not(isteps); 52 istepw = zeros(M,N); istepw(1,:) = 1; stepw = not(istepw); 53 istepn = zeros(M,N); istepn(:,N) = 1; stepn = not(istepn); 54
80 Capítulo B. Códigos de MATLAB 55 % Datos del flujo 56 57 H1 = H; 58 H2 = ys(L); 59 v2 = v1*H1/H2; 60 61 deltaxi = 1/M; 62 deltaeta = 1/N; 63 A=sparse(n_c,n_c); 64 b=sparse(1,n_c); 65 phi_scc = 0; 66 phi_ncc = v1*H1; 67 phi_ecc = @(eta) v2*H2*eta; 68 dphi_e_deta = @(eta) v2*H2; 69 phi_wcc = @(eta) v1*H1*eta; 70 dphi_w_deta = @(eta) v1*H1; 71 72 % Formulación de la Matriz A 73 for j=1:N 74 for i=1:M 75 k = (j-1)*M+i; 76 dphidxi_e = (phi_E - phi_O)/deltaxi*stepe(i,j) + (phi_W + 8*phi_ecc(etae(k))*b_v - 9*phi_O)/(3*deltaxi)*istepe(i,j); 77 dphideta_e = (phi_NE + phi_N - phi_SE - phi_S)/(4*deltaeta)*stepe(i,j)*stepn(i,j )*steps(i,j) ... 78 + (phi_ncc*b_v-(phi_S + phi_SE + phi_O + phi_E)/4)/deltaeta*istepn(i,j)*stepe (i,j) ... 79 + ((phi_N + phi_NE + phi_O + phi_E)/4 - phi_scc*b_v)/(deltaeta)*isteps(i,j)* stepe(i,j) ... 80 + dphi_e_deta(etae(k))*istepe(i,j)*b_v; 81 dphidxi_w = (phi_O - phi_W)/deltaxi*stepw(i,j) + (9*phi_O - phi_E - 8*phi_wcc( etaw(k))*b_v)/(3*deltaxi)*istepw(i,j); 82 dphideta_w = (phi_NW + phi_N - phi_SW - phi_S)*steps(i,j)/(4*deltaeta)*stepw(i,j )*stepn(i,j)*steps(i,j) ... 83 + (phi_ncc*b_v-(phi_S + phi_SW + phi_O + phi_W)/4)/deltaeta*istepn(i,j)*stepw (i,j) ... 84 + ((phi_N + phi_NW + phi_O + phi_W)/4 - phi_scc*b_v)/(deltaeta)*isteps(i,j)* stepw(i,j) ... 85 + dphi_w_deta(etaw(k))*istepw(i,j)*b_v; 86 dphidxi_s = (phi_E + phi_SE - phi_SW - phi_W)/(4*deltaeta)*steps(i,j)*stepe(i,j) *stepw(i,j) ... 87 + (phi_ecc(etas(k))*b_v-(phi_SW + phi_W + phi_O + phi_S)/4)/(deltaeta)*istepe (i,j)*steps(i,j) ... 88 + ((phi_SE + phi_E + phi_O + phi_S)/4 - phi_wcc(etas(k))*b_v)/(deltaeta)* istepw(i,j)*steps(i,j); 89 dphideta_s = (phi_O - phi_S)/deltaxi*steps(i,j) + (9*phi_O - phi_N - 8*phi_scc* b_v)/(3*deltaeta)*isteps(i,j); 90 dphidxi_n = (phi_NE + phi_E - phi_NW - phi_W)/(4*deltaeta)*stepn(i,j)*stepe(i,j) *stepw(i,j) ... 91 + (phi_ecc(etan(k))*b_v-(phi_NW + phi_W + phi_O + phi_N)/4)/(deltaeta)*istepe (i,j)*stepn(i,j) ... 92 + ((phi_NE + phi_E + phi_O + phi_N)/4 - phi_wcc(etan(k))*b_v)/(deltaeta)* istepw(i,j)*stepn(i,j); 93 dphideta_n = (phi_N - phi_O)/deltaxi*stepn(i,j) + (phi_S + 8*phi_ncc*b_v - 9* phi_O)/(3*deltaeta)*istepn(i,j); 94 95 qEcc = (beta11(xie(k),etae(k))*dphidxi_e + beta12(xie(k),etae(k))*dphideta_e);
B.1 Mallado estructurado 81 96 qE2 = (beta21(xie(k),etae(k))*dphidxi_e + beta22(xie(k),etae(k))*dphideta_e); 97 qWcc = - (beta11(xiw(k),etaw(k))*dphidxi_w + beta12(xiw(k),etaw(k))*dphideta_w); 98 qW2 = - (beta21(xiw(k),etaw(k))*dphidxi_w + beta22(xiw(k),etaw(k))*dphideta_w); 99 qNcc = (beta21(xin(k),etan(k))*dphidxi_n + beta22(xin(k),etan(k))*dphideta_n); 100 qN2 = (beta11(xin(k),etan(k))*dphidxi_n + beta12(xin(k),etan(k))*dphideta_n); 101 qScc = - (beta21(xis(k),etas(k))*dphidxi_s + beta22(xis(k),etas(k))*dphideta_s); 102 qS2 = - (beta11(xis(k),etas(k))*dphidxi_s + beta12(xis(k),etas(k))*dphideta_s); 103 104 eq = ( - qEcc*beta11(xie(k),etae(k)) - qE2*beta21(xie(k),etae(k)))*deltaeta ... 105 - ( qWcc*beta11(xiw(k),etaw(k)) + qW2*beta21(xiw(k),etaw(k)))*deltaeta ... 106 + ( - qNcc*beta22(xin(k),etan(k)) - qN2*beta12(xin(k),etan(k)))*deltaxi ... 107 - ( qScc*beta22(xis(k),etas(k)) + qS2*beta12(xis(k),etas(k)))*deltaxi; 108 109 index = [k-M-1 k-M k-M+1 zeros(1,k-1-(k-M+1)-1) k-1 k k+1 zeros(1,k-1-(k-M+1)-1) k+M-1 k+M k+M+1]; 110 starts_index = index(index > 0 & index <= M*N); 111 starts = starts_index(1); ends = starts_index(end); 112 113 eq_start = find(index == starts_index(1)); 114 eq_end = find(index == starts_index(end)); 115 116 A(k, starts:ends) = eq(eq_start:eq_end); 117 b(k) = eq(end); 118 end 119 end 120 121 phi1D = A\b’; 122 k=reshape((0:M*N-1), M, N) + 1; 123 124 phi2D = phi1D(k); 125 xiplot = xinodo(k); 126 xplot = xinodo(k)*L; 127 etaplot = etanodo(k); 128 yplot = etanodo(k).*ys(xplot); 129 130 vx_f = zeros(M*N,4); 131 vy_f = zeros(M*N,4); 132 133 for j=1:N 134 for i=1:M 135 k = (j-1)*M+i; 136 index = [k-M-1 k-M k-M+1 k-1 k k+1 k+M-1 k+M k+M+1]; 137 starts_index = index(index > 0); 138 end_index = index(index <= M*N & index > 0); 139 140 eq_start = find(index == starts_index(1)); 141 eq_end = find(index == end_index(end)); 142 143 a=zeros(size(index)); 144 a(1,eq_start:eq_end) = phi1D(index(eq_start:eq_end))’; 145 146 phi_SW = a(1,1); 147 phi_S = a(1,2); 148 phi_SE = a(1,3); 149 phi_W = a(1,4); 150 phi_O = a(1,5); 151 phi_E = a(1,6);
88 Capítulo B. Códigos de MATLAB 134 + qscc*isteps(i,j)*bfimp_s*b_v; 135 qS2 = - Gamma/detJ(xis(k),etas(k))*(beta11(xis(k),etas(k))*dphidxi_s + beta12( xis(k),etas(k))*dphideta_s); 136 137 phi_re = upwind(vx_f(k,2), phi_O, phi_E)*stepe(i,j) + phi_O*istepe(i,j)*bfimp_e + phi_ecc*istepe(i,j)*b_v*bphiimp_e; 138 phi_rw = upwind(vx_f(k,4), phi_W, phi_O)*stepw(i,j) + phi_O*istepw(i,j)*bfimp_w + phi_wcc*istepw(i,j)*b_v*bphiimp_w; 139 phi_rn = upwind(vy_f(k,3), phi_O, phi_N)*stepn(i,j) + phi_O*istepn(i,j)*bfimp_n + phi_ncc*istepn(i,j)*b_v*bphiimp_n; 140 phi_rs = upwind(vy_f(k,1), phi_S, phi_O)*steps(i,j) + phi_O*isteps(i,j)*bfimp_s + phi_scc*isteps(i,j)*b_v*bphiimp_s; 141 142 eq = (rho*U1e*Cv*phi_re + qEcc*beta11(xie(k),etae(k)) + qE2*beta21(xie(k),etae(k )))*deltaeta ... 143 - (rho*U1w*Cv*phi_rw + qWcc*beta11(xiw(k),etaw(k)) + qW2*beta21(xiw(k),etaw(k )))*deltaeta ... 144 + (rho*U2n*Cv*phi_rn + qNcc*beta22(xin(k),etan(k)) + qN2 + qncc*b_v*istepn(i, j)*bfimp_n)*deltaxi ... 145 - (rho*U2s*Cv*phi_rs + qScc*beta22(xis(k),etas(k)) + qS2*beta12(xis(k),etas(k )))*deltaxi; 146 147 index = [k-M-1 k-M k-M+1 zeros(1,k-1-(k-M+1)-1) k-1 k k+1 zeros(1,k-1-(k-M+1)-1) k+M-1 k+M k+M+1]; 148 starts_index = index(index > 0 & index <= M*N); 149 starts = starts_index(1); ends = starts_index(end); 150 151 eq_start = find(index == starts_index(1)); 152 eq_end = find(index == starts_index(end)); 153 154 A(k, starts:ends) = eq(eq_start:eq_end); 155 b(k) = eq(end); 156 end 157 end 158 159 phi1D = A\b’; 160 161 k=reshape((0:M*N-1), M, N) + 1; 162 163 phi2D = phi1D(k); 164 xiplot = xinodo(k); 165 xplot = xinodo(k)*L; 166 etaplot = etanodo(k); 167 yplot = etanodo(k).*ys(xplot); 168 169 %% Representaciones 170 t_solucion=toc; 171 horas = floor(t_solucion/3600); 172 minuts = floor(mod(t_solucion,3600)/60); 173 seg = mod(t_solucion,60); 174 fprintf("Tiempo en calcular al solución: %d h %d min %4.1f s\n",horas, minuts, seg); 175 176 minphi = min(phi1D)-273.15; 177 maxphi = max(phi1D)-273.15; 178 % Curvas de nivel 179 figure(2)
B.2 Mallado no estructurado 89 180 contour(xplot,yplot,phi2D-273.15,minphi:(maxphi-minphi)/11:maxphi,’Linewidth’,2,’ ShowText’,’on’); 181 colormap(turbo) 182 xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 183 ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 184 zlabel(’$T\,(^{\circ}C)$’,’Interpreter’,’latex’) 185 title(’Curvas de nivel del Campo de Temperaturas’,’Interpreter’,’latex’) 186 187 figure(1) 188 h=surf(xplot,yplot,phi2D-273.15); 189 colormap(turbo) 190 cb=colorbar; 191 colorbar(’Ticks’,round(linspace(minphi+0.1,maxphi-0.1,10),1)) 192 view([0 90]) 193 axis equal 194 axis([0 L 0 H]) 195 shading interp 196 % set(h,’EdgeColor’,’k’) 197 grid off 198 xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 199 ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 200 zlabel(’$T\,^{\circ} C$’,’Interpreter’,’latex’) 201 titulo = sprintf(’Campo de Temperaturas para $v_1=%g \\, \\, m/s$’, v1); 202 title(titulo,’Interpreter’,’latex’) 203 204 function [phi_r] = upwind(v, phi_1, phi_2) 205 phi_r = phi_1*(sign(v) + 1)/2 + phi_2*(1 - sign(v))/2; 206 end B.2 Mallado no estructurado B.2.1 Problema difusivo con mallado cuadrangular 1clear;close all;clc; 2 3phi_scc = 40 + 273.15; bphiimp_s = 0; 4phi_ecc = 15 + 273.15; bphiimp_e = 0; 5phi_ncc = 40 + 273.15; bphiimp_n = 0; 6phi_wcc = 30 + 273.15; bphiimp_w = 1; 7 8% variable bfimp_x 1 si se impone el flujo en la cara ’x’, 0 si se impone 9% temperatura 10 qwcc = 0; bfimp_w = not(bphiimp_w); 11 qncc = 15; bfimp_n = not(bphiimp_n); 12 qecc = 0; bfimp_e = not(bphiimp_e); 13 qscc = 0; bfimp_s = not(bphiimp_s); 14 15 %% Mallado 16 17 H=2; 18 L=2; 19 Nx=30; % Número de vértices en la dirección x 20 Ny=30; % Número de vértices en la dirección y 21 M=Nx-1; % Número de celdas en x
90 Capítulo B. Códigos de MATLAB 22 N=Ny-1; % Número de celdas en y 23 24 x_ini=linspace(0,L,Nx); 25 eta=linspace(0,1,Ny); 26 27 %% Definición de la geometría 28 % Límite superior y = s2(x) 29 30 % Elipse 31 a = L*1.2; b = H; 32 S2 = @(x) b*sqrt(1-x.^2/a^2); 33 s2=S2(x_ini); 34 35 xvert = zeros(Nx,Ny); 36 yvert = zeros(Nx,Ny); 37 for i=1:Nx 38 xvert(i,1:Ny) = x_ini(i); 39 yvert(i,1:Ny) = eta*s2(i); 40 end 41 42 %% Información a nivel Global 43 % Número total de nodos 44 n_c = M*N; 45 46 % Número total de vértices 47 48 n_v = Nx*Ny; 49 50 % Coordenadas de los vértices 51 52 iv = 0; 53 xv=zeros(1,n_v); 54 yv=zeros(1,n_v); 55 for j=1:Ny 56 for i=1:Nx 57 iv = iv + 1; 58 xv(iv) = xvert(i,j); 59 yv(iv) = yvert(i,j); 60 end 61 end 62 63 % Coordenadas de los nodos 64 65 ic = 0; 66 xc=zeros(1,n_c); 67 yc=zeros(1,n_c); 68 for j=1:N 69 for i=1:M 70 ic = ic + 1; 71 iv = ic + (j-1); % Vértice asociado a la esquina inferior izquierda de la celda 72 xc(ic) = (xv(iv) + xv(iv+1))/2; 73 yc(ic) = (yv(iv) + yv(iv+Nx) + yv(iv+1) + yv(iv+Nx+1))/4; 74 end 75 end 76 77 %% Conectividad 78 % Conectividad Celda-A-Cara
B.2 Mallado no estructurado 91 79 80 ic = 0; 81 Lcf = zeros(n_c,4); 82 for j=1:N 83 for i=1:M 84 ic = ic + 1; 85 Lcf(ic,1) = ic; 86 Lcf(ic,2) = M*Ny + N*i + j; 87 Lcf(ic,3) = ic + M; 88 Lcf(ic,4) = M*Ny + N*(i-1) + j; 89 end 90 end 91 92 n_f = Lcf(ic,2); 93 94 % Conectividad Cara-A-Celda 95 Lfc = zeros(n_f,2); 96 for iface=1:n_f 97 c1 = 0; c2 = 0; 98 for ic=1:n_c 99 for k=1:4 100 if Lcf(ic,k) == iface 101 if c1 == 0 102 c1 = ic; 103 else 104 c2 = ic; 105 end 106 end 107 end 108 end 109 Lfc(iface,1) = c1; 110 Lfc(iface,2) = c2; 111 if Lfc(iface,1) == 0 112 Lfc(iface,1) = Lfc(iface,2); 113 elseif Lfc(iface,2) == 0 114 Lfc(iface,2) = Lfc(iface,1); 115 end 116 end 117 118 119 % Conectividad Celda-A-Vértice 120 ic = 0; 121 Lcv = zeros(n_c,4); 122 for j=1:N 123 for i=1:M 124 ic = ic + 1; 125 Lcv(ic,1) = ic + (j-1); 126 Lcv(ic,2) = ic + j; % Es uno más que el primer vértice de la celda 127 Lcv(ic,3) = Nx*j + (i+1); 128 Lcv(ic,4) = Nx*j + i; 129 end 130 end 131 132 % Conectividad Cara-A-Vértice 133 Lfv = zeros(n_f,2); 134 for iface=1:M*N 135 ic = Lfc(iface,2);
92 Capítulo B. Códigos de MATLAB 136 Lfv(iface,1) = Lcv(ic,1); 137 Lfv(iface,2) = Lcv(ic,2); 138 end 139 for iface=(M*N+1):M*Ny 140 ic = Lfc(iface,2); 141 Lfv(iface,1) = Lcv(ic,4); 142 Lfv(iface,2) = Lcv(ic,3); 143 end 144 for iface=(M*Ny+1):(M*Ny+N*M) 145 ic = Lfc(iface,2); 146 Lfv(iface,1) = Lcv(ic,1); 147 Lfv(iface,2) = Lcv(ic,4); 148 end 149 for iface=(M*Ny+N*M+1):n_f 150 ic = Lfc(iface,2); 151 Lfv(iface,1) = Lcv(ic,2); 152 Lfv(iface,2) = Lcv(ic,3); 153 end 154 155 % Centros de las caras 156 xf=zeros(n_f,1); 157 yf=zeros(n_f,1); 158 for iface=1:n_f 159 xf(iface) = (xv(Lfv(iface,1))+xv(Lfv(iface,2)))/2; 160 yf(iface) = (yv(Lfv(iface,1))+yv(Lfv(iface,2)))/2; 161 end 162 163 164 %% Indexación de Caras Frontera 165 n_f_b = 2*(M+N); % Total en Contorno 166 bface = zeros(1,n_f); % 1 SI cara en Contorno, 0 SI cara Interior 167 f_to_bf = zeros(1,n_f); 168 bf_to_f = zeros(1,n_f_b); 169 bvert = zeros(1,n_v); 170 iface_b = 0; 171 bphiimp_bf = zeros(1,n_f_b); 172 phi_bf = zeros(1,n_f_b); 173 q_bf = zeros(1,n_f_b); 174 bfimp_f = zeros(1,n_f); 175 % Cara Sur 176 for i=1:M 177 ic = i; 178 iface = ic; 179 iface_b = iface_b + 1; 180 bface(iface) = 1; 181 f_to_bf(iface) = iface_b; 182 bf_to_f(iface_b) = iface; 183 bvert(Lfc(iface,1)) = 1; 184 phi_bf(iface_b) = phi_scc; 185 q_bf(iface_b) = qscc; 186 bphiimp_bf(iface_b) = bphiimp_s; 187 bfimp_f(iface) = not(bphiimp_bf(iface_b)); 188 end 189 % Cara Este 190 for j=1:N 191 ic = i + M*(j-1); 192 iface = Lcf(ic,2);
B.2 Mallado no estructurado 93 193 iface_b = iface_b + 1; 194 bface(iface) = 1; 195 f_to_bf(iface) = iface_b; 196 bf_to_f(iface_b) = iface; 197 bvert(Lfv(iface,1)) = 1; 198 phi_bf(iface_b) = phi_ecc; 199 q_bf(iface_b) = qecc; 200 bphiimp_bf(iface_b) = bphiimp_e; 201 bfimp_f(iface) = not(bphiimp_bf(iface_b)); 202 end 203 % Cara Norte 204 for i=M:-1:1 205 ic = M*(j-1) + i; 206 iface = Lcf(ic,3); 207 iface_b = iface_b + 1; 208 bface(iface) = 1; 209 f_to_bf(iface) = iface_b; 210 bf_to_f(iface_b) = iface; 211 bvert(Lfv(iface,2)) = 1; 212 phi_bf(iface_b) = phi_ncc; 213 q_bf(iface_b) = qncc; 214 bphiimp_bf(iface_b) = bphiimp_n; 215 bfimp_f(iface) = not(bphiimp_bf(iface_b)); 216 end 217 % Cara Oeste 218 for j=N:-1:1 219 ic = M*(j-1) + 1; 220 iface = Lcf(ic,4); 221 iface_b = iface_b + 1; 222 bface(iface) = 1; 223 f_to_bf(iface) = iface_b; 224 bf_to_f(iface_b) = iface; 225 bvert(Lfv(iface,2)) = 1; 226 phi_bf(iface_b) = phi_wcc; 227 q_bf(iface_b) = qwcc; 228 bphiimp_bf(iface_b) = bphiimp_w; 229 bfimp_f(iface) = not(bphiimp_bf(iface_b)); 230 end 231 232 %% CALCULO DE AREAS Y NORMALES Y TANGENTES UNITARIAS 233 234 lfcx = zeros(n_c,4); 235 lfcy = zeros(n_c,4); 236 area = zeros(1,n_f); 237 tx=zeros(n_c,4); 238 ty=zeros(n_c,4); 239 nx=zeros(n_c,4); 240 ny=zeros(n_c,4); 241 for ic=1:n_c 242 for k=1:4 243 iface = Lcf(ic,k); 244 lfcx(ic,k) = xf(iface) - xc(ic); 245 lfcy(ic,k) = yf(iface) - yc(ic); 246 v1 = Lfv(iface,1); 247 v2 = Lfv(iface,2); 248 dx1 = xv(v2) - xv(v1); 249 dy1 = yv(v2) - yv(v1);
94 Capítulo B. Códigos de MATLAB 250 area(iface) = sqrt(dx1^2 + dy1^2); 251 tx(ic,k) = dx1/area(iface); 252 ty(ic,k) = dy1/area(iface); 253 nx(ic,k) = ty(ic,k); 254 ny(ic,k) = -tx(ic,k); 255 scalar_nl = nx(ic,k)*lfcx(ic,k) + ny(ic,k)*lfcy(ic,k); 256 if scalar_nl < 0 257 tx(ic,k) = -tx(ic,k); 258 ty(ic,k) = -ty(ic,k); 259 nx(ic,k) = -nx(ic,k); 260 ny(ic,k) = -ny(ic,k); 261 end 262 end 263 end 264 265 %% Cálculo de DELTA 266 267 delta = zeros(1,n_f); 268 for iface=1:n_f 269 c1 = Lfc(iface,1); 270 c2 = Lfc(iface,2); 271 for ic=1:n_c 272 for k=1:4 273 if Lcf(ic,k) == iface 274 nx1 = nx(ic,k); 275 ny1 = ny(ic,k); 276 end 277 end 278 end 279 280 if c1 == c2 % Cara frontera 281 dx1 = xf(iface) - xc(c1); 282 dy1 = yf(iface) - yc(c1); 283 delta(iface) = abs(dx1*nx1 + dy1*ny1); 284 else 285 dx1 = xc(c2) - xc(c1); 286 dy1 = yc(c2) - yc(c1); 287 delta(iface) = abs(dx1*nx1 + dy1*ny1); 288 end 289 end 290 291 %% Cálculo del VOLUMEN 292 vol(:) = zeros(1,n_c); 293 for ic=1:n_c 294 v1 = Lcv(ic,1); 295 v2 = Lcv(ic,2); 296 v3 = Lcv(ic,3); 297 v4 = Lcv(ic,4); 298 dx1 = xv(v2) - xv(v1); 299 dy1 = yv(v2) - yv(v1); 300 vector1 = [dx1; dy1]; 301 dx2 = xv(v4) - xv(v1); 302 dy2 = yv(v4) - yv(v1); 303 vector2 = [dx2; dy2]; 304 semiVol1 = 0.5*abs(dx1*dy2 - dy1*dx2); % semiArea = 1/2*|t1 x t2| 305 306 dx3 = xv(v4) - xv(v3);
B.2 Mallado no estructurado 95 307 dy3 = yv(v4) - yv(v3); 308 vector3 = [dx3; dy3]; 309 dx4 = xv(v2) - xv(v3); 310 dy4 = yv(v2) - yv(v3); 311 vector4 = [dx4; dy4]; 312 semiVol2 = 0.5*abs(dx3*dy4 - dy3*dx4); % semiArea = 1/2*|t1 x t2| 313 314 vol(ic) = semiVol1 + semiVol2; 315 end 316 317 318 %% Funciones de interpolación Nodo-A-Cara 319 320 wfun = zeros(1,n_f); 321 gamma = 0.6*ones(1,n_c); 322 gammaf = zeros(1,n_f); 323 for iface=1:n_f 324 c1 = Lfc(iface); 325 c2 = Lfc(iface); 326 dx1 = xc(c1) - xf(iface); 327 dy1 = yc(c1) - yf(iface); 328 d1=sqrt(dx1^2 + dy1^2); 329 dx2 = xc(c2) - xf(iface); 330 dy2 = yc(c2) - yf(iface); 331 d2=sqrt(dx2^2 + dy2^2); 332 wfun(iface) = d2/(d1+d2); 333 gammaf(iface) = wfun(iface)*gamma(c1) + (1-wfun(iface))*gamma(c2); 334 end 335 336 337 %% Funciones interpolación Nodo-A-Vértice 338 339 dtotal = zeros(n_v,4); 340 di=zeros(n_c,4); 341 for ic=1:n_c 342 for iv=1:4 343 v1 = Lcv(ic,iv); 344 di(ic,iv) = sqrt((xv(v1) - xc(ic))^2 + (yv(v1) - yc(ic))^2); 345 dtotal(v1) = dtotal(v1) + 1/di(ic,iv); 346 end 347 end 348 cwfun = zeros(n_c,4); 349 for ic=1:n_c 350 for iv=1:4 351 v1 = Lcv(ic,iv); 352 cwfun(ic,iv) = (1/di(ic,iv))/(dtotal(v1)); 353 end 354 end 355 356 357 %% Cálculo de la solución 358 359 err = 10^3; % Error inicial para entrar en el bucle 360 TOL = 10^-3; % Tolerancia 361 phi = zeros(n_c,1); 362 363 flag = 0;
96 Capítulo B. Códigos de MATLAB 364 while max(err,[],’all’) > TOL 365 366 skew = zeros(1,n_c); 367 b=zeros(1,n_c); 368 sc=zeros(n_c); 369 370 if flag==0 371 A=zeros(n_c, n_c); 372 for ic=1:n_c 373 for k=1:4 374 iface = Lcf(ic,k); 375 if bface(iface) == 0 % Cara interior 376 icvec = Lfc(iface,:); 377 A(ic,icvec(icvec ~= ic)) = -gammaf(iface)*area(iface)/delta(iface); 378 end 379 end 380 A(ic,ic) = sum(gammaf(Lcf(ic,:)).*area(Lcf(ic,:))./delta(Lcf(ic,:)).*not( bfimp_f(Lcf(ic,:)))); 381 end 382 end 383 384 % Cálculo de phi en los vértices 385 if flag == 1 386 phiv = zeros(1,n_v); 387 for ic=1:n_c 388 for k=1:4 389 iv = Lcv(ic,k); 390 % Cara conectado al vértice 391 for iface=1:4 392 if Lfv(Lcf(ic,k), 1) == iv 393 ifacecond = iface; 394 elseif Lfv(Lcf(ic,k), 2) == iv 395 ifacecond = iface; 396 end 397 end 398 if bvert == 0 399 phiv(iv) = phiv(iv) + phi_ini(ic)*cwfun(ic,k); 400 else 401 xif = xf(iface)/L; 402 etaf = yf(iface)/S2(xif*L); 403 if bphiimp_bf(ifacecond) == 1 404 phiv(iv) = phi_bf(ifacecond); 405 else 406 phiv(iv) = phiv(iv) + phi_ini(ic)*cwfun(ic,k); 407 end 408 end 409 end 410 end 411 412 % Cálculo término skew 413 414 sign =zeros(n_c,4); 415 for ic=1:n_c 416 for k=1:4 417 iface = Lcf(ic,k); 418 vb = Lfv(iface,1); 419 va = Lfv(iface,2);
B.2 Mallado no estructurado 97 420 ic1 = Lfc(iface,1); 421 ic2 = Lfc(iface,2); 422 Lf(1) = xc(ic2) - xc(ic1); 423 Lf(2) = yc(ic2) - yc(ic1); 424 tf(1) = xv(va) - xv(vb); 425 tf(2) = yv(va) - yv(vb); 426 utf = tf(:)/area(iface); 427 utf_dot_Lf = utf(1)*Lf(1) + utf(2)*Lf(2); 428 scalar_nl = nx(ic,k)*Lf(1) + ny(ic,k)*Lf(2); 429 if scalar_nl > 0 430 sign(ic,k) = 1; 431 else 432 sign(ic,k) = -1; 433 end 434 skew(ic) = skew(ic) - gammaf(iface)*utf_dot_Lf*sign(ic,k)*(phiv(va) - phiv(vb))/delta(iface); 435 end 436 b(ic) = skew(ic) - sc(ic)*vol(ic); 437 end 438 end 439 440 b(Lfc(bf_to_f(1:n_f_b),1)) = b(Lfc(bf_to_f(1:n_f_b),1)) + area(bf_to_f).*q_bf.*not( bphiimp_bf) + gammaf(bf_to_f).*area(bf_to_f).*phi_bf./delta(bf_to_f).* bphiimp_bf; 441 442 % Solución 443 444 phi = A\b’; 445 if flag == 1 446 err = abs(phi - phi_ini); 447 end 448 phi_ini = phi; 449 flag = 1; 450 451 end 452 453 454 %% REPRESENTACION GRAFICA 455 if M*N<20 456 figure(2) 457 axis off 458 hold on 459 vlabels = cell(1,n_v); 460 flabels = cell(1,n_f); 461 clabels = cell(1,n_c); 462 463 % Draw faces 464 for i=1:n_f 465 plot([xvert(Lfv(i,1)) xvert(Lfv(i,2))], [yvert(Lfv(i,1)) yvert(Lfv(i,2))], ’ Color’,’#126D70’,’LineWidth’, 1.2); 466 plot(xf(i), yf(i), ’x’,’Color’,’#212F3D’,’MarkerSize’,9, ’LineWidth’,1.1); 467 flabels(i) = {int2str(i)}; 468 end 469 text(xf,yf,flabels,’VerticalAlignment’,’bottom’,’HorizontalAlignment’,’right’,’ Color’,’#212F3D’,’Interpreter’,’latex’,’FontSize’,12) 470 471 % Draw vertices
104 Capítulo B. Códigos de MATLAB 238 v1 = Lcv(ic,iv); 239 d1(ic,iv) = sqrt((xv(v1) - xc(ic))^2 + (yv(v1) - yc(ic))^2); 240 % cwfun(ic,iv) = 1/d1; 241 sumwt(v1) = sumwt(v1) + 1/d1(ic,iv); 242 end 243 end 244 245 for ic=1:n_c 246 for iv=1:3 247 v1 = Lcv(ic,iv); 248 cwfun(ic,iv) = (1/d1(ic,iv))/sumwt(v1); 249 end 250 end 251 252 %% Cálculo de la solución 253 254 err = 10^3; % Error inicial para entrar en el bucle 255 TOL = 10^-3; % Tolerancia 256 phi = zeros(n_c,1); 257 258 flag = 0; 259 while max(err,[],’all’) > TOL 260 261 skew = zeros(1,n_c); 262 b=zeros(1,n_c); 263 sc=zeros(n_c); 264 265 if flag==0 266 A=zeros(n_c, n_c); 267 for ic=1:n_c 268 for k=1:3 269 iface = Lcf(ic,k); 270 if bface(iface) == 0 % Cara interior 271 icvec = Lfc(iface,:); 272 A(ic,icvec(icvec ~= ic)) = -gammaf(iface)*area(iface)/delta(iface); 273 end 274 end 275 A(ic,ic) = sum(gammaf(Lcf(ic,:)).*area(Lcf(ic,:))./delta(Lcf(ic,:))); 276 end 277 end 278 279 % Cálculo de phi en los vértices 280 if flag == 1 281 phiv = zeros(1,n_v); 282 for ic=1:n_c 283 for k=1:3 284 iv = Lcv(ic,k); 285 if bvert(iv) == 0 286 phiv(iv) = phiv(iv) + phi_ini(ic)*cwfun(ic,k); 287 else 288 phiv(iv) = phi_bv(iv); 289 end 290 end 291 end 292 293 % Cálculo término skew 294
B.2 Mallado no estructurado 105 295 sign =zeros(n_c,4); 296 for ic=1:n_c 297 for k=1:3 298 iface = Lcf(ic,k); 299 vb = Lfv(iface,1); 300 va = Lfv(iface,2); 301 ic1 = Lfc(iface,1); 302 ic2 = Lfc(iface,2); 303 Lf(1) = xc(ic2) - xc(ic1); 304 Lf(2) = yc(ic2) - yc(ic1); 305 tf(1) = xv(va) - xv(vb); 306 tf(2) = yv(va) - yv(vb); 307 utf = tf(:)/area(iface); 308 utf_dot_Lf = utf(1)*Lf(1) + utf(2)*Lf(2); 309 scalar_nl = nx(ic,k)*Lf(1) + ny(ic,k)*Lf(2); 310 if scalar_nl > 0 311 sign(ic,k) = 1; 312 else 313 sign(ic,k) = -1; 314 end 315 skew(ic) = skew(ic) - gammaf(iface)*utf_dot_Lf*sign(ic,k)*(phiv(va) - phiv(vb))/delta(iface); 316 end 317 b(ic) = skew(ic) - sc(ic)*vol(ic); 318 end 319 end 320 321 for ibf=1:n_f_b 322 iface = bf_to_f(ibf); 323 ic1 = Lfc(iface,1); 324 b(ic1) = b(ic1) + gammaf(iface).*area(iface).*phi_bf(ibf)./delta(iface); 325 end 326 327 % Solución 328 329 phi = A\b’; 330 if flag == 1 331 err = abs(phi - phi_ini); 332 end 333 phi_ini = phi; 334 flag = 1; 335 336 end 337 338 %% Cálculo de la velocidad 339 340 phiv = zeros(1,n_v); 341 for ic=1:n_c 342 for k=1:3 343 iv = Lcv(ic,k); 344 if bvert(iv) == 0 345 phiv(iv) = phiv(iv) + phi(ic)*cwfun(ic,k); 346 else 347 phiv(iv) = phi_bv(iv); 348 end 349 end 350 end
106 Capítulo B. Códigos de MATLAB 351 352 gphi_t = zeros(n_f,1); 353 gphi_n = zeros(n_f,1); 354 vx_f = zeros(n_f,1); 355 vy_f = zeros(n_f,1); 356 signo = zeros(n_c,3); 357 vn=zeros(n_c,3); 358 359 for ic=1:n_c 360 for k=1:3 361 iface = Lcf(ic,k); 362 va = Lfv(iface,1); 363 vb = Lfv(iface,2); 364 ic12 = Lfc(iface,:); 365 ic2 = ic12(ic12~=ic); 366 signo(ic,k) = ((xv(va)-xv(vb))*tx(ic,k) + (yv(va)-yv(vb))*ty(ic,k))/(area(iface) ); 367 if not(bface(iface)) 368 tdotl = signo(ic,k)*((xv(va)-xv(vb))*(xc(ic2)-xc(ic)) + (yv(va)-yv(vb))*(yc( ic2)-yc(ic))); 369 gphi_t(iface) = signo(ic,k)*(phiv(va)-phiv(vb))/area(iface); 370 gphi_n(iface) = (phi(ic2)-phi(ic))/delta(iface) - (signo(ic,k)*(phiv(va)-phiv (vb)))/area(iface)*tdotl/(area(iface)*delta(iface)); 371 else 372 tdotl = signo(ic,k)*((xv(va)-xv(vb))*(xf(iface)-xc(ic)) + (yv(va)-yv(vb))*(yf (iface)-yc(ic))); 373 gphi_t(iface) = signo(ic,k)*(phiv(va)-phiv(vb))/area(iface); 374 gphi_n(iface) = (phi_bf(f_to_bf(iface))-phi(ic))/delta(iface) - (signo(ic,k) *(phiv(va)-phiv(vb)))*tdotl/(area(iface)^2*delta(iface)); 375 end 376 vx_f(iface) = gphi_t(iface)*nx(ic,k) - gphi_n(iface)*tx(ic,k); 377 vy_f(iface) = gphi_t(iface)*ny(ic,k) - gphi_n(iface)*ty(ic,k); 378 vn(ic,k) = gphi_t(iface); 379 end 380 end 381 382 vx_c = zeros(n_c,1); 383 vy_c = zeros(n_c,1); 384 for ic=1:n_c 385 d = 1./sqrt((xc(ic) - xf(Lcf(ic,:))).^2 + (yc(ic) - yf(Lcf(ic,:))).^2); 386 vx_c(ic) = sum(vx_f(Lcf(ic,:)).*d./sum(d)); 387 vy_c(ic) = sum(vy_f(Lcf(ic,:)).*d./sum(d)); 388 end 389 390 391 %% Representación 392 393 % Representación del Mallado 394 395 if n_c < 30 396 figure(4) 397 hold on 398 vlabels = {}; 399 400 % Draw faces 401 for i=1:n_f 402 plot([xface(i,1) xface(i,2)], [yface(i,1) yface(i,2)], ’Color’,’#212F3D’);
B.2 Mallado no estructurado 107 403 flabels(i) = {int2str(i)}; 404 end 405 text(xf,yf,flabels,’VerticalAlignment’,’bottom’,’HorizontalAlignment’,’right’,’ Color’,’#212F3D’) 406 407 % Draw vertices 408 for i=1:length(xv) 409 scatter(xv(i),yv(i),’filled’,’MarkerFaceColor’,’#FF5733’); 410 vlabels(i) = {int2str(i)}; 411 end 412 text(xv,yv,vlabels,’VerticalAlignment’,’bottom’,’HorizontalAlignment’,’right’,’ Color’,’#FF5733’) 413 414 % Draw incenters 415 for i=1:length(xc) 416 plot(xc(i),yc(i),’o’,’Color’,’#5903A1’) 417 clabels(i) = {int2str(i)}; 418 end 419 text(xc,yc,clabels,’VerticalAlignment’,’bottom’,’HorizontalAlignment’,’right’,’ Color’,’#5903A1’) 420 end 421 %% Representación de la Solución 422 423 % Velocidades 424 425 Mv = 8; 426 Nv = 6; 427 xrepr = zeros(Mv*Nv,1); 428 yrepr = zeros(Mv*Nv,1); 429 ic_list = zeros(Mv*Nv,1); 430 for j=1:Nv 431 for i=1:Mv 432 k = (j-1)*Mv+i; 433 xrepr(k) = L/(Mv+1)/2 + (i-1)*(L/(Mv+1)); 434 yrepr(k) = S2(xrepr(k))/(Nv)/2 + (j-1)*S2(xrepr(k))/Nv; 435 end 436 end 437 438 for j=1:Nv 439 for i=1:Mv 440 k = (j-1)*Mv+i; 441 [~, ic_list(k)] = min(sqrt((xc-xrepr(k)).^2 + (yc-yrepr(k)).^2)); 442 end 443 end 444 % figure(5) 445 % hold on 446 % axis([-0.1*L L*1.1 -0.1*H H*1.1]); 447 % axis(’equal’) 448 % tckness = 1.5; 449 % linecolor = ’#CE9FFF’; % ’#5bd966’ 450 % plot([0 L], [0 0], ’Color’,linecolor,’LineWidth’,tckness); 451 % plot([0 0],[0 H1], ’Color’,linecolor,’LineWidth’,tckness); 452 % plot([0:0.01:L], S2([0:0.01:L]), ’Color’,linecolor,’LineWidth’,tckness); 453 % plot([L L], [0 H2], ’Color’,linecolor,’LineWidth’,tckness); 454 % title(’Campo de velocidades’, ’Interpreter’,’latex’) 455 % 456 % % Dibujar el campo de vectores
108 Capítulo B. Códigos de MATLAB 457 % quiver(xrepr, yrepr, vx_c(ic_list), vy_c(ic_list), ’Color’,’#388ADC’,’linewidth’,1.2); 458 % 459 % 460 % 461 % % Ajustar los ejes para que se vean correctamente 462 % axis equal; 463 % title(’Campo de velocidades’); 464 % xlabel(’X’); 465 % ylabel(’Y’); 466 % hold off; 467 468 % Phi 469 470 % figure(2) 471 % tri = delaunay(xc, yc); 472 % h=trisurf(tri, xc, yc, phi); 473 % colormap(turbo) 474 % colorbar 475 % shading interp 476 % % set(h,’EdgeColor’,’k’) 477 % xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 478 % ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 479 % zlabel(’$T\,(^{\circ}C)$’,’Interpreter’,’latex’) 480 % title(’Campo de Temperaturas’,’Interpreter’,’latex’) 481 482 483 %% Aproximación 484 M=50; 485 N=50; 486 Nx = M+1; % Número de vértices en x 487 Ny = N+1; % Número de vértices en y 488 xi=zeros(M,N); eta = zeros(M,N); 489 yvert = zeros(M,N); 490 491 % Distribución xi y eta. 492 for j=1:Ny 493 xi(1:Nx, j) = linspace(0,1,Nx); 494 end 495 for i=1:Nx 496 eta(i,1:Ny) = linspace(0,1,Ny); 497 yvert(i,1:Ny) = linspace(0,S2(xi(i)*L),Ny); 498 end 499 xvert = xi*L; 500 501 n_ce = M*N; % Número total de caras 502 503 xi1 = reshape(xi(1:M,1:N), 1, n_ce); eta1 = reshape(eta(1:M,1:N), 1, n_ce); 504 xi2 = reshape(xi(1:M,2:N+1), 1, n_ce); eta2 = reshape(eta(1:M,2:N+1), 1, n_ce); 505 xi3 = reshape(xi(2:M+1,2:N+1), 1, n_ce); eta3 = reshape(eta(2:M+1,2:N+1), 1, n_ce); 506 xi4 = reshape(xi(2:M+1,1:N), 1, n_ce); eta4 = reshape(eta(2:M+1,1:N), 1, n_ce); 507 % Nodos 508 xinodo = (xi1+xi4)/2; 509 etanodo = (eta1 + eta2)/2; 510 511 k=reshape((0:M*N-1), M, N) + 1; 512 xplot = xinodo(k)*L;
B.2 Mallado no estructurado 109 513 yplot = etanodo(k).*S2(xplot); 514 515 F = scatteredInterpolant(xc,yc,phi); 516 phi2 = F(xinodo*L, etanodo.*S2(xinodo*L)); 517 518 519 for j=1:50 520 for i=1:50 521 k = (j-1)*50 + i; 522 phi2D(i,j) = phi2(k); 523 end 524 end 525 minphi = min(phi); 526 maxphi = max(phi); 527 % 528 % figure(1) 529 % h=surf(xplot,yplot,phi2D); 530 % colormap(turbo) 531 % cb = colorbar; 532 % nticks = 11; 533 % % set(cb,’YTick’,linspace(minphi,maxphi,nticks)); 534 % % numberOfDecimals = 1; 535 % 536 % % Crear las etiquetas de las marcas con el número de decimales especificado 537 % % tickLabels = arrayfun(@(x) sprintf(’%4.1f’, round(x,1)), linspace(minphi, maxphi, nticks), ’UniformOutput’, false); 538 % 539 % % Establecer las etiquetas en la barra de colores 540 % % set(cb, ’YTickLabel’, tickLabels); 541 % shading interp 542 % % set(h,’EdgeColor’,’k’) 543 % xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 544 % ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 545 % zlabel(’$T\,^{\circ} C$’,’Interpreter’,’latex’) 546 % titulo = sprintf(’Campo de Temperaturas’); 547 % title(titulo,’Interpreter’,’latex’) B.2.3 Problema convectivo-difusivo con una malla externa 1%% STL Read MATLAB Script 2clear;close all;clc 3 4config; 5 6%% ------------------------ Mallado ------------------------ %% 7 8% TR = stlread("exponencial.stl"); 9 10 % Número total de nodos y Coordenadas de los nodos 11 % Incenter 12 incenter = TR.incenter; 13 xc = incenter(:,1); 14 yc = incenter(:,2); 15 n_c = length(incenter); 16
110 Capítulo B. Códigos de MATLAB 17 % Número total de vértices y Coordenadas de los vértices 18 vertices = TR.Points; 19 xv = vertices(:,1); 20 yv = vertices(:,2); 21 n_v = length(vertices); 22 L=max(xv); 23 idH = find(xv == 0); H = max(yv(idH)); 24 25 % Número total de caras y Centros de las caras 26 % Edges (vert1 --> vert2) 27 edges = TR.edges; 28 v1face = edges(:,1); 29 v2face = edges(:,2); 30 n_f = length(edges); 31 % Faces coordinates 32 xface = [xv(v1face(1:n_f)) xv(v2face(1:n_f))]; 33 yface = [yv(v1face(1:n_f)) yv(v2face(1:n_f))]; 34 35 xf = (xface(:,1) + xface(:,2))/2; 36 yf = (yface(:,1) + yface(:,2))/2; 37 38 %% Conectividades 39 % Conectividad Celda-A-Vértice 40 Lcv = TR.ConnectivityList; 41 42 % Conectividad Cara-A-Vértice 43 Lfv = [v1face v2face]; 44 45 % Conectividad Cara-A-Celda 46 for iface=1:n_f 47 face_con = edgeAttachments(TR,v1face(iface),v2face(iface)); 48 Lfc(iface,1:2) = face_con{:}; 49 end 50 51 % Conectividad Celda-A-Cara 52 for ic=1:n_c 53 k=1; 54 for iface=1:n_f 55 if ismember(ic, Lfc(iface,:)) 56 Lcf(ic,k) = iface; 57 k = k+1; 58 end 59 end 60 end 61 62 %% Indexación de Caras Frontera 63 n_f_b = nnz(Lfc(:,1)==Lfc(:,2)); % Total en Contorno 64 bface = zeros(1,n_f); % 1 SI cara en Contorno, 0 SI cara Interior 65 f_to_bf = zeros(1,n_f); 66 bf_to_f = zeros(1,n_f_b); 67 bvert = zeros(1,n_v); 68 iface_b = 0; 69 bphiimp_bf = zeros(1,n_f_b); 70 phi_bf = zeros(1,n_f_b); 71 q_bf = zeros(1,n_f_b); 72 bfimp_f = zeros(1,n_f); 73
B.2 Mallado no estructurado 111 74 for iface=1:n_f 75 if Lfc(iface,1) == Lfc(iface,2) 76 iface_b = iface_b + 1; 77 bface(iface) = 1; 78 f_to_bf(iface) = iface_b; 79 bf_to_f(iface_b) = iface; 80 xif = xf(iface)/L; 81 etaf = yf(iface)/S2(xif*L); 82 if etaf < 0.01 % Cara sur 83 phi_bf(iface_b) = phi_scc; 84 q_bf(iface_b) = qscc; 85 bphiimp_bf(iface_b) = bphiimp_s; 86 bfimp_f(iface) = bfimp_s; 87 end 88 if xif > 0.99 % Cara este 89 q_bf(iface_b) = qecc; 90 phi_bf(iface_b) = phi_ecc; 91 bphiimp_bf(iface_b) = bphiimp_e; 92 bfimp_f(iface) = bfimp_e; 93 end 94 if xif < 0.01 % Cara oeste 95 phi_bf(iface_b) = phi_wcc; 96 q_bf(iface_b) = qwcc; 97 bphiimp_bf(iface_b) = bphiimp_w; 98 bfimp_f(iface) = bfimp_w; 99 end 100 if etaf > 0.99 % Cara norte 101 q_bf(iface_b) = qncc; 102 phi_bf(iface_b) = phi_ncc; 103 bphiimp_bf(iface_b) = bphiimp_n; 104 bfimp_f(iface) = bfimp_n; 105 end 106 107 end 108 end 109 110 %% CALCULO DE AREAS Y NORMALES Y TANGENTES UNITARIAS 111 112 lfcx = zeros(n_c,3); 113 lfcy = zeros(n_c,3); 114 area = zeros(1,n_f); 115 tx=zeros(n_c,3); 116 ty=zeros(n_c,3); 117 nx=zeros(n_c,3); 118 ny=zeros(n_c,3); 119 for ic=1:n_c 120 for k=1:3 121 iface = Lcf(ic,k); 122 lfcx(ic,k) = xf(iface) - xc(ic); 123 lfcy(ic,k) = yf(iface) - yc(ic); 124 v1 = Lfv(iface,1); 125 v2 = Lfv(iface,2); 126 dx1 = (xv(v2) - xv(v1)); 127 dy1 = (yv(v2) - yv(v1)); 128 area(iface) = sqrt(dx1^2 + dy1^2); 129 tx(ic,k) = dx1/area(iface); 130 ty(ic,k) = dy1/area(iface);
112 Capítulo B. Códigos de MATLAB 131 nx(ic,k) = ty(ic,k); 132 ny(ic,k) = -tx(ic,k); 133 scalar_nl = nx(ic,k)*lfcx(ic,k) + ny(ic,k)*lfcy(ic,k); 134 if scalar_nl < 0 135 tx(ic,k) = -tx(ic,k); 136 ty(ic,k) = -ty(ic,k); 137 nx(ic,k) = -nx(ic,k); 138 ny(ic,k) = -ny(ic,k); 139 end 140 end 141 end 142 143 %% Cálculo de DELTA 144 145 delta = zeros(1,n_f); 146 for iface=1:n_f 147 c1 = Lfc(iface,1); 148 c2 = Lfc(iface,2); 149 for ic=1:n_c 150 for k=1:3 151 if Lcf(ic,k) == iface 152 nx1 = nx(ic,k); 153 ny1 = ny(ic,k); 154 end 155 end 156 end 157 if c1 == c2 % Cara frontera 158 dx1 = xf(iface) - xc(c1); 159 dy1 = yf(iface) - yc(c1); 160 delta(iface) = abs(dx1*nx1 + dy1*ny1); 161 else 162 dx1 = xc(c2) - xc(c1); 163 dy1 = yc(c2) - yc(c1); 164 delta(iface) = abs(dx1*nx1 + dy1*ny1); 165 end 166 end 167 168 %% Cálculo del VOLUMEN 169 170 vol(:) = 0; 171 for ic=1:n_c 172 v1 = Lcv(ic,1); 173 v2 = Lcv(ic,2); 174 v3 = Lcv(ic,3); 175 dx1 = xv(v2) - xv(v1); 176 dy1 = yv(v2) - yv(v1); 177 dx2 = xv(v3) - xv(v1); 178 dy2 = yv(v3) - yv(v1); 179 vol(ic) = 0.5*abs(dx1*dy2 - dy1*dx2); % A = 1/2*|t1 x t2| 180 end 181 182 %% Funciones de interpolación Nodo-A-Cara 183 184 wfun = zeros(1,n_f); 185 gamma = Gamma*ones(1,n_c); 186 gammaf = zeros(1,n_f); 187 for iface=1:n_f
B.2 Mallado no estructurado 113 188 c1 = Lfc(iface,1); 189 c2 = Lfc(iface,2); 190 dx1 = xf(iface) - xc(c1); 191 dy1 = yf(iface) - yc(c1); 192 d1=sqrt(dx1^2+dy1^2); 193 dx2 = xf(iface) - xc(c2); 194 dy2 = yf(iface) - yc(c2); 195 d2=sqrt(dx2^2+dy2^2); 196 wfun(iface) = (d2)*(d1 + d2); 197 gammaf(iface) = wfun(iface)*gamma(c1) + (1-wfun(iface))*gamma(c2); 198 end 199 200 %% Funciones interpolación Nodo-A-Vértice 201 202 cwfun = zeros(n_c,3); 203 sumwt = zeros(1,n_v); 204 d1=zeros(n_c,4); 205 for ic=1:n_c 206 for iv=1:3 207 v1 = Lcv(ic,iv); 208 d1(ic,iv) = sqrt((xv(v1) - xc(ic))^2 + (yv(v1) - yc(ic))^2); 209 % cwfun(ic,iv) = 1/d1; 210 sumwt(v1) = sumwt(v1) + 1/d1(ic,iv); 211 end 212 end 213 214 for ic=1:n_c 215 for iv=1:3 216 v1 = Lcv(ic,iv); 217 cwfun(ic,iv) = (1/d1(ic,iv))/sumwt(v1); 218 end 219 end 220 221 %% Cálculo de la solución 222 223 vn = Flujo_Potencial_UNST_ext_mesh(v_e); 224 Gf = rho*Cv*vn.*area(Lcf); 225 err = 10^3; % Error inicial para entrar en el bucle 226 TOL = 10^-3; % Tolerancia 227 phi = zeros(n_c,1); 228 229 flag = 0; 230 while max(err,[],’all’) > TOL 231 232 skew = zeros(1,n_c); 233 b=zeros(1,n_c); 234 sc=zeros(n_c); 235 236 A=zeros(n_c, n_c); 237 for ic=1:n_c 238 for k=1:3 239 iface = Lcf(ic,k); 240 if bface(iface) == 0 % Cara interior 241 icvec = Lfc(iface,:); 242 A(ic,icvec(icvec ~= ic)) = -gammaf(iface)*area(iface)/delta(iface) + 0.5*(Gf(ic,k) - abs(Gf(ic,k))); 243 A(ic,ic) = A(ic,ic) + 0.5*(Gf(ic,k) + abs(Gf(ic,k)));
120 Bibliografía [16] J. M. GORDILLO ARIAS DE SAAVEDRA and G. RIBOUX ACHER, Introducción a la Aerodinámica Potencial. Ediciones Paraninfo, S.A, 2012. [17] V. D. LISEIKIN, Grid Generation Methods, 3rd ed. Cham: Springer, 2017. [Online]. Available: https://link.springer.com/book/10.1007/978-3-319-57846-0 [18] T. K. D. S. W. CHENG and J. R. SHEWCHUK, Delaunay Mesh Generation. Boca Raton, FL: CRC Press, 2013. [19] P. MAJUMDAR, Computational Fluid Dynamics and Heat Transfer, 2nd ed. Boca Raton, FL: CRC Press, 2021. [20] W. Y. TEY, R. Y. W. HONG, Y. ASAKO, H. S. KANG, and K. C. NG, “Analysis on computational efficiency of convection discretisation schemes in simple algorithm,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 32, no. 3, pp. 989–1010, 2022. [21] S. BHATTACHARYYA, L. C. J. P. ABRAHAM, and J. GORMAN, “Introductory chapter: A brief history of and introduction to computational fluid dynamics,” in Applications of Computational Fluid Dynamics Simulation and Modeling, S. BHATTACHARYYA, Ed. Rijeka: IntechOpen, 2021, ch. 1. [Online]. Available: https://doi.org/10.5772/intechopen.97235 [22] M. KAWAGUTI, “Numerical solution to the navier-stokes equations for the flow around a circular cylinder at reynolds number 40,” Journal of the Physical Society of Japan, vol. 8, no. 6, pp. 747–757, 1953. [23] F. P. CARBALLO, Teoría de la Elasticidad, 1st ed. Sevilla, España: Editorial Científica, 2020. [24] AMWEL, “History of CFD,” Disponible en http://www.amwel.com/history.html.