Análisis y valoración del FVM en problemas de convección-difusión térmica
Abstract
El FVM (Finite Volume Method) es un método extendidamente utilizado en CFD (Computational Fluid Dynamics) para la resolución de problemas en Mecánica de Fluidos. Este proyecto pretende constituir un primer paso para el desarrollo de un código CFD en MATLAB que, aunque simplificado, contenga los ingredientes básicos para permitir entender a los estudiantes el funcionamiento interno de códigos CFD comerciales mucho más complejos (Fluent, Start CCM, etc.), que tienen que ser utilizados como “cajas negras” por el usuario. Con este objetivo, en este Trabajo Fin de Grado se analizará y valorará el método FVM para problemas de convección-difusión térmica en geometrías bidimensionales, usando mallas tanto estructuradas como no estructuradas. De particular interés resultará la implementación del mallado y la subsecuente obtención a partir del mismo de la información geométrica necesaria para resolver el problema, así como de los diferentes métodos upwind de diferentes órdenes necesarios para tratar el término convectivo de forma estable.
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 Análisis y valoración del FVM en problemas de convección-difusión térmica Autor: Alejandro Cumplido Díaz 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, 2022
Trabajo Fin de Grado Grado en Ingeniería Aeroespacial. Mención en Vehículos Aeroespaciales Análisis y valoración del FVM en problemas de convección-difusión térmica Autor: Alejandro Cumplido Díaz 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, 2022
Trabajo Fin de Grado: Análisis y valoración del FVM en problemas de conveccióndifusión térmica Autor: Alejandro Cumplido Díaz 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 En primer lugar, me gustaría agradecer a D. Miguel Pérez-Saborid Sánchez-Pastor por haber tutorizado de manera soberbia este Trabajo Fin de Grado. Su gran implicación, su ejemplar docencia y su extenso conocimiento sobre la materia han promovido enormemente mi motivación por el proyecto. Ha sabido ofrecerme las referencias literarias idóneas, que me han hecho sentir guiado en todo momento, al igual que ha sabido adaptar el alcance del proyecto según los obstáculos y limitaciones que he ido encontrando en el camino. En segundo lugar, también he de agradecer a toda la población de la Línea de Vuelo del A400M por la oportunidad brindada, pues éste ha sido el lugar donde me he iniciado en el mundo laboral. Aunque no tienen una relación directa con el proyecto, debido a que el desempeño de dicho trabajo ha tenido lugar en simultaneidad con la elaboración de este documento, veo necesario dedicar este espacio a todas esas personas que tanto me han enseñado, tanto en lo profesional como en lo personal. Especial mención merece mi equipo de Coordination Support, que me acogieron con los brazos abiertos y me formaron excelentemente. Por último, pero no menos importante, tengo que agradecer a todas esas personas que de una manera o de otra, directa o indirectamente, me han ofrecido su ayuda durante el trascurso de este trabajo. A toda mi familia, amigos y pareja, daros las gracias de corazón por haber estado siempre ahí, ofreciéndome un apoyo incondicional y resultando ser una vía de escape en los momentos más difíciles. Estoy seguro de que sin vosotros jamás lo hubiera conseguido. Alejandro Cumplido Díaz Escuela Técnica Superior de Ingeniería Sevilla, 2022 I
Resumen El FVM (Finite Volume Method) es un método extendidamente utilizado en CFD (Computational Fluid Dynamics) para la resolución de problemas en Mecánica de Fluidos. Este proyecto pretende constituir un primer paso para el desarrollo de un código CFD en MATLAB que, aunque simplificado, contenga los ingredientes básicos para permitir entender a los estudiantes el funcionamiento interno de códigos CFD comerciales mucho más complejos (Fluent, Start CCM, etc.), que tienen que ser utilizados como “cajas negras” por el usuario. Con este objetivo, en este Trabajo Fin de Grado se analizará y valorará el método FVM para problemas de convección-difusión térmica en geometrías bidimensionales, usando mallas tanto estructuradas como no estructuradas. De particular interés resultará la implementación del mallado y la subsecuente obtención a partir del mismo de la información geométrica necesaria para resolver el problema, así como de los diferentes métodos upwind de diferentes órdenes necesarios para tratar el término convectivo de forma estable. III
Índice de Figuras 1.1 Ejemplo de una simulación CFD para el análisis de una aeronave (imagen extraída de [7]) 5 1.2 Ejemplo de mallado estructurado para un perfil aerodinámico (imagen tomada de [3]) 6 1.3 Ejemplo de mallado no estructurado para un perfil aerodinámico (imagen tomada de [4]) 7 1.4 Definición del volumen de control 11 2.1 Mallado no homogéneo según el procedimiento A 14 2.2 Mallado no homogéneo según el procedimiento B 14 2.3 Representación del problema descrito 17 2.4 Solución proporcionada por el FVM para el problema descrito 19 2.5 Comparación entre la solución exacta y la obtenida mediante el FVM empleando 5 nodos equiespaciados 20 2.6 Comparación entre la solución exacta y la obtenida mediante el FVM empleando 25 nodos equiespaciados 20 2.7 Comparación entre la solución exacta y la obtenida mediante el FVM empleando 5 nodos ubicados aleatoriamente 21 2.8 Comparación entre la solución exacta y la obtenida mediante el FVM empleando 50 nodos con mayor concentración en los extremos 22 2.9 Comparación entre la solución exacta y la obtenida mediante el FVM para u=0.1m/s empleando 5 nodos equiespaciados y usando el método de las diferencias centradas 26 2.10 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 5 nodos equiespaciados y usando el método de las diferencias centradas 27 2.11 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 20 nodos equiespaciados y usando el método de las diferencias centradas 27 2.12 Problema paradigmático para ejemplificar la propiedad de conservatividad 28 2.13 Distribución de φ en la vecindad de dos fuentes. En la figura superior, el caso límite Pe →0. En la figura inferior, el caso límite Pe →∞30 2.14 Comparación entre la solución exacta y la obtenida mediante el FVM para u=0.1m/s empleando 5 nodos equiespaciados y usando el método upwind de primer orden 32 2.15 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 5 nodos equiespaciados y usando el método upwind de primer orden 32 2.16 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 25 nodos (con mayor densidad en los extremos) y usando el método upwind de primer orden 33 2.17 Nomenclatura de nodos en el método upwind de segundo orden 34 2.18 Comparación entre la solución exacta y la obtenida mediante el FVM para u=0.1m/s empleando 5 nodos equiespaciados y usando el método upwind de segundo orden 36 XI
XII Índice de Figuras 2.19 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 20 nodos equiespaciados y usando el método upwind de segundo orden 36 2.20 Comparación entre la solución exacta y la aproximación polinómica propuesta por Patankar 38 2.21 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 20 nodos equiespaciados y usando el esquema de ley de potencia de Patankar 38 2.22 Comparación entre la solución exacta y los distintos métodos de discretización para u=0.1m/sempleando 5 nodos equiespaciados 39 2.23 Comparación entre la solución exacta y los distintos métodos de discretización para u=3m/sempleando 5 nodos equiespaciados 39 2.24 Comparación entre la solución exacta y los distintos métodos de discretización para u=3m/sempleando 20 nodos equiespaciados 40 2.25 Geometría del problema que ejemplifica el caso no estacionario para t>042 2.26 Comparación entre la solución analítica y la solución mediante el FVM empleando 5 nodos, usando el esquema explícito para el término no estacionario y con un paso ∆t=2s43 2.27 Comparación entre la solución analítica y la solución mediante el FVM empleando 5 nodos equiespaciados, usando el esquema explícito para el término no estacionario y con un paso ∆t=10 sen t=40 s44 2.28 Comparación entre la solución analítica y la solución mediante el FVM empleando 20 nodos equiespaciados, usando el esquema explícito para el término no estacionario y con un paso ∆t=2sen t=40 s44 2.29 Comparación entre la solución analítica y la solución mediante el FVM empleando 20 nodos equiespaciados, usando el esquema explícito para el término no estacionario y con un paso ∆t=0.25 sen t=40 s45 2.30 Comparación entre la solución analítica y la solución mediante el FVM empleando 5 nodos equiespaciados, usando el esquema implícito para el término no estacionario y con un paso ∆t=2s46 2.31 Comparación entre la solución analítica y la solución mediante el FVM empleando 20 nodos equiespaciados, usando el esquema implícito para el término no estacionario y con un paso ∆t=2s46 3.1 Mallado 2D cartesiano 50 3.2 Representación del problema descrito 52 3.3 Mallado del dominio de computación representado en MATLAB. Los puntos rojos son los nodos donde se calcularán los valores de temperatura 55 3.4 Valores del campo de temperaturas en los nodos para un mallado uniforme M×N con M=4yN=455 3.5 Valores del campo de temperaturas en los nodos para un mallado uniforme M×N con M=50 yN=50 56 3.6 Curvas de nivel (isotermas) para un mallado uniforme M×Ncon M=50 yN=50 57 3.7 Representación del problema descrito 58 3.8 Valores del campo de temperaturas en los nodos para un mallado M×N con M=50 y N=10 60 3.9 Curvas de nivel (isotermas) para un mallado M×Ncon M=50 yN=10 61 4.1 Malla ortogonal cartesiana para aproximar un problema de flujo en un semicilindro 63 4.2 Indexación global en malla no estructurada 66 4.3 Representación de los vectores unitarios tangente y normal en una cara localmente con respecto a una celda 71 4.4 Esquema paradigmático para la obtención de las funciones interpoladoras wfywv,i74
Índice de Figuras XIII 4.5 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=3yN=384 4.6 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=50 yN=50 85 4.7 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=50 yN=50 86 4.8 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=30 yN=30 en el problema analizado en [6] y [17]-[20] 87 4.9 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=30 yN=30 en el problema analizado en [6] y [17]-[20] 89 4.10 Representación del problema descrito en [4] 89 4.11 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=30 yN=40 en el problema analizado en [4] 90 4.12 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=30 yN=40 en el problema analizado en [4] 91 4.13 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=50 yN=10 utilizando el esquema upwind 94 4.14 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=50 yN=10 95 4.15 Representación del problema descrito en [6] 96 4.16 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=50 yN=10 utilizando el esquema upwind en el problema analizado en [6] 97 4.17 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=50 yN=10 en el problema analizado en [6] 98
Índice de Tablas 2.1 Valores del campo de temperaturas y error relativo cometido por el FVM empleando 5 nodos equiespaciados 19 2.2 Valores del campo de temperaturas y error relativo cometido por el FVM empleando 5 nodos ubicados aleatoriamente 21 2.3 Valores de φ y comparación del error relativo cometido por el FVM usando el método FOU y SOU, para u=3m/sy empleando 7 nodos equiespaciados 36 2.4 Valores del campo de temperaturas y error relativo cometido por el FVM usando el esquema explícito, para ∆t=2sy empleando 5 nodos equiespaciados en t=40 s43 XV
Índice Abreviado Resumen III Abstract V Notación VII Índice de Figuras XI Índice de Tablas XV Índice Abreviado XVII 1 Introducción 1 1.1 Fundamentos del CFD 1 1.2 Historia del CFD 3 1.3 Estado del arte en CFD 4 1.4 Estructura de los códigos CFD 5 1.5 Ecuación de transporte 9 2 FVM para problemas convectivo-difusivos en 1D 13 2.1 Problema difusivo, estacionario y unidimensional 13 2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 24 2.3 Problema difusivo, no estacionario y unidimensional 40 3 FVM para problemas convectivo-difusivos en 2D. Mallado estructurado 49 3.1 Problema difusivo, estacionario y bidimensional 49 3.2 Problema convectivo-difusivo, estacionario y bidimensional 57 4 FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado 63 4.1 Problema difusivo, estacionario y unidimensional 64 4.2 Problema convectivo-difusivo, estacionario y unidimensional 91 5 Conclusión y trabajos futuros 101 Apéndice A Códigos de MATLAB 103 A.1 Convección-difusión estacionaria 1D 103 A.2 Difusión no estacionaria 1D 110 A.3 Difusión estacionaria 2D con mallado estructurado 115 A.4 Convección-difusión estacionaria 2D con mallado estructurado 120 XVII
XVIII Índice Abreviado Bibliografía 127
Índice Resumen III Abstract V Notación VII Índice de Figuras XI Índice de Tablas XV Índice Abreviado XVII 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 Historia del CFD 3 1.3 Estado del arte en CFD 4 1.4 Estructura de los códigos CFD 5 1.4.1 Pre-procesador 5 1.4.2 Solucionador 7 1.4.3 Post-procesador 8 1.5 Ecuación de transporte 9 2 FVM para problemas convectivo-difusivos en 1D 13 2.1 Problema difusivo, estacionario y unidimensional 13 2.1.1 Mallado del dominio unidimensional 13 2.1.2 Discretización de las ecuaciones 14 2.1.3 Resolución del sistema 16 2.1.4 Aplicación práctica 1D 17 2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 24 2.2.1 Esquema de las diferencias centradas 25 2.2.2 Propiedades de los esquemas de discretización 28 2.2.3 Esquema upwind de primer orden 31 2.2.4 Esquema upwind de segundo orden 34 2.2.5 Esquema de ley de potencia de Patankar 37 2.3 Problema difusivo, no estacionario y unidimensional 40 2.3.1 Esquema explícito 42 2.3.2 Esquema completamente implícito 45 XIX
6Capítulo 1. Introducción Es importante destacar que la variable en cuestión que se pretenda calcular como solución al problema de flujo objeto de estudio (temperatura, velocidad...) se obtiene en nodos ubicados en cada volumen de control definido. Así, la precisión del método estará supeditada al número de celdas en que se divide el domino, lo que es equivalente a decir que la solución será generalmente más precisa cuanto más fino sea el mallado propuesto. No obstante, un mallado más refinado requiere un mayor número de operaciones y por tanto un mayor coste computacional. Para llegar a un compromiso entre tiempo de computación y precisión de la solución será necesario definir un mallado adecuado. De forma general, existen dos maneras de realizar el mallado: de forma estructurada y de forma no estructurada. • Mallado estructurado. Consiste en emplear elementos típicamente cuadriláteros que conformen un diseño simple y adaptado a la geometría del problema en la división del dominio computacional. Normalmente, las celdas de los mallados estructurados pueden ser fácilmente identificadas por la terna {i,j,k} . Esto facilita enormemente los algoritmos al estar las celdas vecinas igualmente identificadas de forma sencilla. Otra ventaja es que la manera en que se guarda la información suele ser más directa y sencilla que en el caso no estructurado. Además, el hecho de disponer los elementos de forma alineada con la dirección del flujo facilita la convergencia del solver. El principal problema radica en que un mallado estructurado puede no adaptarse correctamente a la geometría cuando esta es compleja. En la Figura 1.2 puede encontrarse un ejemplo de un mallado estructurado, fácilmente reconocible por presentar un patrón que se repite en toda la malla. Figura 1.2 Ejemplo de mallado estructurado para un perfil aerodinámico (imagen tomada de [3]). • Mallado no estructurado. Cuando resulta difícil realizar un mallado estructurado debido a la complejidad de la geometría, es necesario recurrir a mallados no estructurados. Típicamente están conformados por celdas triangulares que no siguen un patrón fijo, y cuyas vecindades locales se disponen de manera arbitraria, inhabilitando la posibilidad de definirlas mediante la terna {i,j,k} . Esto provoca que el almacenamiento de la información sea más costoso, pues aparece la necesidad de guardar en la memoria la conectividad de las celdas. Además, al no
1.4 Estructura de los códigos CFD 7 estar alineadas con la dirección del flujo, la convergencia suele ser más lenta. En la Figura 1.3 puede observarse un mallado no estructurado. Son los más usuales en las aplicaciones CFD gracias a su adaptabilidad frente a geometrías complejas. Figura 1.3 Ejemplo de mallado no estructurado para un perfil aerodinámico (imagen tomada de [4]). Existe otro tipo de mallado que busca el compromiso entre los estructurados y los no estructurados. Son los conocidos como mallados estructurados por bloques, en los que el dominio se divide arbitrariamente en sub-regiones que, internamente, se comportan como mallados estructurados. Aunque se trata de una disposición de gran utilidad, queda fuera del alcance de esta memoria. En definitiva, algo común a cualquier estructura es que el mallado óptimo debe ser no uniforme, con una mayor densidad de nodos en las regiones donde la variable fluida de interés tiene variaciones grandes, y más espaciados cuando dichas variaciones son pequeñas. La tendencia actual está enfocada al desarrollo de códigos con una capacidad de adaptabilidad automática en función de las condiciones del problema, delegando en el propio código algunas de las decisiones que el usuario debe tomar en la definición del sistema, y que normalmente abarcan la mayor parte del tiempo en la industria del CFD. 1.4.2 Solucionador Usualmente conocido por su voz inglesa, solver, es el núcleo de la estructura de los códigos CFD y donde se encuentran embebidos los métodos numéricos que permiten lidiar con los problemas de flujo. Dentro de las técnicas de resolución numérica pueden distinguirse tres vertientes, ya descritas en el Capítulo 1.1: método de las diferencias finitas, método de los volúmenes finitos y método de los elementos finitos. Como ya se ha comentado, se centrarán los esfuerzos en el método de los volúmenes finitos (FVM), también denominado método de los volúmenes de control (CVM), un tipo de formulación en diferencias finitas de extendido uso en los códigos CFD habituales y cuyo algoritmo consta de los siguientes pasos:
8Capítulo 1. Introducción • Integración de las ecuaciones que gobiernan el problema en todos los volúmenes de control en los que se ha dividido el dominio fluido. Es este paso el que distingue al método de los volúmenes finitos de otras técnicas empleadas en CFD. La diferencia radica, como se ha comentado en otras ocasiones, en la relación entre el algoritmo numérico y los principios de conservación de la física del problema. • Discretización de las ecuaciones integrales para generar un sistema de ecuaciones algebraico. Como se verá más adelante, existen numerosas técnicas que implementan los códigos CFD para el tratamiento de los términos convectivos, difusivos y transitorios. • Resolución del sistema de ecuaciones algebraico mediante el uso de métodos iterativos. Existe toda una rama de las matemáticas dedicada al estudio de algoritmos que optimizan la resolución de sistemas de ecuaciones algebraicos. No obstante, este análisis queda fuera del alcance de este documento. En efecto, la prevalencia de este método frente a otras técnicas CFD estriba en la integración sobre los volúmenes de control, estableciendo la conservación de las magnitudes de interés en cada célula que compone el mallado. Así, el algoritmo empleado no es más que una manifestación de los principios físicos de conservación subyacentes, lo que facilita la comprensión de este método en comparación con los métodos diferenciales o de elementos finitos. La conservación de una variable fluídica genérica (que en adelante se denominará por la letra φ ) en un volumen de control puede entenderse como un balance neto entre los procesos que tienden a incrementar la variable y los que tienden a disminuirla. Sin pérdida de generalidad, este balance puede expresarse (con palabras) de la siguiente forma (véase [4]): Tasa de cambio de φcon respecto al tiempo en el volumen de control = Incremento neto de φdebido a convección en el volumen de control + Incremento neto de φdebido a difusión en el volumen de control + Incremento neto de φdebido a generación en el volumen de control Dentro de los códigos CFD se encuentran embebidos las técnicas y modelos apropiados para el tratamiento de los fenómenos de transporte, a decir: convección (transporte asociado al flujo), difusión (transporte asociado a las variaciones de φ entre diferentes puntos), generación (transporte asociado a la creación o destrucción de la magnitud φ ) y el transporte asociado a la tasa de cambio con respecto al tiempo. Estos modelos usualmente requieren resolver sistemas lineales mediante métodos iterativos. Son muchos y variados los algoritmos que proporcionan resultados adecuados a la hora de resolver este tipo de sistemas, y su desarrollo abarcaría un capítulo completo. Por ello, no se ahondará en el presente informe en el estudio de dichos algoritmos, sino que se empleará el comando \ de MATLAB para resolver los sistemas del tipo Ax =b⇒x=A\b . Esto puede hacerse sin cavilaciones puesto que MATLAB posee algoritmos internos para detectar automáticamente la estructura de la matriz Ay aplicar el método más apropiado. 1.4.3 Post-procesador Una vez llevada a cabo la simulación, el paso último del análisis CFD consiste en la adquisición, tratamiento e interpretación de la información obtenida. Para poder extraer conclusiones relevantes en relación con el estudio llevado a cabo es de vital importancia la elección de las variables que se quieren graficar. De forma parecida a lo comentado en el pre-procesamiento, el post-procesador deberá contar con aptitudes similares, disponiendo de potentes herramientas gráficas de visualización
1.5 Ecuación de transporte 9 de datos que faciliten la interpretación de los resultados y la exportación de los mismos para su uso más allá de los códigos. Así, es imprescindible que el post-procesador cuente con la capacidad de manipular las superficies 2D o 3D generadas y mostrar su mallado o representar los resultados mediante campos vectoriales apoyados en códigos de colores o isolíneas, entre otras funciones. Aunque, como ya se ha comentado, más del 50% del tiempo empleado en la industria de los códigos CFD se consume en la generación del mallado y la definición de la geometría, esta actividad compete enteramente al usuario y estará condicionada por la experiencia y habilidad que éste presente en el manejo de pre-procesadores y post-procesadores CFD. Como este asunto queda fuera de alcance, en este documento se estará particularmente interesado en el solver, por ser el elemento en el que se encuentran embebidos los algoritmos que permiten la resolución numérica del problema. 1.5 Ecuación de transporte Las ecuaciones que gobiernan el flujo de fluidos se infieren de las ecuaciones de conservación paradigmáticas de la física: •Conservación de la masa. La masa de fluido se conserva. • Conservación de la cantidad de movimiento. La variación de la cantidad de movimiento es igual a la suma de todas las fuerzas en la partícula fluida (Segunda Ley de Newton). • Conservación de la energía. La variación de la energía es igual a la suma de las variaciones de calor y trabajo en la partícula fluida (Primera Ley de la Termodinámica). El desarrollo físico-matemático de dichas expresiones aplicado a una partícula fluida lleva a la formulación del sistema de ecuaciones que gobierna el problema tridimensional de transferencia de flujo y calor en fluidos newtonianos compresibles, las extendidamente conocidas ecuaciones de Navier-Stokes. Dicha formulación es válida siempre que se conciba el fluido como continuo. Esto es, la descripción del comportamiento del fluido según sus propiedades macroscópicas (densidad, velocidad, presión, temperatura...), ignorando los efectos de estructura y movimiento molecular. Así, una partícula fluida será el elemento de fluido más pequeño posible cuyas propiedades macroscópicas no se ven influenciadas por el comportamiento molecular. Sin entrar en su demostración, las ecuaciones de Navier-Stokes presentan la siguiente forma para un fluido newtoniano compresible (véase [2]): Ecuación de continuidad: ∂ρ ∂t+∇·(ρv) = 0(1.1) Ecuación de cantidad de movimiento: ρ∂v ∂t+ρv·∇v=−∇p+∇·(µ∇v) +ρ fm(1.2) Ecuación de la energía: ρ∂e ∂t+ρv·∇e=−p∇·v+τ′:∇v+Qr+Qq+∇·(k∇T)(1.3) Adicionalmente bajo la suposición de fluido ideal y equilibrio termodinámico, se tiene:
10 Capítulo 1. Introducción Ecuación de los gases ideales: p=p(ρ,T)e=e(ρ,T)→p ρ=RgT e =cvT(1.4) En estas expresiones, t denota el tiempo, fm es el vector de fuerzas másicas, cv es la capacidad calorífica del gas a volumen constante, k su conductividad térmica, τ′:∇v es el término de disipación de energía cinética interna, Qr y Qq son las potencias caloríficas que, por unidad de volumen, recibe el fluido por radiación y reacción química, respectivamente, Rg=R Mm=cp−cv=cv(γ−1) , R es la constante universal de los gases, cp la capacidad calorífica del gas a presión constante, e es la energía interna del fluido, Mm la masa molecular del gas y γ el coeficiente de dilatación adiabática (para el aire, Rg≈287 J KgK yγ=1.4). Así, se conforma un sistema de siete ecuaciones (debido a que la ecuación de cantidad de movimiento es una ecuación vectorial que equivale a tres ecuaciones escalares) y siete incógnitas ( p,T,ρ,v,e ) a las que, adicionalmente, habrá que imponer unas condiciones iniciales y de contorno adecuadas. Debido a la similitud que presentan entre sí las ecuaciones de Navier-Stokes, es posible reescribirlas en una única ecuación que representa la evolución de una magnitud o propiedad genérica, denotada por φ . Se llega así a lo que se conoce como ecuación de transporte para la propiedad φ : ∂(ρφ) ∂t+∇·(ρφv) = ∇·(Γ∇φ)+Sφ(1.5) En el primer miembro aparecen los términos de variación temporal y transporte convectivo, y en el segundo los términos de transporte difusivo y generación, respectivamente. Esta ecuación constituye un punto de partida para los algoritmos computacionales convencionales en el FVM. El paso clave de este método consiste en integrar dicha expresión en un volumen de control tridimensional (en el caso más genérico), resultando: ZVC ∂(ρφ) ∂tdV +ZVC ∇·(ρφv)dV =ZVC ∇·(Γ∇φ)dV +ZVC SφdV (1.6) Reescribiendo los términos convectivo y difusivo como integrales de superficie cerrada (aquélla que contiene al volumen de control) mediante la aplicación del Teorema de la Divergencia de Gauss, la ecuación (1.6) puede reescribirse como: ∂ ∂tZVC ρφ dV +ZA n·(ρφv)dA =ZA n·(Γ∇φ)dA +ZVC SφdV (1.7) Esta expresión tiene un claro significado físico que puede ser fácilmente interpretado analizando individualmente cada uno de los términos. El primero de ellos representa la tasa de cambio total de la cantidad de propiedad φ que hay en el volumen de control. El producto escalar n·(ρφv) expresa la componente de flujo de φ en la dirección de la normal exterior de la superficie que encierra al volumen de control. Así, el segundo término está asociado al decremento neto de φ en el elemento fluido debido a convección. El flujo difusivo es por definición positivo en la dirección negativa del gradiente, por lo que el producto n·(Γ∇φ) = −n·(−Γ∇φ) puede interpretarse como la componente del flujo difusivo positivo en la dirección del vector normal interior de la superficie del volumen de control. Por tanto, el primer término del segundo miembro representa el incremento neto de φ en el elemento fluido debido a difusión. El último término de la expresión simboliza el incremento de φ debido a las fuentes de generación de φ internas en el elemento fluido. Esto pone de manifiesto con mayor claridad lo que ya se había declarado con anterioridad: la integración de la ecuación de transporte equivale a una ecuación de conservación de la propiedad fluida φ en
1.5 Ecuación de transporte 11 un volumen de control macroscópico y finito. Este elemento fluido puede caracterizarse según la siguiente figura: E W N T B S Figura 1.4 Definición del volumen de control. El caso más general posible es un volumen de control cúbico acotado por seis caras. Siguiendo la nomenclatura de obras como [ 1 ]o[ 4 ], se denotarán estas caras por N , S , E , W , T y B , haciendo referencia a North, South, East, West, Top yBottom, respectivamente. Esta denominación se hace según la dirección positiva de los ejes cartesianos, también representada en la Figura 1.4. Cuando se haga referencia a una celda dentro de un mallado, se usará P para nombrar al nodo de la celda en cuestión, N , S , E , W , T y B para denotar a los nodos de las celdas vecinas, y n , s , e , w , t y b para las caras entre ellas. Es interesante resaltar que P= (x,y,z) estará ubicado normalmente en el baricentro del volumen de control de forma que la magnitud de dicha celda sea lo más representativa posible. Cabe destacar además que todas las propiedades fluidas son función del tiempo y del espacio, por lo que estrictamente hablando habría que escribir φ=φ(x,y,z,t) en todo momento. No obstante, por economía de la escritura, se dará por hecho dicha dependencia funcional sin necesidad de expresarla. En los capítulos venideros se aplicarán las expresiones desarrolladas a problemas convectivo-difusivos prácticos y se obtendrán y analizarán los resultados. Es necesario destacar la fuerte influencia que en adelante tendrán las obras de Versteeg ([ 4 ]) y Mazumder ([ 6 ]) a la hora de explanar los principales capítulos de este proyecto. Sin embargo, una aportación importante de este TFG ha sido el desarrollo de programas en MATLAB propios que permiten refinar arbitrariamente el mallado en los problemas y ejemplos mostrados. Los códigos más representativos se han incluido en la memoria, mientras que los que pueden resultar redundantes o captan menos fielmente la esencia del FVM están contenidos en el Apéndice ??. Además, para poder comparar los resultados de los programas propios con los de la literatura, se han reproducido en ocasiones algunos de los problemas (con los mismos datos) que se incluyen en estas obras.
2 FVM para problemas convectivo-difusivos en 1D Una vez definida la naturaleza de la ecuación de transporte se va a desarrollar el método numérico basado en su integración, a decir, el FVM. Primeramente va a considerarse el proceso de transporte más simple posible: difusión estacionaria pura en el caso unidimensional. Posteriormente y de forma progresiva, se irán introduciendo los términos convectivos y no estacionarios hasta llegar al estudio de un problema de carácter general para el caso unidimensional. Es de vital importancia resaltar la fuerte influencia que ha tenido la obra de Versteeg y Malalasekera ([ 4 ]) en el orden de exposición del capítulo, las ilustraciones mostradas y los ejemplos exhibidos. 2.1 Problema difusivo, estacionario y unidimensional En el caso estacionario y en ausencia de convección, la integración de la ecuación de transporte queda: ZA n·Γdφ dx exdA +ZVC SφdV =0(2.1) Para caracterizar y resolver el problema hay que realizar fundamentalmente tres pasos: generación del mallado, discretización de las ecuaciones y resolución del sistema. A continuación se describirán sucintamente cada una de estas etapas. Aunque para ilustrar esta explicación servirá de ejemplo el caso estacionario sin convección, este desarrollo es extrapolable al resto de problemas. La razón es que la estructura, que no es más que un mero reflejo de la que se expuso en el Capítulo 1.4, es idéntica en todos los problemas. 2.1.1 Mallado del dominio unidimensional En primer lugar, es necesario dividir el dominio en volúmenes de control finitos en los cuales discretizar las ecuaciones. En el caso unidimensional, cada uno de ellos quedará completamente caracterizado por un nodo y dos caras laterales. En casos multidimensionales, de forma tradicional bastará con un nodo y cuatro o seis caras laterales si se utilizan cuadriláteros o hexaedros, respectivamente. Es habitual definir los volúmenes de control ubicados en los extremos del dominio de forma que sus caras laterales coincidan con los límites físicos del dominio, facilitando así la implementación de las condiciones de contorno. Existen distintas prácticas para disponer estos elementos según se quiera generar un mallado uniforme o no uniforme. Como ya se ha comentado anteriormente, lo más usual en la resolución 13
14 Capítulo 2. FVM para problemas convectivo-difusivos en 1D de la mayoría de problemas en ingeniería es el empleo de mallas no estructuradas, con mayor refinamiento en regiones donde las variaciones de la variable en estudio son más acusadas. Aún en mallas estructuradas, la disposición de las caras y los nodos puede prescindir de homogeneidad en busca de una malla más eficiente, refinada en zonas con gradientes elevados. Citando [ 1 ], son dos los procedimientos más habituales para la generación de mallas no homogéneas: • Procedimiento A. Se definen primero los nodos y después se localizan las caras de los volúmenes de control en los puntos medios entre nodos. Este procedimiento se ilustra en la Figura 2.1. W P E e w Volumen de control Figura 2.1 Mallado no homogéneo según el procedimiento A. • Procedimiento B. Se ubican primero las caras de los volúmenes de control y después se colocan los nodos en los puntos medios entre caras. Este procedimiento se ilustra en la Figura 2.2. W P E e w Volumen de control Figura 2.2 Mallado no homogéneo según el procedimiento B. Es más recomendable optar por la práctica B en detrimento del primer procedimiento cuando se quiere trabajar con mallas no uniformes. La justificación de dicha elección, junto con un análisis detallado de las ventajas y desventajas de ambos métodos, puede encontrarse en [ 1 ]. Sin pormenorizar los detalles, puede decirse que en la práctica A los métodos de aproximación mediante diferencias centradas para computar los valores de una propiedad en las caras son más precisos que en la práctica B, pues las caras se encuentran en los puntos medios entre nodos. Por otro lado, en la práctica B los valores de φ en los nodos son más representativos al localizarse el nodo en el centro de la celda, cosa que no sucede en la práctica A. No obstante, por simplicidad, los capítulos venideros se desarrollarán en base a una malla uniforme (nodos equiespaciados). 2.1.2 Discretización de las ecuaciones El siguiente paso consiste en particularizar la ecuación integral de transporte para cada uno de los volúmenes de control, resultando:
2.1 Problema difusivo, estacionario y unidimensional 15 Z∆V d dx Γdφ dx dV +Z∆V SdV =ΓeAene dφ dx e +ΓwAwnw dφ dx w +¯ S∆V=0(2.2) Esta expresión pone de manifiesto la característica antológica y de fácil interpretación física que diferencia al método de los volúmenes finitos del resto de métodos numéricos y que ya ha sido citada en ocasiones anteriores; la expresión (2.2) representa una ecuación de balance en el volumen de control. Nótese que, aunque la solución del problema se obtiene en los nodos, éstas ecuaciones discretizadas requieren evaluar los parámetros característicos del sistema y los gradientes en las caras de los volúmenes de control. Para ello, es necesario desarrollar un modelo que permita aproximar correctamente estos valores. La forma más simple de acometer esta tarea pasa por utilizar una interpolación lineal entre nodos para determinar los valores en las caras. Definiendo las funciones de interpolación gwyge, estos valores pueden expresarse de la siguiente manera: Γw=gwΓP+(1−gw)ΓWdonde gw=xw−xW xP−xW (2.3) Γe=geΓP+(1−ge)ΓEdonde ge=xE−xe xE−xP (2.4) En el caso de trabajar con una malla uniforme, gw=ge=0.5 y el valor en las caras no será más que la semisuma de los valores en los nodos adyacentes a la cara. Es decir: Γw=ΓP+ΓW 2Γe=ΓP+ΓE 2(2.5) Asimismo, para el caso de los términos difusivos se empleará lo que se conoce como método de las diferencias centradas, una aproximación de segundo orden (solo si la malla es uniforme) que da buenos resultados gracias a la propia naturaleza de la difusión, que se transporta homogéneamente en todas direcciones. Así, estos términos pueden expresarse como: dφ dx w =φP−φW xP−xW dφ dx e =φE−φP xE−xP (2.6) Introduciendo estas expresiones en la ecuación (2.2) se obtiene: −ΓeAene φE−φP xE−xP −ΓwAwnw φP−φW xP−xW =¯ S∆V(2.7) Definiendo De=ΓeAene/(xE−xP) y Dw=ΓwAwnw/(xP−xW) como términos asociados al flujo difusivo y reordenando se llega a: DwφW+(De−Dw)φP−DeφE=¯ S∆V(2.8) Que puede expresarse como: aWφW+aPφP+aEφE=aSP (2.9) El término aSP depende del modelo de fuente que se utilice. Es importante resaltar que esta ecuación discretizada representa la expresión típica alrededor de la cual girarán todos los modelos que se desarrollarán a lo largo de la memoria.
22 Capítulo 2. FVM para problemas convectivo-difusivos en 1D reproducir correctamente posibles gradientes elevados de φ en dichos límites (por ejemplo, su uso es fundamental en el estudio de la capa límite). Esto se ejemplifica a continuación: 0 0.005 0.01 0.015 0.02 0.025 0.03 50 100 150 200 250 300 350 400 Figura 2.8 Comparación entre la solución exacta y la obtenida mediante el FVM empleando 50 nodos con mayor concentración en los extremos. Para que sirva de ejemplo al lector, se introduce a continuación el código MATLAB desarrollado para resolver el problema descrito. En dichas líneas puede apreciarse fácilmente la construcción de la matriz A según lo establecido en la Sección 2.1.3, así como la generación de la malla, el post-procesamiento y la gestión de la información. 1%% PROBLEMA DE DIFUSION PURA 2 3clear all;close all;clc 4 5TablaErrores = []; 6 7% MALLADO 8 9L = 0.03; % Longitud placa (m) 10 N = 5; % Número de Nodos 11 12 caras = linspace(0,L,N+1); % Nodos Equiespaciados 13 14 % caras_theta = linspace(0,pi,N+1); % Nodos Concentrados en Extremos 15 % caras = L/2*(1 - cos(caras_theta)); 16 17 % rng = [0,L*rand(1,N-1),L]; % Nodos Aleatorios 18 % caras = sort(rng); 19 20 for i = 2:length(caras) 21 nodos(i-1) = (caras(i) - caras(i-1))/2 + caras(i-1); 22 end 23
2.1 Problema difusivo, estacionario y unidimensional 23 24 % DATOS 25 26 Gamma_w = 0.75*ones(N,1); % Conductividad Térmica (W/(m*K)) -> CTE 27 Gamma_e = 0.75*ones(N,1); % Conductividad Térmica (W/(m*K)) -> CTE 28 Aw = 1*ones(N,1); % Ãrea Transversal (m^2) -> CTE 29 Ae = 1*ones(N,1); % Ãrea Transversal (m^2) -> CTE 30 nw = -ones(N,1); % Vector Normal (nw=-1 cte para caras west) 31 ne = ones(N,1); % Vector Normal (ne=1 cte para caras east) 32 Fw=zeros(N,1); % Flujo (F=0 por ser problema de difusión pura) 33 Fe=zeros(N,1); % Flujo (F=0 por ser problema de difusión pura) 34 q = 1500e3*ones(N,1); % Generación Volumétrica de Calor (W/m^3) 35 TA = 50+273; % CdC Izquierda: phi impuesto (K) 36 TB = 250+273; % CdC Derecha: phi impuesto (K) 37 38 % COEFICIENTES 39 40 for I = 2:N-1 41 gw(I) = (caras(I) - nodos(I-1))/(nodos(I) - nodos(I-1)); % Factor Interpolación 42 ge(I) = (caras(I+1) - nodos(I+1))/(nodos(I) - nodos(I+1)); % Factor Interpolación 43 Dw(I) = Gamma_w(I)*nw(I)*Aw(I)/(nodos(I) - nodos(I-1)); 44 De(I) = Gamma_e(I)*ne(I)*Ae(I)/(nodos(I+1) - nodos(I)); 45 aW(I) = Fw(I)*(1 - gw(I)) + Dw(I); % Coeficiente phi_W 46 aP(I) = Fw(I)*gw(I) + Fe(I)*ge(I) - Dw(I) + De(I); % Coeficiente phi_P 47 aE(I) = Fe(I)*(1 - ge(I)) - De(I); % Coeficiente phi_E 48 aSP(I) = q(I)*Ae(I)*(caras(I+1) - caras(I)); % Término Asociado a la Generación 49 A(I,I-1) = aW(I); 50 A(I,I) = aP(I); 51 A(I,I+1) = aE(I); 52 b(I) = aSP(I); 53 end 54 55 % CONDICIONES DE CONTORNO 56 57 % LEFT 58 59 phiL = TA; % Temperatura en Cara Izquierda (K) 60 gw(1) = 0; 61 ge(1) = (caras(2) - nodos(2))/(nodos(1) - nodos(2)); 62 Dw(1) = Gamma_w(1)*nw(1)*Aw(1)/(nodos(1) - caras(1)); 63 De(1) = Gamma_e(1)*ne(1)*Ae(1)/(nodos(2) - nodos(1)); 64 aW(1) = Fw(1)*(1 - gw(1)) + Dw(1); 65 aP(1) = Fw(1)*gw(1) + Fe(1)*ge(1) - Dw(1) + De(1); 66 aE(1) = Fe(1)*(1 - ge(1)) - De(1); 67 aSP(1) = q(1)*Ae(1)*(caras(2) - caras(1)); 68 A(1,1) = aP(1); 69 A(1,2) = aE(1); 70 b(1) = aSP(1) - aW(1)*phiL; 71 72 % RIGHT 73 74 phiR = TB; % Temperatura en Cara Derecha (K) 75 gw(N) = (caras(N) - nodos(N-1))/(nodos(N) - nodos(N-1)); 76 ge(N) = 0; 77 Dw(N) = Gamma_w(N)*nw(N)*Aw(N)/(nodos(N) - nodos(N-1)); 78 De(N) = Gamma_e(N)*ne(N)*Ae(N)/(caras(end) - nodos(N)); 79 aW(N) = Fw(N)*(1 - gw(N)) + Dw(N); 80 aP(N) = Fw(N)*gw(N) + Fe(N)*ge(N) - Dw(N) + De(N); 81 aE(N) = Fe(N)*(1 - ge(N)) - De(N); 82 aSP(N) = q(N)*Ae(N)*(caras(end) - caras(end-1)); 83 A(N,N-1) = aW(N); 84 A(N,N) = aP(N); 85 b(N) = aSP(N) - aE(N)*phiR;
24 Capítulo 2. FVM para problemas convectivo-difusivos en 1D 86 87 % SOLUCION 88 89 phi = A\b’ - 273 90 phiexact = @(x)(((TB - TA)/L + q(1)/2/Gamma_e(1)*(L-x)).*x+TA - 273); 91 plot([0,nodos,L],[TA - 273;phi;TB - 273],’ro-’,’MarkerIndices’,2:length(nodos) + 1,’ linewidth’,1.2) 92 hold on 93 x = 0:0.0002:L; 94 plot(x,phiexact(x),’k’,’linewidth’,1.5) 95 hold on 96 xlabel(’$m$’,’Interpreter’,’latex’) 97 ylabel(’$T(^\circ C$)’,’Interpreter’,’latex’) 98 title(’Campo de Temperaturas’,’Interpreter’,’latex’) 99 legend({’Soluci\’’on con FVM’,’Soluci\’’on exacta’},’Interpreter’,’latex’,’Location’,’ Southeast’) 100 101 for i = 1:N 102 TablaErrores = [TablaErrores;i,phi(i),phiexact(nodos(i)),... 103 abs((phi(i) - phiexact(nodos(i)))/phiexact(nodos(i)))*100]; 104 end 105 106 disp(’ Nodo Sol. FVM Sol. Exacta Error(%)’) 107 disp(’ ’) 108 disp(TablaErrores) En el siguiente capítulo, se introducirán los efectos de advección (más comúnmente conocido como convección) a los problemas de flujo y se implementarán y discutirán las diferentes técnicas de modelización que existen a la hora de implementarlos en el FVM. 2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización Hasta ahora se han estudiado únicamente los problemas difusivos, aquellos en los que no existe flujo de fluido. Sin embargo, en la naturaleza son de mucha mayor abundancia los sistemas en los que difusión y convección ocurren simultáneamente, por lo que surge la necesidad de desarrollar nuevos métodos que permitan predecir correctamente ambos efectos en consonancia y aumentar la aplicabilidad del FVM. Por ello, en adelante se analizarán sistemas en los que el flujo de fluido juega un papel significativo y se describirán las distintas técnicas que pueden emplearse para tratar los términos convectivos que, como se verá, dan lugar en ocasiones a inestabilidades en la solución. La integración de la ecuación de transporte para el caso convectivo-difusivo estacionario y unidimensional sobre el volumen de control (que representa un balance de flujo sobre éste) resulta: ZA n·(ρφv)dA =ZA n·Γdφ dx exdA +ZVC SφdV (2.26) El principal problema en la discretización de los términos convectivos de la ecuación radica en la determinación de la propiedad φ en las caras de los volúmenes de control y en el cálculo de los flujos a través de éstas. En el análisis del problema difusivo se introdujo el método de las diferencias centradas para evaluar los términos difusivos en las caras y se demostró que funcionaba bien incluso para un bajo número de celdas, por lo que puede resultar lógico extender su uso a los términos convectivos. No obstante, existe una diferencia importante que limita la extensión de este método y que se describe a continuación.
2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 25 Mientras que en los procesos difusivos la distribución de la propiedad φ es transportada de forma homogénea en todas las direcciones, en los procesos de convección la influencia se limita a la dirección del flujo. Como se comprobará más adelante, esta disimilitud se materializa en una restricción en el número de celdas necesarias para obtener una solución válida si se usa el método de las diferencias centradas para modelar la convección. Dicho número está sujeto a la importancia relativa que tengan los términos convectivos frente a los difusivos, siendo necesario un mallado más fino para tratar fuertes convecciones. Es por ello que se presentarán métodos alternativos para obtener soluciones estables sin necesidad de un mallado tan restrictivo, que no siempre es posible. Adicionalmente, se supondrán conocidos los campos de velocidades en las caras de los volúmenes de control y, por tanto, los flujos. La forma de resolver el problema es conceptualmente idéntica al caso de difusión pura, por lo que se omitirá el desarrollo exhaustivo de las fases para evitar ser redundantes. Se compone de los mismos tres pasos: mallado, discretización y resolución. Una vez definidos los volúmenes de control, integrando sobre cada uno de ellos se llega a la siguiente ecuación: ρwuwAwnwφw+ρeueAeneφe−ΓeAene φE−φP xE−xP −ΓwAwnw φP−φW xP−xW =¯ S∆V(2.27) Donde se ha utilizado la diferenciación central para tratar los gradientes difusivos. Además, el flujo debe cumplir la ecuación de continuidad (ecuación de conservación). Es decir: ZA n·(ρv)dA =0⇒Fw+Fe=0(2.28) Para obtener el sistema de ecuaciones discretizadas que gobierna el problema de forma más compacta resulta conveniente definir Fw=ρwuwAwnw y Fe=ρeueAene como términos asociados al flujo convectivo y Dw=ΓwAwnw/(xP−xW) y De=ΓeAene/(xE−xP) asociados al transporte difusivo. Así, la ecuación que debe satisfacerse en cada nodo será: Feφe+Fwφw−De(φE−φP)−Dw(φP−φW) = ¯ S∆V(2.29) Que puede expresarse fácilmente como la ecuación (2.18) sin más que reagrupar términos. Las próximas secciones se dedicarán al desarrollo, análisis, implementación y valoración de algunos de los distintos métodos que eminentemente se utilizan en los códigos comerciales para evaluar la propiedad φ en las paredes de los volúmenes de control. Se señalarán las ventajas y desventajas de cada uno de ellos, y las limitaciones que han llevado a la necesidad de recurrir a nuevos métodos. Para dar cierta linealidad a la exposición y siendo fiel a [ 4 ], éstos serán presentados en el orden lógico en que fueron surgiendo. 2.2.1 Esquema de las diferencias centradas Este método ya ha sido presentado anteriormente, pues es el que se está usando por defecto para evaluar los términos difusivos. Su fundamento no es otro que el de interpolar linealmente los valores de los nodos adyacentes a la cara que se quiere evaluar. Para una malla no uniforme, las expresiones serían análogas a las expuestas en (2.3) y (2.4), siendo necesario el cálculo de las funciones de interpolación gw y gw . Para un mallado uniforme, las expresiones que resultan son nuevamente la semisuma de los valores: φw=φW+φP 2φe=φP+φE 2(2.30)
26 Capítulo 2. FVM para problemas convectivo-difusivos en 1D Sustituyendo dichas expresiones en la ecuación (2.29), se obtiene: Fw 2+DwφW+Fw+Fe 2+De−DwφP+Fe 2−DeφE=¯ S∆V(2.31) Puede reconocerse fácilmente la misma estructura general que en (2.18). A modo ilustrativo y para que sirva de comparación entre los diferentes métodos, en adelante se tomará como problema de convección-difusión paradigmático el que se describe a continuación, esencialmente idéntico al presentado en [ 4 ]. Sea una propiedad φ transportada por medio de difusión de cantidad de movimiento (difusión viscosa) y convección a través de un dominio unidimensional de longitud L=1m. La ecuación que gobierna este sistema será: d dx(ρuφ) = d dx Γdφ dx (2.32) Sujeto a las condiciones de contorno φ0=1 en x=0 y φL=0 en x=L . La velocidad se considerará uniforme y de valor conocido u∈[0.1,3](m/s) . Como se verá más adelante, el valor de u será determinante para caracterizar la aplicabilidad de los distintos métodos, pues a mayor velocidad, mayor será la convección. La difusión vendrá caracterizada por Γ=0.1kg/(ms) . Además, se considerará que el fluido en cuestión es agua ( ρ=1kg/m3 ). Utilizando las herramientas adecuadas de la teoría de resolución de ecuaciones diferenciales ordinarias, puede demostrarse que la solución analítica de dicha ecuación es: φ(x) = φ0+(φL−φ0)eρux/Γ−1 eρuL/Γ−1(2.33) Se muestra a continuación una comparación entre la solución exacta y la solución que proporciona el FVM (utilizando diferenciación central para tratar los términos convectivos y difusivos en las caras) bajo distintas condiciones. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2.9 Comparación entre la solución exacta y la obtenida mediante el FVM para u=0.1m/s empleando 5 nodos equiespaciados y usando el método de las diferencias centradas.
2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 27 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.5 0 0.5 1 1.5 2 2.5 3 Figura 2.10 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 5 nodos equiespaciados y usando el método de las diferencias centradas. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2.11 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 20 nodos equiespaciados y usando el método de las diferencias centradas. Son varios los comentarios interesantes que pueden hacerse acerca de estas representaciones. En primer lugar nótese cómo, para un valor bajo de la velocidad del campo fluido, un mallado poco refinado es suficiente para obtener una buena aproximación de la solución. A medida que la velocidad se incrementa, el término convectivo prepondera sobre el difusivo, apareciendo unas
28 Capítulo 2. FVM para problemas convectivo-difusivos en 1D oscilaciones en la solución aproximada que la hacen diferir notablemente de la solución analítica. Para resolver esta problemática, puede comprobarse cómo un aumento en el número de celdas es suficiente para que la solución numérica sea razonablemente parecida a la exacta. La conclusión que puede inferirse de este análisis es que a medida que el flujo convectivo cobra importancia sobre el difusivo, que es equivalente a decir que la relación F/D aumenta, la diferenciación centrada produce resultados cada vez menos precisos. Esta disparidad se materializa en oscilaciones sobre la solución exacta, que pueden atenuarse si se disminuye el ratio F/D , por ejemplo, refinando la malla. La razón de que aparezcan estas oscilaciones y la influencia que tiene en ello el factor F/Dserá tratado a continuación. 2.2.2 Propiedades de los esquemas de discretización Para poder desarrollar el resto de métodos de discretización de forma inteligible, es necesario estudiar en profundidad las razones por las que la diferenciación centrada no proporciona buenos resultados bajo determinadas condiciones en problemas convectivo-difusivos. Más aún, aunque los resultados numéricos deben coincidir exactamente con la solución analítica para un número de celdas infinitamente grande (independientemente del método de discretización utilizado), esto no es posible en la práctica. La necesidad de obtener resultados físicamente realistas incluso con un número limitado de celdas induce a analizar ciertas propiedades que subyacen a los métodos de discretización y que éstos deben tener. Algunas de estas propiedades fundamentales son: conservatividad,delimitación ytransportividad. Estas características serán explicadas a continuación. Destacar que la exposición y las ilustraciones han sido fuertemente influenciadas por [4]. Conservatividad Tal y como se ha comentado en ocasiones anteriores, la conservatividad es uno de los pilares que sustenta la formulación del FVM. La integración de la ecuación de transporte sobre cada volumen de control lleva a un conjunto de ecuaciones de conservación de la propiedad φ . Para garantizar dicha conservación, tanto global como localmente, el flujo de φ saliente a través de una determinada 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. Esto se consigue representando los valores de flujo de manera consistente y unívoca. Para ejemplificar esto, considérese el siguiente problema de difusión pura sin generación, estacionario y unidimensional: 1234 Figura 2.12 Problema paradigmático para ejemplificar la propiedad de conservatividad. Empleando las diferencias centradas para tratar los gradientes difusivos y sumando las ecuaciones discretizadas en cada volumen de control, se llega a la ecuación de balance neto global: Γe1 φ2−φ1 δx−qA+Γe2 φ3−φ2 δx−Γw2 φ2−φ1 δx+ Γe3 φ4−φ3 δx−Γw3 φ3−φ2 δx+qB−Γw4 φ4−φ3 δx=qB−qA (2.34)
2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 29 Teniendo en cuenta que Γe1=Γw2 , Γe2=Γw3 y Γe3=Γw4 , los flujos a través de las caras se cancelan dos a dos en el balance neto global, quedando únicamente los términos de flujo en los límites del dominio fluido ( qA y qB ). Esta consistencia en la formulación asegura la conservación de la propiedad φ cuando se usan las diferencias centradas. Puede ocurrir que una formulación con otro tipo de interpolación (por ejemplo, cuadrática) lleve a inconsistencias en la conservación. Los flujos en las caras pueden ser desiguales si los gradientes de las parábolas en dichas caras son distintos entre sí, de forma que los términos no se cancelarían y la conservación global no se satisfaría. No obstante, existen esquemas cuadráticos perfectamente consistentes y que se usan en la actualidad, como el QUICK scheme, aunque su estudio queda fuera del alcance de este proyecto. Delimitación La resolución del sistema de ecuaciones algebraico se lleva a cabo por métodos iterativos que persiguen la convergencia de la solución. En el marco de este proyecto, solo interesa saber que dichos métodos están embebidos en el comando \ de MATLAB. Puede demostrarse que es condición suficiente para la convergencia del método iterativo en cuestión que se cumpla el criterio de Scarborough. Éste puede expresarse de la siguiente forma según el valor de los coeficientes del sistema: ∑aW,E |aP|≤1en todos los nodos <1en al menos un nodo (2.35) En esta condición, aP representa el valor del coeficiente que multiplica al nodo central, y aW.E el del coeficiente que multiplica a todos los nodos de la vecindad (en unidimensional, éstos son únicamente el nodo Este y el nodo Oeste). Si dicha condición se verifica, la matriz de coeficientes resultante estará bien condicionada y será diagonal dominante, de forma que los métodos iterativos convergerán a la solución del sistema. Físicamente, esta condición establece que, en ausencia de fuentes, los valores de la variable φestán delimitados por sus valores en el contorno. Otra condición que debe cumplirse para garantizar una buena delimitación es que todos los coeficientes del sistema tengan el mismo signo (cuando se expresa de la forma aPφP=aEφE+aWφW ). Físicamente esto implica en que un aumento de φ en un nodo se traduce en un aumento de φ en los nodos vecinos. Si el esquema de discretización en cuestión no cumple alguna de estas dos condiciones, la solución divergirá o presentará oscilaciones, proporcionando resultados poco precisos. Transportividad Para analizar el concepto de transportividad de φ , resulta interesante definir el conocido como número de Peclet, un número adimensional que mide la importancia relativa de la convección frente a la difusión en un determinado dominio fluido. Así: Pe =F D=ρuA ΓA/δx(2.36) Donde δxrepresenta la longitud característica (ancho de la celda). A continuación se representan esquemáticamente las líneas de φ constante según diversos valores de Peclet. Para ilustrar este concepto, considérense los casos extremos: difusión pura ( Pe →0 ) y convección pura (Pe →∞).
30 Capítulo 2. FVM para problemas convectivo-difusivos en 1D W P E WPE Figura 2.13 Distribución de φ en la vecindad de dos fuentes. En la figura superior, el caso límite Pe →0. En la figura inferior, el caso límite Pe →∞. Para el caso en el que solo hay difusión ( Pe →0 ), el fluido estará remansado y las líneas de φ constante describirán circunferencias centradas en los nodos, ya que el proceso de difusión se extiende equitativamente de forma radial. Así, puede afirmarse que la influencia que ejercen los nodos E y W sobre P es la misma. A medida que el número de Pe se incrementa, las líneas de φ constante tienden a adoptar una forma elíptica y con centro arrastrado aguas abajo del nodo, tal y como puede observarse en la figura inferior. Así, cada nodo P se verá mayoritariamente influenciado por el nodo que se encuentre aguas arriba ( W o E según el sentido de la corriente). En el límite de Pe →∞ (sin difusión), toda la propiedad φ que emane de los nodos W o E se verá inmediatamente convectada aguas abajo. De esta forma, en el caso de un flujo en sentido positivo, el nodo P no se verá afectado por el nodo E , sino que estará completamente condicionado por el nodo W . Tanto es así que, en ausencia de difusión, puede declararse que φP=φW y φE=φP para un flujo positivo. La consideración conjunta de la direccionalidad del flujo, y por tanto de la influencia, y la magnitud del número de Peclet es lo que se conoce como transportividad, y será una característica de vital importancia en la caracterización de los esquemas de discretización. Una vez presentadas las propiedades que todo esquema de discretización debe tener para garantizar una convergencia adecuada, se está en condiciones de establecer una valoración objetiva de la precisión de cualquiera de los métodos. Comenzando por el método de las diferencias centradas, el único presentado hasta ahora, puede decirse que: a) Conservatividad. El método de las diferencias centradas emplea expresiones consistentes para evaluar los flujos en las caras y asegurar que el balance neto de flujo se cumple, por lo que esta propiedad se satisface. b) Delimitación. Si la expresión (2.31) se reordena de la forma aPφP=aWφW+aEφE y se tiene en cuenta que en el caso unidimensional ne=1 y nw=−1 , puede comprobarse que aP=aE+aW , por lo que el criterio de Scarborough se cumple. La segunda condición para garantizar la delimitación es que todos los coeficientes de la ecuación discretizada tengan el mismo signo. Teniendo en cuenta las consideraciones anteriores, el valor de los coeficientes es el que sigue:
2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 31 aW=Dw+Fw 2aP=Dw+Fw 2+De−Fe 2+(Fe−Fw)aE=De−Fe 2(2.37) Por ser un método conservativo, Fe−Fw=0 , luego el único coeficiente que corre riesgo de ser negativo será aE . Esto dependerá de cuán grande sea la contribución convectiva en dicho coeficiente. Como De>0yFe>0, para que aEsea positivo debe satisfacerse la siguiente condición: Fe De =ρeueAe Γe/δx=Pee<2(2.38) Así, valores de Pee mayores que 2 pueden llevar a soluciones que carezcan de sentido físico. Sujeto a los parámetros que caracterizan el número de Peclet de la celda, puede afirmarse que la condición (2.38) se cumple solo si se trata de un flujo difusivo-dominante a bajo número de Re (velocidad u pequeña) o si el mallado es suficientemente fino ( δx pequeño). Esto último no siempre será posible, lo que limita la aplicación de este método a un número reducido de problemas. c) Transportividad. En la formulación del método de las diferencias centradas la influencia que ejercen sobre P los flujos convectivo-difusivos es la misma para todos los nodos vecinos. De esta forma, no se está considerando la dirección del flujo ni la importancia relativa de la convección frente a la difusión. Puede concluirse por tanto que la propiedad de transportividad no es una característica del método de las diferencias centradas para altos Pe. En compendio, este análisis pone de manifiesto que las propiedades de delimitación y transportividad pueden no cumplirse para números de Pe lo suficientemente elevados en el método de las diferencias centradas. A pesar de ser un método de segundo orden, no es un esquema de discretización favorable para numerosos problemas convectivo-difusivos de carácter general. Así, surge la necesidad de desarrollar nuevos métodos que permitan superar estas limitaciones. A continuación se presentan algunos de los más conocidos, y que están implementados en softwares comerciales de extendido uso como OpenFOAM o Ansys Fluent (véase [7]). 2.2.3 Esquema upwind de primer orden El esquema upwind o FOU (acrónimo de su nombre en inglés, First Order Upwind) es un método de discretización que nace como solución a una de las mayores incompetencias del método de las diferencias centradas: su incapacidad de identificar la dirección del flujo. En dicho método el valor de la propiedad φ en la cara de una celda está igualmente influenciado por el valor en los nodos colindantes, independientemente de las características del flujo. Por contra, el método upwind establece el valor convectivo de φ en las caras igual al valor de φ en el nodo aguas arriba. Así, se rechaza la influencia equitativa de ambos nodos y se tiene en cuenta la influencia dominante del nodo aguas arriba de la corriente. Esta forma de proceder es precisamente la que da nombre al método. Matemáticamente, este modelo se expresa de la siguiente manera (destacar que para las condiciones en los flujos se considera Fw=ρwuwAwyFe=ρeueAe, sin tener en cuenta las normales): φw=φW,φe=φPsi Fw>0,Fe>0 φw=φP,φe=φEsi Fw<0,Fe<0(2.39) Esto lleva a la siguiente ecuación discretizada (para un flujo positivo): (Fw+Dw)φW+(Fe+De−Dw)φP+ (−De)φE=¯ S∆V(2.40)
38 Capítulo 2. FVM para problemas convectivo-difusivos en 1D φw=φW−(1−0.1Pew)5 Pew (φP−φW)si 0<Pew<10 φw=φWsi Pew>10 φe=φP−(1−0.1Pee)5 Pee (φE−φP)si 0<Pee<10 φe=φPsi Pee>10 (2.50) -20 -15 -10 -5 0 5 10 15 20 -10 -5 0 5 10 15 20 25 30 Figura 2.20 Comparación entre la solución exacta y la aproximación polinómica propuesta por Patankar. Si se aplica este método al problema (2.32) puede observarse que los errores obtenidos son con diferencia los más reducidos frente al resto de métodos. En efecto, obsérvese la siguiente figura: 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2.21 Comparación entre la solución exacta y la obtenida mediante el FVM para u=3m/s empleando 20 nodos equiespaciados y usando el esquema de ley de potencia de Patankar.
2.2 Problema convectivo-difusivo, estacionario y unidimensional. Esquemas de discretización 39 En conclusión, las propiedades del esquema de ley de potencia lo hacen sobresalir frente a los demás, pues pretende representar la solución exacta a un coste computacional reducido. Es una alternativa de extendido uso y códigos comerciales como Ansys Fluent lo implementan como opción de discretización para el usuario (véase [7]). En último lugar, resulta bastante representativo graficar en una misma figura todos los métodos empleados para varios valores del número de nodos y número de Peclet. Los resultados son congruentes con lo que se ha ido comentando a lo largo del capítulo, quedando claramente representadas las limitaciones de cada método desarrollado y la mejora que progresivamente se ha ido consiguiendo con cada uno de ellos. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2.22 Comparación entre la solución exacta y los distintos métodos de discretización para u=0.1m/sempleando 5 nodos equiespaciados. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.5 0 0.5 1 1.5 2 2.5 3 Figura 2.23 Comparación entre la solución exacta y los distintos métodos de discretización para u=3m/sempleando 5 nodos equiespaciados.
40 Capítulo 2. FVM para problemas convectivo-difusivos en 1D 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 2.24 Comparación entre la solución exacta y los distintos métodos de discretización para u=3m/sempleando 20 nodos equiespaciados. En el próximo capítulo se dará un paso más en el estudio de los problemas convectivo-difusivos relajando la hipótesis de estacionariedad. Se expondrán y analizarán varios métodos para tratar numéricamente el término de variación temporal de la ecuación general, así como ventajas y limitaciones de los mismos, y se aplicarán a un problema característico. 2.3 Problema difusivo, no estacionario y unidimensional Una vez estudiados en detalle los diferentes métodos para resolver analíticamente problemas convectivo-difusivos en régimen estacionario, se está en condiciones de abordar el régimen transitorio o no estacionario. Para ello y siguiendo la exposición de [ 4 ], se recuperará el término de variación temporal de la ecuación general (1.5), que representa la tasa de cambio de la propiedad φ con respecto del tiempo. Así, retener la variable temporal requiere una nueva integración respecto al tiempo en la ecuación (1.6), resultando la siguiente expresión general: ZCV Zt+∆t t ∂ ∂t(ρφ)dt dV+Zt+∆t tZA n·(ρφv)dAdt = Zt+∆t tZA n·(Γ∇φ)dAdt +Zt+∆t tZCV SφdV dt (2.51) La integración a lo largo del volumen de control es esencialmente idéntica a lo que se ha desarrollado en el caso estacionario, así como los métodos expuestos para el tratamiento de los términos convectivo-difusivos. Por ello, en este capítulo se centrarán los esfuerzos en la integración temporal a lo largo de un paso ∆t . La forma más sencilla de abordar esta tarea es partiendo nuevamente de un problema paradigmático que permita focalizar la atención en lo que verdaderamente se está interesado. Es buen ejemplo la transferencia de calor no estacionaria en una barra unidimensional. Este es un problema puramente difusivo y está gobernado por la siguiente ecuación: ρc∂T ∂t=∂ ∂xk∂T ∂x+Sφ(2.52)
2.3 Problema difusivo, no estacionario y unidimensional 41 Donde c es el calor específico del fluido y se mide en J/(KgK) . Integrando con respecto al tiempo y al volumen de control, se tiene: Zt+∆t tZCV ρc∂T ∂tdV dt =Zt+∆t tZCV ∂ ∂xk∂T ∂xdV dt +Zt+∆t tZVC SdV dt (2.53) Si se considera la temperatura en los nodos como valor representativo del volumen de control, el primer miembro de la ecuación (2.53) puede aproximarse por: ZCV Zt+∆t t ρc∂T ∂tdtdV =ρc(TP−T0 P)∆V(2.54) El valor T0 P hace referencia a la temperatura en el instante t , mientras que TP se refiere al valor de la temperatura en el instante t+∆t . En adelante, se mantendrá esta notación para distinguir los diferentes instantes temporales. Nótese que para llegar a este resultado se ha discretizado el término ∂T/∂t usando la diferenciación central ( ∂T/∂t= (TP−T0 P)/∆t ). De igual forma, si se aplica el método de las diferencias centradas para modelar los términos difusivos tal y como se ha venido haciendo a lo largo de la memoria, se llega a la siguiente expresión: ρc(TP−T0 P)∆V=Zt+∆t tkeAene TE−TP xE−xP+kwAwnw TP−TW xP−xWdt +Zt+∆t t ¯ S∆V dt (2.55) Para evaluar las integrales del segundo miembro de la ecuación (2.55) debe modelarse la variación de la temperatura con respecto al tiempo. Puede hacerse la suposición de que el valor que prevalece es el que toma en el instante t , aunque también podría suponerse que el valor representativo es el del instante t+∆t . Para resolver esta problemática sin pérdida de generalidad, lo más útil es aproximar el valor de la integral mediante una ponderación de ambos valores, caracterizada por el parámetro θ∈[0,1]. Así: IT=Zt+∆t t TPdt = [θTP+(1−θ)T0 P]∆t(2.56) Aplicando la discretización de (2.56) en (2.55) y teniendo en cuenta que De=keAene/(xE−xP) yDw=kwAwnw/(xP−xW), se obtiene: −θDeTE+ [θ(De−Dw)+ ρc∆V ∆tTP+θDwTW= (1−θ)DeT0 E+ (1−θ)(Dw−De)+ ρc∆V ∆tT0 P−(1−θ)DwT0 W+¯ S∆V (2.57) Que puede escribirse en la forma común: aETE+aPTP+aWTW=a0 ET0 E+a0 PT0 P+a0 WT0 W+¯ S∆V=aSP (2.58) La forma final de la ecuación discretizada dependerá del valor que se asigne a θ . Para θ=0 , IT=T0 P∆t , por lo que en la expresión únicamente se usarán las temperaturas en el instante anterior t para obtener el valor de TP . El esquema de discretización resultante se denomina explícito. Cuando 0<θ≤1 , aparecen las temperaturas en el instante nuevo t+∆t en ambos miembros de la ecuación, y los esquemas resultantes se conocen como implícitos. El caso límite de θ=1 , donde IT=TP∆t , se designa como esquema completamente implícito. El caso correspondiente a θ=1/2 se conoce como esquema de Crank-Nicolson, y se verifica que IT=1/2(TP+T0 P)∆t.
42 Capítulo 2. FVM para problemas convectivo-difusivos en 1D En los capítulos venideros se desarrollarán minuciosamente los esquemas explícito y completamente implícito y se ilustrará su aplicación tomando como ejemplo el problema descrito en la Figura 2.25 y gobernado por (2.52). El esquema de Crank-Nicolson queda fuera del alcance de este documento por ser su uso menos extendido que los casos límite de θ=0yθ=1. 2.3.1 Esquema explícito Para analizar y valorar en profundidad las propiedades de los esquemas de discretización explícito y completamente implícito y poder establecer una comparación entre ambos, considérese el siguiente problema descrito en [ 4 ]. Sea una placa metálica y delgada de espesor L=2cm que se encuentra inicialmente a una temperatura constante de T=200 ◦C . En un determinado instante ( t=0 ), la temperatura de la cara este se reduce instantáneamente a TB=0◦C . La cara oeste permanece aislada térmicamente ( ∂T/∂x|A=0 ). Además, la placa presenta una conductividad térmica k=10W/(mK) y un factor ρc=10 ·106J/(m3K) . Puede demostrarse que la solución analítica de este problema es la que se muestra en (2.59). 12345 Figura 2.25 Geometría del problema que ejemplifica el caso no estacionario para t>0. T(x,t) T0 =4 π ∞ ∑ n=1 (−1)n+1 2n−1e−αλ2 ntcos(λnx)(2.59) Donde λn=(2n−1)π 2L y α=k/ρc . Bajo estas condiciones, el esquema explícito lleva a la siguiente expresión: ρc∆V ∆tTP=DeT0 E+Dw−De+ρc∆V ∆tT0 P−DwT0 W(2.60) Tal y como se explicó con detalle en el Capítulo 2.2.2, para que la solución tenga una convergencia adecuada todos los coeficientes del sistema deben tener el mismo signo. El coeficiente a0 P puede entenderse como un coeficiente del nodo vecino T0 P (en este caso un vecino temporal), pues conecta los valores en P del instante t con los del instante t+∆t . En tal caso, reescribiendo (2.60) como aPTP=a0 ET0 E+a0 PT0 P+a0 WT0 W , considerando un mallado uniforme ( xE−xP=xP−xW=∆x ), área constante y teniendo en cuenta que ne=1 y nw=−1 en el caso unidimensional, deberá cumplirse que:
2.3 Problema difusivo, no estacionario y unidimensional 43 a0 P>0⇒ρcA∆x ∆t−kA ∆x−kA ∆x>0⇒∆t<ρc∆x2 2k(2.61) Esta condición establece una cota superior en el paso de integración y representa una limitación considerable para el esquema explícito. Un paso que no verifique la desigualdad puede llevar a soluciones completamente alejadas de la solución real, mientras que escoger un paso dentro del intervalo impuesto puede ser computacionalmente muy costoso. En particular, el paso se reduce como ∆x2 , por lo que el refinamiento del mallado necesario puede llegar a ser inabordable. Es por todo esto que el método explícito no suele ser una buena recomendación para resolver problemas transitorios de carácter general. Todo esto se ilustra a continuación: 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 20 40 60 80 100 120 140 160 180 200 Solución exacta (t=40) Solución con FVM (t=40) Solución exacta (t=80) Solución con FVM (t=80) Solución exacta (t=120) Solución con FVM (t=120) Figura 2.26 Comparación entre la solución analítica y la solución mediante el FVM empleando 5 nodos, usando el esquema explícito para el término no estacionario y con un paso ∆t=2s. Tabla 2.4 Valores del campo de temperaturas y error relativo cometido por el FVM usando el esquema explícito, para ∆t=2sy empleando 5 nodos equiespaciados en t=40 s. Nodo 1 2 3 4 5 Solución FVM (◦C) 188.6386 176.4132 148.2926 100.7597 35.9418 Solución Exacta (◦C) 188.3845 175.7649 147.1303 99.5043 35.3836 Error Relativo (%) 0.1349 0.3688 0.7900 1.2616 1.5777 Nótese cómo a medida que transcurre el tiempo la temperatura disminuye, como era de esperar. Para el caso particular del problema en cuestión y con una discretización de 5 nodos, se tiene que la condición a cumplir es ∆t<8s . Los resultados obtenidos en la Figura 2.26 son extraordinariamente precisos debido a que ∆t=2s. Escogiendo un paso fuera del intervalo:
44 Capítulo 2. FVM para problemas convectivo-difusivos en 1D 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 -100 -50 0 50 100 150 200 250 Solución exacta (t=40) Solución con FVM (t=40) Figura 2.27 Comparación entre la solución analítica y la solución mediante el FVM empleando 5 nodos equiespaciados, usando el esquema explícito para el término no estacionario y con un paso ∆t=10 sen t=40 s. En efecto, con un paso ∆t=10 s la solución oscila violentamente en torno a la solución analítica y claramente diverge de ésta. El problema de este método radica en que para mejorar la precisión de la solución es necesario aumentar el número de nodos, lo que disminuye el paso cuadráticamente. La exactitud de la Figura 2.26 se debe a la simplicidad del problema descrito, pero imagínese el escenario en que la solución mostrada en la Tabla 2.4 no es válida porque el máximo error relativo permitido sea del 0.1 % . En tal caso, el único camino sería aumentar el número de nodos. Por ejemplo, usando 20 nodos: 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 -2 -1.5 -1 -0.5 0 0.5 1 1.5 21018 Solución exacta (t=40) Solución con FVM (t=40) Figura 2.28 Comparación entre la solución analítica y la solución mediante el FVM empleando 20 nodos equiespaciados, usando el esquema explícito para el término no estacionario y con un paso ∆t=2sen t=40 s.
2.3 Problema difusivo, no estacionario y unidimensional 45 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 20 40 60 80 100 120 140 160 180 200 Solución exacta (t=40) Solución con FVM (t=40) Figura 2.29 Comparación entre la solución analítica y la solución mediante el FVM empleando 20 nodos equiespaciados, usando el esquema explícito para el término no estacionario y con un paso ∆t=0.25 sen t=40 s. La Figura 2.28 ilustra a la perfección la limitación del esquema explícito. Con 20 nodos y un mallado homogéneo debe cumplirse que ∆t<0.5s y, efectivamente, con ∆t=2s los resultados obtenidos son completamente desorbitados y no se acercan para nada a la solución. Bajando el paso a ∆t=0.25 s , la solución numérica se adapta perfectamente a la analítica, logrando un error relativo máximo menor al 0.1 % tal y como se muestra en la Figura 2.29. 2.3.2 Esquema completamente implícito Para solucionar la problemática computacional que presenta el esquema explícito, se desarrollan a continuación las expresiones sobre las que se sustenta el esquema completamente implícito aplicado al problema descrito. Haciendo θ=1en la ecuación (2.57), se llega a la siguiente expresión: −DeTE+De−Dw+ρc∆V ∆tTP+DwTW=ρc∆V ∆tT0 P(2.62) Si ahora se reescribe la ecuación (2.62) como aPTP=aETE+aPTP+aWTW+a0 PT0 P , considerando un mallado uniforme ( xE−xP=xP−xW=∆x ), área constante y teniendo en cuenta que ne=1 y nw=−1en el caso unidimensional, se tiene: De+Dw+ρc∆V ∆tTP=DwTW+DeTE+ρc∆V ∆tT0 P(2.63) Puede comprobarse que este esquema tiene la particularidad de que todos los coeficientes del sistema son positivos independientemente de su valor, lo que convierte a esta práctica en un método incondicionalmente estable para cualquier paso. Son su robustez y estabilidad los que hacen que este método sea el recomendado para cálculos transitorios de carácter general. A continuación se ilustra este esquema mediante su aplicación al problema descrito en (2.59).
46 Capítulo 2. FVM para problemas convectivo-difusivos en 1D 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 20 40 60 80 100 120 140 160 180 200 Solución exacta (t=40) Solución con FVM (t=40) Solución exacta (t=80) Solución con FVM (t=80) Solución exacta (t=120) Solución con FVM (t=120) Figura 2.30 Comparación entre la solución analítica y la solución mediante el FVM empleando 5 nodos equiespaciados, usando el esquema implícito para el término no estacionario y con un paso ∆t=2s. Para el mismo paso y número de nodos, puede comprobarse como la solución, aunque ligeramente menos precisa que en el caso explícito, se aproxima con bastante exactitud a la solución analítica. Si se quiere mejorar la precisión puede aumentarse el número de nodos, esta vez sin necesidad de disminuir el paso. 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 20 40 60 80 100 120 140 160 180 200 Solución exacta (t=40) Solución con FVM (t=40) Solución exacta (t=80) Solución con FVM (t=80) Solución exacta (t=120) Solución con FVM (t=120) Figura 2.31 Comparación entre la solución analítica y la solución mediante el FVM empleando 20 nodos equiespaciados, usando el esquema implícito para el término no estacionario y con un paso ∆t=2s.
2.3 Problema difusivo, no estacionario y unidimensional 47 Llegados a este punto se han desarrollado todos y cada uno de los esquemas de discretización que se emplearán a lo largo de la memoria para modelar los términos difusivos, convectivos y transitorios. Todos ellos, salvando las distancias con ciertas variaciones, son extrapolables al caso bidimensional. Su estudio será abordado en el próximo capítulo.
54 Capítulo 3. FVM para problemas convectivo-difusivos en 2D. Mallado estructurado [Dn−Dw−Ds]φP+DwφW−DnφN+DsφS=qEAe(3.17) En el caso en que se impusiese una temperatura constante de valor TE=φright , se tendría que ∂T ∂xw=φright −TP ∆x/2y se llegaría a la siguiente expresión: [2De+Dn−Dw−Ds]φP+DwφW−DnφN+DsφS=2Deφright (3.18) Cara norte (k=14,k=15) En la cara norte se impone un flujo saliente de calor por unidad de superficie de valor qN= −150 kW /m2 . El área que atraviesa dicho flujo será An=t∆x para ambos volúmenes de control. Por lo tanto, en la ecuación (3.11) habrá que imponer knAnnn∂T ∂yn=qNAn. Así: [De−Dw−Ds]φP−DeφE+DwφW+DsφS=qNAn(3.19) Por último, en la situación en la que se impusiese una temperatura constante de valor TN=φup , se tendría que ∂T ∂yn=φup−TP ∆y/2y se llegaría a la siguiente expresión: [De+2Dn−Dw−Ds]φP−DeφE+DwφW+DsφS=2Dnφup (3.20) Una vez caracterizadas todas las caras limítrofes, se pasa a modelar las celdas ubicadas en las esquinas del dominio computacional. En ellas deberán verificarse simultáneamente dos condiciones de contorno. Esquina sudoeste (k=1) En la esquina sudoeste deben imponerse las condiciones ∂T ∂ys=TP−TS ∆y/2 y kwAwnw∂T ∂xw=0 . La ecuación resultante será: [De+Dn−2Ds]φP−DeφE−DnφN=−2DsTS(3.21) Esquina sudeste (k=4) En la esquina sudeste deben imponerse las condiciones ∂T ∂ys=TP−TS ∆y/2 y keAene∂T ∂xe=qEAe . La ecuación resultante será: [Dn−Dw−2Ds]φP+DwφW−DnφN=−2DsTS+qEAe(3.22) Esquina noreste (k=16) En la esquina noreste deben imponerse las condiciones knAnnn∂T ∂yn=qNAn y keAene∂T ∂xe=qEAe . La ecuación resultante será: [−Dw−Ds]φP+DwφW+DsφS=qEAe+qNAn(3.23) Esquina noroeste (k=13) En la esquina noroeste deben imponerse las condiciones knAnnn∂T ∂yn=qNAn y kwAwnw∂T ∂xw=0 . La ecuación resultante será: [De−Ds]φP−DeφE+DsφS=qNAn(3.24)
3.1 Problema difusivo, estacionario y bidimensional 55 Una vez obtenido el conjunto de ecuaciones discretizadas que deben satisfacerse en cada uno de los nodos de la malla, se resuelve el sistema lineal mediante el comando \ de MATLAB. A continuación, se presenta la solución calculada. 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Figura 3.3 Mallado del dominio de computación representado en MATLAB. Los puntos rojos son los nodos donde se calcularán los valores de temperatura. Figura 3.4 Valores del campo de temperaturas en los nodos para un mallado uniforme M×N con M=4yN=4.
56 Capítulo 3. FVM para problemas convectivo-difusivos en 2D. Mallado estructurado Figura 3.5 Valores del campo de temperaturas en los nodos para un mallado uniforme M×N con M=50 yN=50. Son varios los comentarios que pueden hacerse acerca de las representaciones anteriores. En primer lugar, el problema que se ha resuelto mediante el FVM no cuenta con una solución analítica, por lo que los resultados no pueden ser contrastados contra una solución exacta. No obstante, éstos han sido contrastados con los que se obtienen en la literatura y han demostrado ser perfectamente válidos. En cualquier caso, para dar mayor credibilidad a los resultados, puede realizarse una valoración cualitativa de lo obtenido. La solución es congruente con lo que cabría esperar por varias razones: primeramente, cabía esperar que el valor máximo de temperatura se obtuviese cerca de la mitad de la cara este, pues existe un flujo de calor entrante por dicha cara, mientras que en las caras adyacentes se tiene un flujo de calor saliente (cara norte) y una temperatura impuesta (cara sur) menor a la temperatura que produce qE . Por otro lado, la región más fría de la placa es la suroeste, pues la cara sur se mantiene a una temperatura constante de 150◦C y la cara oeste está aislada. La
3.2 Problema convectivo-difusivo, estacionario y bidimensional 57 temperatura aumenta gradualmente a medida que se avanza hacia la derecha en la placa, como es lógico. Además, las caras aisladas térmicamente se caracterizan porque las curvas isotermas son perpendiculares a ellas. Esto se aprecia claramente en la siguiente figura. 160 160 160 160 170 170 170 170 180 180 180 190 190 190 200 200 210 210 220 220 230 230 240 240 250 260 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Figura 3.6 Curvas de nivel (isotermas) para un mallado uniforme M×Ncon M=50 yN=50. En segundo lugar, resulta interesante comentar que con un mayor número de nodos se obtiene una representación mucho más realista debido a que la ubicación de los mismos se aproxima a los límites, logrando que la implementación de las condiciones de contorno sea más fidedigna. En la siguiente sección, se resolverá y analizará un problema bidimensional convectivo-difusivo estacionario. Se recuperarán en las ecuaciones los términos de convección y se aplicará alguno de los modelos estudiados en el Capítulo 2 para su tratamiento. 3.2 Problema convectivo-difusivo, estacionario y bidimensional Para continuar con el análisis de problemas convectivo-difusivos en 2D se procede a resolver un problema típico que incluye ambos fenómenos de transporte, similar a uno de los problemas que se encuentran en [ 6 ]. El tratamiento del término referente al flujo convectivo permanece inalterado en problemas multidimensionales, siempre y cuando se haga un uso correcto de las componentes de la velocidad, que ahora puede no estar alineada con los ejes coordenados (dando lugar al fenómeno de falsa difusión). Sirva de ejemplo el problema que se considera a continuación, una situación típica dentro del campo de la transferencia de calor. Para simplificar, el problema se resolverá haciendo uso únicamente del esquema upwind de primer orden para el término advectivo. Considérese un conducto plano bidimensional por el que circula agua, de dimensiones L=0.5m y H=0.1m . El agua entra en el canal por la cara izquierda, que se encuentra a una temperatura constante de TW=30 ◦C . Una resistencia calienta el conducto por la cara superior con un flujo de calor por unidad de superficie de valor constante qN=−k∂T ∂yN=15W/m2 . Por otro lado, la cara inferior está completamente aislada (qS=−k∂T ∂yS=0) . En la cara de salida se tiene un flujo de calor por unidad de superficie negativo (sale del dominio) de valor (qE=−k∂T ∂xE=−20W/m2) . Bajo estas condiciones, las propiedades termodinámicas del agua presentan los siguientes valores constantes: ρ=1000 kg/m3 , cp=4200 J/kgK y k=0.613W/mK . También se asumirá que el flujo
58 Capítulo 3. FVM para problemas convectivo-difusivos en 2D. Mallado estructurado es laminar y está completamente desarrollado. En este escenario, las componentes de la velocidad presentan las siguientes expresiones: u=u(y) = 6umy H−y H2v=0(3.25) Donde um es la velocidad media en dirección x a la entrada del canal. Su valor será tal que el Pe =ρumH/Γ=200 , siendo Γ=k/ρcp . Así, se obtiene que um=2.92 ·10−7m/s El objetivo es obtener la distribución de temperaturas dentro del conducto. Esta situación se describe en la Figura 3.7. La ecuación que gobierna este problema de convección-difusión estacionaria es la siguiente: ∂ ∂x(ρucpT) + ∂ ∂y(ρvcpT) = ∂ ∂xk∂T ∂x+∂ ∂yk∂T ∂y(3.26) Figura 3.7 Representación del problema descrito. Para una malla uniforme de tamaño genérico M×N se tiene que en cada volumen de control debe satisfacerse la siguiente expresión: ρcpueneAeφe+ρcpuwnwAwφw+ρcpvnnnAnφn+ρcpvsnsAsφs= keAene ∂T ∂xe +kwAwnw ∂T ∂xw +knAnnn ∂T ∂yn +ksAsns ∂T ∂ys (3.27) Teniendo en cuenta que para este problema en particular se verifica que kw=ke=kn=ks=k , Aw=Ae=∆y , As=An=∆x , ne=nn=1 , nw=ns=−1 , u=u(y) y v=0 , pueden definirse los siguientes términos difusivos: Dw=−k∆y/∆x , Ds=−k∆x/∆y , De=k∆y/∆x y Dn=k∆x/∆y . Por otro lado, también pueden definirse los siguientes términos convectivos: Fw=−ρcpu(y)∆y , Fe=ρcpu(y)∆y , Fn=0 y Fs=0 . Así, usando el método de las diferencias centradas para aproximar los gradientes, se llega a la siguiente expresión: [Fe+De+Dn−Dw−Ds]φP−DeφE+(Fw+Dw)φW−DnφN+DsφS=0(3.28) Esta ecuación discretizada servirá únicamente para los nodos intermedios ( i=2 : M−1 , j= 2 : N−1 ). El resto de celdas estará sujeta a las condiciones de contorno que presenta el problema descrito. Se exponen a continuación de manera sucinta las expresiones que deben satisfacerse en cada volumen de control fronterizo y en las celdas ubicadas en las esquinas.
3.2 Problema convectivo-difusivo, estacionario y bidimensional 59 Cara oeste En la cara oeste se impone la condición de temperatura conocida TW=30 ◦C . Debe verificarse que ∂T ∂xw=TP−TW ∆x/2. La ecuación resultante será: [Fe+De+Dn−2Dw−Ds]φP−DeφE−DnφN+DsφS=−[2Dw+Fw]TW(3.29) Cara norte En la cara norte se impone la condición de flujo conocido, de valor qN=15W/m2 . El área que atraviesa dicho flujo es An=∆x . Debe verificarse que knAnnn∂T ∂yn=qNAn . La ecuación resultante será: [Fe+De+Dn−Dw]φP−DeφE+[Fw+Dw]φW−DnφN=qNAn(3.30) Cara este En la cara este se impone un flujo saliente de valor qE=−20W/m2 . El área que atraviesa dicho flujo es Ae=∆y . Como consecuencia habrá que imponer keAene∂T ∂xe=qEAe . La ecuación resultante será: [Fe+Dn−Dw−Ds]φP+[Fw+Dw]φW−DnφN+DsφS=qEAe(3.31) Cara sur La cara sur se encuentra aislada térmicamente. Como consecuencia habrá que imponer ksAsns∂T ∂ys= 0. La ecuación resultante será: [Fe+De+Dn−Dw]φP−DeφE+[Fw+Dw]φW−DnφN=0(3.32) Esquina sudoeste En la esquina sudoeste deben imponerse las condiciones ksAsns∂T ∂ys=0 y ∂T ∂xw=TP−TW ∆x/2 . La ecuación resultante será: [Fe+De+Dn−2Dw]φP−DeφE−DnφN=−[2Dw+Fw]TW(3.33) Esquina sudeste En la esquina sudeste deben imponerse las condiciones ksAsns∂T ∂ys=0 y keAene∂T ∂xe=qEAe . La ecuación resultante será: [Fe+Dn−Dw]φP+[Fw+Dw]φW−DnφN=qEAe(3.34) Esquina noreste En la esquina noreste deben imponerse las condiciones knAnnn∂T ∂yn=qNAn y keAene∂T ∂xe=qEAe . La ecuación resultante será: [Fe−Dw−Ds]φP+[Fw+Dw]φW+DsφS=qNAn+qEAe(3.35) Esquina noroeste En la esquina noroeste deben imponerse las condiciones knAnnn∂T ∂yn=qNAn y ∂T ∂xw=TP−TW ∆x/2 . La ecuación resultante será:
60 Capítulo 3. FVM para problemas convectivo-difusivos en 2D. Mallado estructurado [Fe+De−2Dw−Ds]φP−DeφE+DsφS=−[2Dw+Fw]TW+qNAn(3.36) Una vez obtenido el conjunto de ecuaciones discretizadas que deben satisfacerse en cada uno de los nodos de la malla, se resuelve el sistema lineal mediante el comando \ de MATLAB. A continuación, se presenta la solución calculada. Figura 3.8 Valores del campo de temperaturas en los nodos para un mallado M×N con M=50 y N=10.
3.2 Problema convectivo-difusivo, estacionario y bidimensional 61 31 31 32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 40 40 41 41 42 42 42 43 43 43 43 44 44 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Figura 3.9 Curvas de nivel (isotermas) para un mallado M×Ncon M=50 yN=10. Puede observarse cómo los resultados obtenidos son congruentes con lo esperado. En efecto, como el canal se calienta por la parte superior, la mayor temperatura se alcanza en dicha zona. La convección de izquierda a derecha se manifiesta con claridad en la curvatura que presentan las isotermas representadas en la Figura 3.9. Por último, la fronteras aislada térmicamente se identifica fácilmente al ser las curvas de nivel ortogonales a dicha pared. Una vez analizado el caso bidimensional con mallas estructuradas para problemas de conveccióndifusión térmica, se procede en el próximo capítulo a desarrollar el FVM para geometrías 2D con mallas no estructuradas. El estudio del término no estacionario se ha omitido para el caso bidimensional por ser su implementación y los modelos en los que se sustenta esencialmente idénticos a los presentados en la Sección 2.3.
4 FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado Hasta ahora las técnicas de resolución para las ecuaciones de flujo estaban basadas en esquemas de discretización sustentados en malladas estructuradas cartesianas. Este escenario es el más simple de todos, y ha sido particularmente útil para introducir y desarrollar los fundamentos del método de los volúmenes finitos de una forma fácil de entender. No obstante, la mayor parte de problemas que aparecen en ingeniería presentan complejas geometrías que difícilmente encajan con este tipo de mallas. Una posible forma de abordar situaciones en las que las fronteras del dominio fluido no coinciden con el sistema de coordenadas de la malla sería aproximar la geometría, por ejemplo usando una aproximación escalonada, tal y como se muestra en la siguiente figura extraída de [ 4 ]. Figura 4.1 Malla ortogonal cartesiana para aproximar un problema de flujo en un semicilindro. Sin embargo, una aproximación así es laboriosa e ineficiente. Las celdas que quedan fuera del dominio objeto de estudio contienen una información innecesaria, lo que se traduce en una pérdida de almacenamiento y un derroche de recursos, computacionalmente hablando. Además, la aproximación introduce errores al método de resolución que solo pueden ser salvados minimizando el tamaño de la malla, lo que implica un mayor refinamiento también en las zonas que no son de interés y por tanto una mayor ineficiencia. Esto pone de manifiesto la clara limitación que tienen este tipo de mallas y la necesidad de desarrollar nuevas estructuras que puedan lidiar con complejas geometrías de forma natural. 63
70 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado 22 v1 = Lfv(iface,1); 23 v2 = Lfv(iface,2); 24 xf(iface) = 0.5*(xv(v1) + xv(v2)); 25 yf(iface) = 0.5*(yv(v1) + yv(v2)); 26 end Son varios los comentarios que pueden hacerse en relación con el almacenamiento de la conectividad que se ha expuesto en las líneas de código anteriores. En primer lugar, uno podría preguntarse sobre la necesidad de este procesamiento. La respuesta a dicha pregunta se ejemplificará más adelante, pero fundamentalmente la exigencia de esto radica en la propia naturaleza del método. En los cálculos posteriores se requerirá la evaluación de ciertas magnitudes en nodos, caras y vértices de forma continuada, y la única forma eficiente de realizar esta tarea es disponer de un almacenamiento compacto de números enteros (ocupan poco espacio en la memoria) que sea fácilmente accesible. En segundo lugar, si uno computa los códigos adjuntos comprobará que la asignación local de caras y vértices se ha realizado siguiendo siempre un sentido antihorario. Esto se hace con la intención de generar matrices claras y fácilmente depurables en caso de encontrarse fallos. Por último, en la matriz de conectividad Cara-A-Celda, aunque la asignación k=1 ó k=2 se ha hecho de forma completamente arbitraria, se ha establecido por norma que aquellas caras limítrofes (que solo tienen una celda colindante) presenten el mismo índice global para cada índice local; es decir, si fb es una cara frontera, entonces L f c(fb,1) = L f c(fb,2) . El interés de esto se justificará posteriormente. A modo ilustrativo y tomando como referencia la malla representada en la Figura 4.2, puede decirse que Lc f (9,1) = 6 , Lc f (9,2) = 20 y Lc f (9,3) = 30 para ejemplificar la conectividad Celda-A-Cara, L f c(32,1) = 14 y Lc f (32,2) = 17 para la conectividad Cara-A-Celda, Lcv(7,1) = 5 , Lcv(7,2) = 6 y Lcv(7,3) = 10 para la conectividad Celda-A-Vértice y L f v(20,1) = 8 y L f v(20,2) = 12 para la conectividad Cara-A-Vértice. Todo esto puede comprobarse ejecutando secuencialmente los códigos. Prosiguiendo con el procesamiento de la geometría, se procede a continuación a desarrollar los algoritmos que permiten la obtención de los vectores tangentes y normales a cada cara, así como el área de las mismas (que en 2D no es más que la longitud del segmento). Particular interés tiene el cálculo de los vectores normales, pues debe asegurarse que éstos apuntan en la dirección exterior con respecto a cada celda para que el Teorema de la Divergencia de Gauss aplicado a la expresión (4.1) sea válido. El cálculo de las áreas es trivial una vez obtenida la matriz de conectividad que relaciona las caras y los vértices. Matemáticamente, puede expresarse como sigue: Af=q(x2−x1)2+(y2−y1)2(4.4) Donde los subíndices 1 y 2 hacen referencia a cada vértice de la cara. Una vez obtenida la longitud del segmento, el cálculo del vector tangente unitario es inmediato: tx,f=x2−x1 p(x2−x1)2+(y2−y1)2=x2−x1 Af ty,f=y2−y1 p(x2−x1)2+(y2−y1)2=y2−y1 Af (4.5) Las componentes del vector normal unitario pueden obtenerse a partir de éste resolviendo el siguiente sistema de ecuaciones: n· t=0 n× t= k⇒ nx,ftx,f+ny,fty,f=0 nx,fty,f+ny,ftx,f=1⇒ nx,f=ty,f ny,f=−tx,f (4.6)
4.1 Problema difusivo, estacionario y unidimensional 71 Una vez determinadas las componentes del vector normal unitario, es necesario desarrollar un método que discrimine la normal exterior en detrimento de la interior. Una posible solución es imponer que el vector normal calculado verifique la condición lc f ·n>0 , donde lc f es el vector que une el nodo de la celda con el centro de la cara en cuestión, tal y como se ilustra en la Figura 4.3. Si dicho producto escalar es positivo, entonces el vector normal es exterior; si es negativo, bastaría con cambiar de signo las componentes de ambos vectores tangente y normal, de forma que la expresión n× t= k se siga cumpliendo. A continuación se expone el código que realiza los cálculos descritos. 1 2 b a f Figura 4.3 Representación de los vectores unitarios tangente y normal en una cara localmente con respecto a una celda. 1%% CALCULO DE AREAS Y NORMALES Y TANGENTES UNITARIAS 2 3lfcx = zeros(n_c,3); 4lfcy = zeros(n_c,3); 5area = zeros(1,n_f); 6tx=zeros(n_c,3); 7ty=zeros(n_c,3); 8nx=zeros(n_c,3); 9ny=zeros(n_c,3); 10 for ic = 1:n_c 11 for k = 1:3 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);
72 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado 30 end 31 end 32 end En la Figura 4.3 también aparece indicada la magnitud δ=n· l , donde l es un vector que une los nodos de dos celdas contiguas. Por tanto, la magnitud escalar δfrepresenta la distancia que existe entre los centros de dos celdas colindantes en la dirección normal de la cara que las separa. Este término aparece en la formulación del FVM por lo que también es necesario calcularlo y almacenarlo. 1%% CALCULO DE DELTA 2 3delta = zeros(1,n_f); 4for iface = 1:n_f 5c1 = Lfc(iface,1); 6c2 = Lfc(iface,2); 7for ic = 1:n_c 8for k = 1:3 9if Lcf(ic,k) == iface 10 nx1 = nx(ic,k); 11 ny1 = ny(ic,k); 12 end 13 end 14 end 15 if c1 == c2 % Cara Frontera 16 dx1 = xf(iface) - xc(c1); 17 dy1 = yf(iface) - yc(c1); 18 delta(iface) = abs(dx1*nx1 + dy1*ny1); 19 else 20 dx1 = xc(c2) - xc(c1); 21 dy1 = yc(c2) - yc(c1); 22 delta(iface) = abs(dx1*nx1 + dy1*ny1); 23 end 24 end El cálculo del volumen de la celda es otro dato geométrico necesario para la resolución de las ecuaciones que gobiernan el problema. En el caso 2D, este valor coincide con el área de la celda. Para el mallado bidimensional triangular en cuestión, puede hacerse uso de la definición geométrica del área de un paralelogramo para realizar éste cálculo. Así, el área del triángulo será la mitad del área del paralelogramo formado por dos vectores directores cualesquiera. Es decir: Vc=1 2| t1× t2|=1 2(tx,1ty,2−ty,1tx,2)(4.7) Es importante puntualizar que en este caso los vectores t1 y t2 no son unitarios. El código que computa estos cálculos se expone a continuación: 1%% CALCULO DEL VOLUMEN 2 3vol(:) = 0; 4for ic = 1:n_c 5v1 = Lcv(ic,1); 6v2 = Lcv(ic,2); 7v3 = Lcv(ic,3); 8dx1 = xv(v2) - xv(v1); 9dy1 = yv(v2) - yv(v1); 10 dx2 = xv(v3) - xv(v1); 11 dy2 = yv(v3) - yv(v1); 12 vol(ic) = 0.5*(dx1*dy2 - dy1*dx2); % A = 1/2*|t1 x t2|
4.1 Problema difusivo, estacionario y unidimensional 73 13 end En la formulación del FVM no estructurado es frecuente la necesidad de calcular los valores de alguna propiedad φ en las caras y vértices a partir de sus valores en los nodos. Por ejemplo, véase la expresión (4.3), donde se requiere conocer Γf . La forma más usual de llevar a cabo este cálculo es mediante interpolaciones lineales. Este interpolante se obtiene normalmente como una función que pondera la distancia. Como la obtención de estas funciones interpoladoras es puramente geométrica, es recomendable realizar este cómputo como parte del pre-procesamiento y almacenar la información para cuando se necesite. Así, el valor de φ en las caras puede expresarse de la siguiente manera: φf=wfφ1+(1−wf)φ2donde wf=1/d1 1/d1+1/d2 (4.8) En la expresión de la función interpolante ( wf ), d1 y d2 son las distancias entre los nodos de las celdas 1 y 2 y el centro de la cara que las separa, respectivamente. Esto se ilustra en la Figura 4.4. De esta forma, cuando d1→0 ó d2→∞ , entonces φf→φ1 . Contrariamente, si d1→∞ ó d2→0 , entonces φf→φ2 . En esencia, la ecuación (4.8) pretende expresar cómo la influencia de un nodo contiguo a la cara decrece con la distancia a dicha cara. Esta interpolación ponderada inversamente proporcional a la distancia no es la única opción que puede tomarse para realizar este cálculo. Existen ciertos generadores de malla que utilizan funciones interpoladoras inversamente proporcionales al cuadrado de la distancia. Este procedimiento alternativo, derivado de ciertas leyes de la física cuyas magnitudes decaen con el cuadrado de la distancia, tiene el mismo comportamiento que el detallado en el párrafo anterior, con la salvedad de que ahora el decaimiento de la influencia con la distancia es mucho más acusado. No obstante, se usará la primera opción por ser la elección mayoritariamente utilizada por los softwares CFD actuales. Se añade subsiguientemente el código que computa los valores de Γ en las caras. 1%% FUNCIONES DE INTERPOLACION NODO-A-CARA 2 3wfun = zeros(1,n_f); 4gammaf = zeros(1,n_f); 5for iface = 1:n_f 6c1 = Lfc(iface,1); 7c2 = Lfc(iface,2); 8dx1 = xf(iface) - xc(c1); 9dy1 = yf(iface) - yc(c1); 10 d1=sqrt(dx1^2 + dy1^2); 11 dx2 = xf(iface) - xc(c2); 12 dy2 = yf(iface) - yc(c2); 13 d2=sqrt(dx2^2 + dy2^2); 14 wfun(iface) = d2 / (d1+d2); 15 gammaf(iface) = wfun(iface)*gamma(ic1) + (1-wfun(iface))*gamma(ic2); 16 end Por otro lado, para el cálculo de una propiedad φ en los vértices, éste puede computarse empleando funciones interpoladoras con un procedimiento similar. La disimilitud con el caso anterior es que ahora un vértice puede estar influenciado por N celdas, y no únicamente por dos como en el caso de las caras. Por tanto, los interpolantes tienen la siguiente forma: φv=φi/di N ∑ i=1 1/di = N ∑ i=1 wv,iφidonde wv,i=1/di N ∑ i=1 1/di (4.9)
74 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado Se expone a continuación el código MATLAB que calcula el valor de φen los vértices. 1%% FUNCIONES DE INTERPOLACION NODO-A-VERTICE 2 3cwfun = zeros(n_c,3); 4sumwt = zeros(1,n_v); 5for ic = 1:n_c 6for iv = 1:3 7v1 = Lcv(ic,iv); 8d1=sqrt((xv(v1) - xc(ic))^2 + (yv(v1) - yc(ic))^2); 9cwfun(ic,iv) = 1 / d1; 10 sumwt(v1) = sumwt(v1) + 1/d1; 11 end 12 end 13 14 for ic = 1:n_c 15 for iv = 1:3 16 v1 = Lcv(ic,iv); 17 cwfun(ic,iv) = cwfun(ic,iv) / sumwt(v1); 18 end 19 end 20 21 % Cálculo de \phi en Vértices 22 23 phiv = zeros(1,n_v); 24 for ic = 1:n_c 25 for k = 1:3 26 iv = Lcv(ic,k); 27 if bvert(iv) == 0 28 phiv(iv) = phiv(iv) + phi_ini(ic)*cwfun(ic,k); 29 end 30 end 31 end 11 2 1 2 f v 2 3 4 Figura 4.4 Esquema paradigmático para la obtención de las funciones interpoladoras wfywv,i. Por último y para terminar el procesamiento de la información relativa a la geometría se detallará el tratamiento de las condiciones de contorno, un concepto inherentemente relacionado con ésta. A la hora de almacenar la información en una formulación no estructurada las condiciones de contorno son tratadas de forma separada e independiente al resto, y también lo es su almacenamiento. Existen esencialmente tres razones de peso para hacer esto. En primer lugar, las condiciones de contorno solo se aplican en las caras fronterizas, y por tanto
4.1 Problema difusivo, estacionario y unidimensional 75 solo deben ser almacenadas en dichas caras. Debido a que el número total de celdas normalmente excede significativamente el número de celdas frontera, una formulación en la que se tengan en cuenta todas las celdas implicaría un derroche de recursos computacionales considerable. En segundo lugar, aislar la información relativa a las condiciones de contorno permite que el almacenamiento sea más compacto y rápido de acceder para la máquina, pues a la hora de barrer la memoria siempre es más eficiente especificar un punto de partida y la longitud de la información que se requiere tratar. Así se evita que el puntero pierda tiempo barriendo celdas interiores. Esta ventaja está directamente relacionada con la arquitectura interna de la máquina. En tercer lugar, se evita la necesidad de comprobar constantemente si la cara en cuestión es frontera o no. La ventaja de esto radica en que en la generación de códigos a gran escala el uso de comprobaciones IF dentro de bucles es prohibitivamente caro. Para resolver toda esta problemática, es necesario almacenar los siguientes arrays de información: •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) Además, para identificar si una cara está en la frontera o no será necesario un vector de valores binarios que al recorrer las caras devuelva 1 si la cara está en la frontera y 0 en caso contrario. Se expone seguidamente el código que lleva a cabo la numeración de las caras limítrofes y su almacenamiento. Esto será clave para la implementación de las condiciones de contorno del problema, en las que se entrará en detalle más adelante. 1%% INDEXACION DE CARAS FRONTERA 2 3n_f_b = 2*(M+N); % Total de Caras en Contorno 4bface = zeros(1,n_f); % 1 SI Cara en Contorno, 0 SI Cara Interior 5f_to_bf = zeros(1,n_f); 6bf_to_f = zeros(1,n_f_b); 7iface_b = 0; 8 9for i = 1:M 10 ic = i; 11 iface_b = iface_b + 1; 12 bf_to_f(iface_b) = Lcf(ic,1); 13 bface(Lcf(ic,1)) = 1; % Cara en Contorno 14 f_to_bf(Lcf(ic,1)) = iface_b; 15 end 16 for j = 1:N 17 ic = M + 2*(j-1)*M; 18 iface_b = iface_b + 1; 19 bf_to_f(iface_b) = Lcf(ic,2); 20 bface(Lcf(ic,2)) = 1; % Cara en Contorno 21 f_to_bf(Lcf(ic,2)) = iface_b; 22 end 23 for i = 1:M 24 ic = n_c - (i-1); 25 iface_b = iface_b + 1; 26 bf_to_f(iface_b) = Lcf(ic,2); 27 bface(Lcf(ic,2)) = 1; % Cara en Contorno 28 f_to_bf(Lcf(ic,2)) = iface_b;
76 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado 29 end 30 for j = 1:N 31 ic = n_c - (M-1) - 2*(j-1)*M; 32 iface_b = iface_b + 1; 33 bf_to_f(iface_b) = Lcf(ic,3); 34 bface(Lcf(ic,3)) = 1; % Cara en Contorno 35 f_to_bf(Lcf(ic,3)) = iface_b; 36 end Llegados a este punto, el procesamiento y almacenamiento de las características geométricas está completo. En la próxima sección se retomará la ecuación principal del FVM en mallas no estructuradas presentada en (4.2) y se buscará la manera de expresar todas las magnitudes en las caras en términos de los valores en los nodos con el fin de completar la formulación. 4.1.2 Discretización de las ecuaciones. Flujos normales y tangenciales Hasta ahora se han centrado los esfuerzos en procesar y almacenar toda la información inherente al mallado del dominio computacional, de forma que tenga que computarse una sola vez y sea accesible para su fin último: aproximar con la mayor rigurosidad y eficiencia posible la ecuación que gobierna el problema objeto de estudio, representada en (4.2). De todos los términos que componen la ecuación, ya se han desarrollado expresiones robustas para la obtención de Γf , Af y ∆V . Puesto que el término relativo a la fuente depende del problema (y en el caso particular del que se está estudiando su valor es nulo), restaría obtener una expresión para (∇φ)f·nf en función de los valores en los nodos. El desarrollo que se lleva a cabo a continuación se hará en pos de dicha finalidad. El primer paso a tomar es tener en cuenta que el gradiente de un campo escalar φ puede reescribirse en términos de cualquier sistema de coordenadas ortogonal sin pérdida de generalidad. Así, la siguiente igualdad se cumple para el caso 2D: ∇φ=∂φ ∂x i+∂φ ∂y j=∂φ ∂t t+∂φ ∂nn(4.10) Donde t y n son los vectores tangentes unitarios a la cara en cuestión, según el criterio especificado en la sección anterior ( n es la normal exterior a la celda y t es el tangente por el cual se verifica que n× t= k). Además, nótese que ∂φ/∂n= (∇φ)·ny∂φ/∂t= (∇φ)· t, luego: (∇φ)f=(∇φ)f·nfnf+(∇φ)f· tf tf(4.11) Si se multiplica la expresión (4.11) escalarmente por l (donde l es el vector representado en la Figura 4.3), se llega a: (∇φ)f· l=(∇φ)f·nfδ+(∇φ)f· tf tf· lf(4.12) Para expresar (∇φ)f· lf en términos de los valores en los nodos es necesario realizar una expansión en serie de Taylor de dos variables (debido a que la malla no es ortogonal). Esto es:
4.1 Problema difusivo, estacionario y unidimensional 77 φ2=φf+∂φ ∂xf (x2−xf)+ ∂φ ∂yf (y2−yf)+ 1 2 ∂2φ ∂x2f (x2−xf)2+1 2 ∂2φ ∂y2f (y2−yf)2+∂2φ ∂x∂yf (x2−xf)(y2−yf)+ ... φ1=φf+∂φ ∂xf (x1−xf)+ ∂φ ∂yf (y1−yf)+ 1 2 ∂2φ ∂x2f (x1−xf)2+1 2 ∂2φ ∂y2f (y1−yf)2+∂2φ ∂x∂yf (x1−xf)(y1−yf)+ ... (4.13) Si se sustraen las dos expresiones de (4.32) se llega a lo siguiente: φ2−φ1=∂φ ∂xf (x2−x1)+ ∂φ ∂yf (y2−y1)+ 1 2 ∂2φ ∂x2f(x2−xf)2−(x1−xf)2+1 2 ∂2φ ∂y2f(y2−yf)2−(y1−yf)2+ ∂2φ ∂x∂yf(x2−xf)(y2−yf)−(x1−xf)(y1−yf)+... (4.14) Si ahora se tiene en cuenta que lf= (x2−x1) i+(y2−y1) j, entonces puede decirse que: φ2−φ1≈(∇φ)f· lf(4.15) En dicha aproximación se comete un error truncado de valor ε que es difícil de estimar, a diferencia de lo que sucedía en mallas cartesianas donde ε∼O(∆x2) ó ε∼O(∆y2) . Sustituyendo (4.15) en (4.12) se obtiene finalmente: (∇φ)f·nf=φ2−φ1 δf −(∇φ)f· tf tf· lf δf =φ2−φ1 δf −Jt,f δf (4.16) La expresión (4.16) representa en esencia una descomposición del flujo normal a la cara f en dos componentes: una en la dirección del vector l y otra en dirección del vector t . En el caso particular en el que l y n sean colineales, entonces t· l=0 y la formulación pasa a ser formalmente idéntica al caso de mallas estructuradas cartesianas ya estudiado. Es importante destacar que una malla cartesiana es condición suficiente pero no necesaria para cancelar el término tangencial; pueden existir mallas no estructuradas donde dicho término se anule (por ejemplo, una malla formada en su totalidad por triángulos equiláteros), la única condición es que los vectores l y n estén alineados. Esto interesa desde el punto de vista computacional, pues puede demostrarse que la existencia de dicho término en las ecuaciones dificulta la convergencia de la solución. Esto ejemplifica el enorme interés que se tiene en una correcta generación de malla y da sentido a la existencia de campos completos de las matemáticas dedicados a la generación óptima de mallas (por ejemplo, las triangulaciones de Delaunay). Para completar el desarrollo de las ecuaciones sin pérdida de generalidad, va a suponerse el escenario en que el término asociado al flujo tangencial es no nulo. En tal caso, su cálculo requiere de un procedimiento completamente análogo al efectuado en (4.32). Así, se llega a: (∇φ)f· tf≈φa−φb | tf|(4.17)
78 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado Donde φa y φb son los valores de la propiedad φ en los vértices de la cara f (véase la Figura 4.3) y | tf|=Af es la longitud de la cara (por ser el problema 2D). De esta forma, y disponiendo de la información geométrica necesaria para computar el producto escalar t· l , el término relativo al flujo tangencial tendrá la siguiente expresión: Jt,f=(∇φ)f· t t· l=φa−φb Af (xa−xb)(x2−x1)+(ya−yb)(y2−y1) Af (4.18) Sustituyendo la ecuación (4.18) en (4.16), y ésta a su vez en (4.3), se llega a la expresión general que gobierna el problema de difusión pura con un mallado no estructurado: Nf,P ∑ f=1 Γfφnb(f)−φP δf −φa(f)−φb(f) δfAf tf· lfAf+¯ Sφ∆V=0(4.19) Donde, para el caso en cuestión, Nf,P=3 en todas las celdas. Además, el subíndice nb hace referencia a los nodos vecinos (del inglés neighbour), especificando nb(f) para hacer referencia al nodo adyacente a la cara fque se está considerando en el sumatorio. Para completar la discretización del problema es necesario expresar de algún modo las variables en los vértices, pues son a priori desconocidas. Existen dos formas de acometer esta tarea. La primera de ellas consiste en expresar los valores en los vértices en términos de los valores en los nodos haciendo uso de las funciones interpoladoras wv,i , descritas en (4.9). Así, las únicas incógnitas que presentarían las ecuaciones discretas serían los valores en los nodos, que es la situación que se andaba buscando. Esto se conoce como procedimiento implícito. La segunda forma de abordar el tratamiento de los valores en los vértices es computarlos según los valores de la iteración anterior, de forma que pueden pasarse al segundo miembro de la ecuación como términos conocidos. Este es el llamado tratamiento explícito. Son varias las ventajas y desventajas de cada uno de estos métodos. El implícito tiene la ventaja de presentar una convergencia más rápida a la solución, pero la matriz de coeficiente perderá parcialmente su cualidad como matriz dispersa (cuanto más dispersa es la matriz, más eficientes son los algoritmos que se emplean en este tipo de matrices y menos memoria se requiere) al aparecer términos relativos a vecinos más lejanos en cada ecuación discreta. Además, dificulta la implementación de las condiciones de contorno. Por otro lado, el método explícito mantiene la matriz mayoritariamente dispersa al incluir únicamente las celdas inmediatamente vecinas, pero tiene una convergencia más lenta. Por las razones expuestas, en la industria del CFD el tratamiento explícito es de los dos el que tiene un uso más extendido, y será por tanto el que se utilice en adelante. Resulta lógico pensar que la convergencia del método explícito estará condicionada por la importancia relativa del término tangencial frente al normal. Esto es así porque la diferencia entre dos iteraciones sucesivas será más notable si el valor de Jt,f juega un papel determinante. De este sencillo análisis puede inferirse que cuanto menos alineados estén los vectores n y l , mayor será Jt,f y más pobre será la convergencia. En la literatura se emplea el término skewness para hacer referencia a la desalineación de estos vectores. Como ya se ha comentado en otras ocasiones, es de vital importancia que los generadores de malla tengan en cuenta este concepto durante el pre-procesamiento para que la computación posterior tenga mejores resultados. Así, el quid de la cuestión consistirá en evitar las celdas con ángulos extremadamente agudos u obtusos y tratar de que los 60◦en mallas triangulares y los 90◦en mallas cuadrangulares prevalezcan. Para terminar con el planteamiento de las ecuaciones, es necesario puntualizar que la expresión general (4.19) se verá modificada por las condiciones de contorno en las celdas limítrofes. Para
4.1 Problema difusivo, estacionario y unidimensional 79 ilustrar este tema, se ejemplificarán los tipos de condiciones de contorno analizando directamente el problema de la Figura 3.2. Cara norte. Flujo de calor impuesto La cara norte del dominio computacional está sujeta a la condición de flujo de calor impuesto de valor conocido (saliente). Matemáticamente, esto se expresa como (∇φ)N·nN=JN . Así, de acuerdo con la Figura 4.2, la ecuación en las celdas 17 y18 queda como sigue: Nf,P ∑ f=1 f=N Γf"φnb(f)−φP δf − φ0 a(f)−φ0 b(f) δfAf! tf· lf#Af+ΓNJNAN+¯ Sφ∆V=0(4.20) Cara oeste. Aislada térmicamente La cara oeste tiene la condición de aislamiento térmico (flujo difusivo nulo a través de esta cara). Este caso es similar al anterior, salvo que ahora (∇φ)W·nW=0 . Esto afecta a las celdas 4 y 10 según la Figura 4.2, y la expresión final a la que se llega es: Nf,P ∑ f=1 f=W Γf"φnb(f)−φP δf − φ0 a(f)−φ0 b(f) δfAf! tf· lf#Af+¯ Sφ∆V=0(4.21) Cara este. Flujo de calor impuesto Caso formalmente idéntico al que se tiene en la cara norte, aunque aquí el flujo es entrante. La condición se expresa como (∇φ)E·nE=JE y la ecuación matemática que se verifica en las celdas 9y15 es: Nf,P ∑ f=1 f=E Γf"φnb(f)−φP δf − φ0 a(f)−φ0 b(f) δfAf! tf· lf#Af+ΓEJEAE+¯ Sφ∆V=0(4.22) Cara sur. Temperatura impuesta En la cara sur se cumple la condición φ=φS . Este caso difiere ligeramente de los presentados anteriormente. De nuevo, se divide el sumatorio entre las caras interiores y la cara frontera, pero esta vez el valor (∇φ)S·nS no es conocido, sino que tiene que ser reescrito en términos del flujo normal y tangencial, tal y como se detalló en el desarrollo de la expresión general. Finalmente, se llega a la siguiente expresión, válida para las celdas 1y2: Nf,P ∑ f=1 f=S Γf"φnb(f)−φP δf − φ0 a(f)−φ0 b(f) δfAf! tf· lf#Af+ ΓS"φS−φP δS − φ0 a(S)−φ0 b(S) δSAS! tS· lS#AS+¯ Sφ∆V=0 (4.23) Es importante resaltar que en esta expresión el vector lS va desde el nodo de la celda P hasta el centro de la cara frontera y no hasta el nodo de la celda vecina, pues no existe una celda colindante en esa dirección. Puesto que el segundo término del primer miembro de la ecuación es conocido, deberá ser tratado dentro del vector de términos independientes. Especial mención merece el caso en que la condición de contorno es constante a lo largo de toda la frontera. En dicha situación,
86 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado La representación de las curvas de nivel, que simbolizan las líneas de temperatura constante, es muy buen indicativo de la validez de la solución. En efecto, se observa cómo la temperatura en la pared inferior anda rondando los 150 ◦C , tal y como debería ser por la condición de contorno impuesta en dicha cara. Por otro lado, la condición de aislamiento térmico (que matemáticamente se expresa como ∇φ·n=0 ) se advierte fácilmente en la cara oeste, pues las isotermas entran en ella de forma perpendicular. Estos comentarios pueden confirmarse con la Figura 4.7. 160 160 160 160 170 170 170 170 180 180 180 190 190 190 200 200 210 210 220 220 230 230 240 240 250 260 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Figura 4.7 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=50 yN=50. Para aumentar la aplicabilidad del código desarrollado y validar su robustez, se pasa a continuación a analizar dos problemas similares al recién estudiado pero modificando las condiciones de contorno. Resulta interesante destacar que en el capítulo de mallas estructuradas las ecuaciones impuestas por las condiciones de contorno debían cambiarse manualmente en cada celda. Al contrario de esto, la gestión y almacenamiento de la información elaborada en el código de mallas no estructuradas permite estudiar problemas fundamentalmente diferentes alterando únicamente unas pocas líneas de código, tal y como se verá a continuación. Por ejemplo, considérese el siguiente problema similar al descrito en la Sección 3.1.4. Se tiene una placa plana de k=1000W/(mK) y dimensiones L×H , con L=1m y H=1m y con las condiciones de contorno siguientes: TW=TN=TE=0◦C y TS=100 ◦C . Este escenario pretende imitar el problema descrito por [ 6 ]y[ 17 ]-[ 20 ]. En este caso, la solución que se obtiene es la que se observa en las Figuras 4.8 y 4.9. En efecto, los resultados son los que cabía esperar según las condiciones de contorno impuestas, y están en perfecta consonancia con los que proporciona la literatura.
4.1 Problema difusivo, estacionario y unidimensional 87 Figura 4.8 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=30 yN=30 en el problema analizado en [6] y [17]-[20]. Obsérvese a continuación el único cambio a realizar dentro del código para obtener los resultados expuestos. El resto de líneas del programa permanecen inalteradas. Notar que no se ha asignado un valor idénticamente igual a cero (sino un valor ínfimo) en las condiciones de contorno de φ=0 para ser consistentes con la arquitectura del código.
88 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado 1%% CAMBIOS PARA MODIFICACION DE CONDICIONES DE CONTORNO 2 3% Condiciones de Contorno (phi impuesto en todas las caras) 4 5phi_south = 100; 6phi_east = 0.001; 7phi_west = 0.001; 8phi_north = 0.001; 9phib = zeros(1,n_f_b); 10 for iface_b = 1:n_f_b 11 iface = bf_to_f(iface_b); 12 if xf(iface) < 0.01*L 13 phib(iface_b) = phi_west; 14 end 15 if xf(iface) > 0.99*L 16 phib(iface_b) = phi_east; 17 end 18 if yf(iface) < 0.01*H 19 phib(iface_b) = phi_south; 20 end 21 if yf(iface) > 0.99*H 22 phib(iface_b) = phi_north; 23 end 24 end 25 26 % Cálculo de phi en Vértices 27 28 if flag == 1 29 phiv = zeros(1,n_v); 30 for ic = 1:n_c 31 for k = 1:3 32 iv = Lcv(ic,k); 33 if bvert(iv) == 0 34 phiv(iv) = phiv(iv) + phi_ini(ic)*cwfun(ic,k); 35 else 36 if xv(iv) < 0.01*L 37 phiv(iv) = phi_west; 38 elseif xv(iv) > 0.99*L 39 phiv(iv) = phi_east; 40 elseif yv(iv) < 0.01*H 41 phiv(iv) = phi_south; 42 elseif yv(iv) > 0.99*H 43 phiv(iv) = phi_north; 44 else 45 phiv(iv) = phiv(iv) + phi_ini(ic)*cwfun(ic,k); 46 end 47 end 48 end 49 end 50 end
4.1 Problema difusivo, estacionario y unidimensional 89 5 5 5 5 10 10 10 10 15 15 15 15 20 20 20 25 25 25 30 30 30 40 40 40 50 50 50 60 60 60 70 70 70 80 80 90 90 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Figura 4.9 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=30 yN=30 en el problema analizado en [6] y [17]-[20]. Por si esto no confiriese suficiente validez al algoritmo desarrollado, se estudiará subsiguientemente uno de los problemas propuestos en [ 4 ] para la exposición del FVM en geometrías bidimensionales. Dicho problema es nuevamente similar al descrito en la Sección 3.1.4. Consiste en una placa bidimensional con k=1000W/(mK) y de dimensiones L×H , con L=0.3m y H=0.4m . Las condiciones de contorno se especifican en la siguiente figura: 12 3 456 789 10 11 12 Figura 4.10 Representación del problema descrito en [4].
90 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado Nuevamente, modificando únicamente las condiciones de contorno, se obtienen lo siguientes resultados: Figura 4.11 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=30 yN=40 en el problema analizado en [4].
4.2 Problema convectivo-difusivo, estacionario y unidimensional 91 110 110 125 125 140 140 155 155 155 170 170 170 185 185 185 200 200 200 215 215 230 230 245 260 0.05 0.1 0.15 0.2 0.25 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Figura 4.12 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=30 yN=40 en el problema analizado en [4]. Los resultados son congruentes con la física del problema (véase la perpendicularidad de las isotermas en las caras aisladas y el cumplimiento de la condición de temperatura impuesta) y coinciden exactamente con los que reporta la literatura en cuestión. Con todo esto, se ha evidenciado que la formulación del FVM con mallas no estructuradas para problemas difusivos es totalmente correcta. Se dará un paso más en la formulación añadiendo el término advectivo en la próxima sección. 4.2 Problema convectivo-difusivo, estacionario y unidimensional En busca de la completitud de la caracterización del FVM en mallas no estructuradas, se dedica esta sección a la generalización de la formulación para problemas convectivo-difusivos. Como ya se ha visto en capítulos anteriores, se retoma a continuación la ecuación general que gobierna un problema convectivo-difusivo no estacionario en forma diferencial, que ya fue expuesta en (1.5): ∂(ρφ) ∂t+∇·(ρvφ) = ∇·(Γ∇φ)+Sφ(4.28) Como ya se sabe, el FVM consiste en integrar la ecuación (4.28) sobre un volumen de control genérico, llegando a la ecuación (1.6). Puesto que ya se han analizado con minuciosidad cada uno de los términos de esta expresión, se tratará a continuación de obtener únicamente una expresión útil para el flujo advectivo, evitando así ser redundantes. El término transitorio se supondrá nulo en este capítulo. No obstante, su caracterización fue ya estudiada en profundidad en la Sección 2.3. A modo de sucinto resumen, puede decirse que su aproximación se sustenta en tomar el valor medio y adopta la siguiente expresión: ZVC ∂ ∂t(ρφ)dV =∂(ρφ) ∂tP ∆V=(ρφ)P−(ρφ)0 P ∂t∆V(4.29)
92 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado Integrando con respecto al tiempo y aplicando alguno de los esquemas conocidos (implícito, explícito o Crank-Nicolson) se tendría completamente caracterizado este término. Volviendo al término convectivo, su integración lleva a la siguiente expresión: ZVC ∇·(ρvφ) = ZS (ρvφ)·ndA = Nf,P ∑ f=1 (ρv·nA)fφf= Nf,P ∑ f=1 Gfφf(4.30) Donde se ha aplicado el Teorema de la Divergencia de Gauss y se ha tenido en cuenta que la malla considerada está constituida por celdas delimitadas por segmentos rectos. Además, se ha definido el gasto másico que atraviesa cada cara como Gf= (ρv·nA)f , donde n es la normal externa. Notar también que nuevamente en el caso de malla triangular se tiene que Nf,P=3 para cualquier celda. El siguiente paso a seguir consiste en computar el valor del flujo Gf . Suponiendo que el campo de velocidades es un dato del problema (si no lo fuera, habría que resolver primero las ecuaciones del flujo), los términos que influyen en el cálculo del gasto son magnitudes conocidas. Por tanto, pueden obtenerse como parte del pre-procesamiento. El código es considerablemente sencillo y se adjunta a continuación. 1%% GASTO MASICO 2 3rho = 1000; 4cp = 4200; 5um = 2.92*10^-7; 6u=zeros(1,n_f); 7G=zeros(n_c,3); 8for ic = 1:n_c 9for k = 1:3 10 iface = Lcf(ic,k); 11 u(iface) = 6*um*((yf(iface)/H)-(yf(iface)/H)^2); 12 G(ic,k) = rho*cp*u(iface)*area(iface)*nx(ic,k); 13 end 14 end Téngase en cuenta que se están tomando los mismos datos de um , ρ y cp que los que se usaron en el problema descrito en la Sección 3.2, así como la expresión ( 3.25 ) para el cómputo de la velocidad. Una vez calculado el gasto másico, la siguiente y última tarea consiste en expresar el valor de φf en términos de los valores en los nodos. La Sección 2.2 estuvo dedicada específicamente a presentar y desarrollar algunos de los métodos que comúnmente se utilizan en la industria para llevar a cabo esta tarea, propia del término advectivo. A modo ilustrativo, en esta sección se discutirán únicamente el método de las diferencias centradas y el esquema upwind, dos de los más utilizados y a la vez más sencillos de implementar, sin descuidar la precisión de los resultados. En el método de las diferencias centradas, el valor de φf se obtiene ponderando los valores de φ en los nodos de las celdas adyacentes a la cara con la función de interpolación wf , cuyo cálculo ya se presentó en la Sección 4.1.1. Así, se verifica la siguiente igualdad: Nf,P ∑ f=1 Gfφf= Nf,P ∑ f=1 GfwfφP+(1−wf)φnb(f)(4.31) Las ventajas y desventas de este método ya han sido descritas en numerosas ocasiones. Por un lado, el error de truncamiento decae cuadráticamente con la distancia entre el nodo y la cara, lo cual es un aspecto positivo; por otro lado, no se tiene en cuenta que el proceso de convección transporta la información aguas abajo acentuando la influencia en uno de los nodos según la dirección del flujo,
4.2 Problema convectivo-difusivo, estacionario y unidimensional 93 lo cual puede llevar a resultados incorrectos u oscilaciones alrededor de la solución. No obstante, si la convección en el problema objeto de estudio no es lo suficientemente fuerte como para penalizar este método frente al upwind, ambos serán buenas opciones a elegir. En el esquema upwind, φf tomará el valor de una de las celdas colindantes según la dirección del flujo. Una forma expresar φfes como una función a trozos: φf= φPsi Gf>0 φnb si Gf<0(4.32) Sin embargo, esta opción obliga al código a realizar numerosas comprobaciones IF para discernir si se está en un caso o en otro, lo cual es, como ya se sabe, una operación computacionalmente costosa. Una alternativa es definir la función de forma que la lógica anteriormente descrita se cumpla automáticamente. Esto se consigue de la siguiente forma: Nf,P ∑ f=1 Gfφf= Nf,P ∑ f=1Gf+|Gf| 2φP+Gf−|Gf| 2φnb(f)(4.33) Nuevamente, este método tiene sus pros y sus contras. Mientras que la física que subyace al método logra ser más realista y su implementación es tremendamente sencilla, la precisión con la que este esquema resuelve el problema es numéricamente peor. La razón de ello es que el error de truncamiento cae con la distancia entre el nodo y la cara, y no con su cuadrado como ocurría en el método de las diferencias centradas. Además, puede llevar a problemas de falsa difusión, como ya se ha comentado en anteriores ocasiones. No obstante, este método resulta perfectamente válido para el caso que nos ocupa, como se verá inmediatamente. Por último, si a la ecuación (4.19) se le añade el término convectivo, se llega a la siguiente ecuación general para problemas convectivo-difusivos estacionarios en mallas no estructuradas: Nf,P ∑ f=1Gf+|Gf| 2φP+Gf−|Gf| 2φnb(f)= Nf,P ∑ f=1 Γfφnb(f)−φP δf −φa(f)−φb(f) δfAf tf· lfAf+¯ Sφ∆V (4.34) El ensamblaje del sistema es muy similar al que se expuso en (4.24). No obstante, el valor de los coeficientes necesita ser modificado para incluir el término relativo a la convección. Así, la forma final de los coeficientes de la matriz será: Ap,P= Nf,P ∑ f=1ΓfAf δf +Gf+|Gf| 2(4.35) Ak,P=−Γk,PAk,P δk,P +Gk,P−|Gk,P| 2(4.36) Askew,P=− Nf,P ∑ f=1 Γf φ0 a(f)−φ0 b(f) δf! tf· lf(4.37) Es interesante notar como el término Askew,P permanece inalterado. Esto es porque el proceso de convección no lleva asociado un flujo tangencial; éste es único del proceso de difusión. El motivo
94 Capítulo 4. FVM para problemas convectivo-difusivos en 2D. Mallado no estructurado de esto radica en que el flujo difusivo depende del gradiente de φ , mientras que el flujo convectivo no lo hace. Llegados a este punto, se está en condiciones de resolver un problema de convección-difusión empleando esta formulación. Para validar los resultados, se resolverá el mismo problema presentado en la Sección 3.2 y se compararán las soluciones. Puede encontrarse una representación esquemática de este problema en la Figura 3.7. Figura 4.13 Valores del campo de temperaturas en los nodos para un mallado no estructurado M×N con M=50 yN=10 utilizando el esquema upwind.
4.2 Problema convectivo-difusivo, estacionario y unidimensional 95 De forma equivalente al caso puramente difusivo, se añade a continuación una representación de las curvas de nivel como herramienta para validar los resultados. Igual que se comentó en el caso de malla estructurada, esta ilustración muestra lo que era de esperar: las isotermas entran perpendicularmente en las caras donde el flujo difusivo es nulo, y presentan cierta curvatura cóncava con respecto al flujo incidente como manifestación del proceso de convección. Además, la temperatura en la pared oeste roza los 30 ◦C, cumpliendo con la condición de contorno. 30 30 31 31 32 32 33 33 34 34 35 35 36 36 37 37 37 37 38 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Figura 4.14 Curvas de nivel (isotermas) del campo de temperaturas para un mallado no estructurado M×Ncon M=50 yN=10. El código ejecutado para obtener estos resultados es prácticamente idéntico al caso anterior. Como el pre-procesamiento es el mismo (añadiendo el cálculo del gasto Gf , que ya ha sido detallado), se presentan a continuación únicamente las líneas que difieren con respecto a la situación puramente difusiva. Destacar que los cálculos se han hecho empleando el esquema upwind para aproximar el término convectivo. 1%% CAMBIOS PARA INTRODUCCION DE TERMINO CONVECTIVO 2 3% Coeficientes de la Matriz A 4 5for ic = 1:n_c 6for k = 1:3 7iface = Lcf(ic,k); 8if bface(iface) == 0 9A0(ic) = A0(ic) + gammaf(iface)*area(iface)/delta(iface) + ... 10 (G(ic,k) + abs(G(ic,k)))/2; 11 Anb(ic,k) = -gammaf(iface)*area(iface)/delta(iface) + ... 12 (G(ic,k) - abs(G(ic,k)))/2; 13 ic1 = Lfc(iface,1); 14 ic2 = Lfc(iface,2); 15 if ic1 == ic 16 A(ic,ic2) = Anb(ic,k); 17 else
102 Capítulo 5. Conclusión y trabajos futuros de cálculo numérico, como son la propagación de errores o el tiempo computacional. Adicionalmente, se listan a continuación una serie de potenciales trabajos futuros que podrían extender el alcance de este proyecto: • Extensión de los métodos descritos a problemas tridimensionales. Generación de mallas hexaédricas estructuradas y tetraédricas no estructuradas, formulación de las ecuaciones del problema, resolución y valoración de los resultados. • Análisis y valoración de problemas convectivo-difusivos no estacionarios de carácter general en geometrías unidimensionales y bidimensionales. • Análisis del FVM en geometrías axilsimétricas. Formulación del problema convectivodifusivo en coordenadas cilíndricas, discretización y resolución. • Estudio de problemas bidimensionales con mallas curvilíneas. Formulación de las ecuaciones en las variables ξyηy cálculo de la solución. • Cómputo del campo de velocidades ( v ) como paso previo a la resolución del problema. Puede calcularse de forma autocontenida a través de la ecuación de cantidad de movimiento con las presiones y velocidades (se conoce como método SIMPLE) o mediante la Teoría Potencial resolviendo ∇2ψ=0, donde u=∂ψ/∂xyv=∂ψ/∂y. • Comparación del FVM con otros métodos similares en relación a determinados problemas convectivo-difusivos. Especial interés cobraría la comparación con el CVFEM, variante del método de los elementos finitos basado en la naturaleza conservativa de los métodos de volúmenes finitos. • Proceso de optimización de mallas. Triangulación de Delaunay. Permite obtener mallas triangulares en las que se maximiza el número de triángulos equiláteros, evitando algunas propiedades no deseadas que presentan los triángulos excesivamente esbeltos. • Implementación de los problemas analizados en algún software CFD usado actualmente en la industria (Ansys Fluent, OpenFOAM, STAR-CCM+, etc.). Comparación de la solución obtenida con los resultados extraídos de MATLAB.
Apéndice A Códigos de MATLAB A lo largo del documento se han introducido ciertas líneas de código que pretendían servir como elemento de ayuda en la exposición del tema en cuestión. La intención es que dichos códigos tuvieran una estructura especialmente clara y compacta para mantener la continuidad de la presentación. Sin embargo, muchos de los resultados que se muestran no aparecen respaldados por ningún código. Por ello, en este apéndice se adjuntan otros códigos MATLAB desarrollados por el autor que también resuelven ciertos problemas térmicos presentados a lo largo de la memoria. En particular, aquí se incluyen aquellos menos representativos, bien porque resultan redundantes o bien porque la claridad de su exposición puede no resultar tan representativa. Así, el lector dispone de una herramienta para verificar que los resultados que en este proyecto se presentan son acordes a los obtenidos. Además, pueden ser un buen punto de partida para lectores interesados en profundizar más en el tema principal de este proyecto. A.1 Convección-difusión estacionaria 1D 1%% COMPARACION DE TODOS LOS METODOS 2 3clear all;close all;clc 4 5% MALLADO 6 7L = 1; % Longitud Barra (m) 8N = 20; % Número de Nodos 9 10 caras = linspace(0,L,N+1); % Nodos Equiespaciados 11 12 % caras_theta = linspace(0,pi,N+1); % Nodos Concentrados en Extremos 13 % caras = L/2*(1 - cos(caras_theta)); 14 15 % rng = [0,L*rand(1,N-1),L]; % Nodos Aleatorios 16 % caras = sort(rng); 17 18 for i = 2:length(caras) 19 nodos(i-1) = (caras(i)-caras(i-1))/2 + caras(i-1); 20 end 21 22 % DATOS 23 24 Gamma_w = 0.1*ones(N,1); % Conductividad Térmica (W/(m*K)) -> CTE 103
104 Capítulo A. Códigos de MATLAB 25 Gamma_e = 0.1*ones(N,1); % Conductividad Térmica (W/(m*K)) -> CTE 26 Aw = 1*ones(N,1); % Ãrea Transversal (m^2) -> CTE 27 Ae = 1*ones(N,1); % Ãrea Transversal (m^2) -> CTE 28 nw = -ones(N,1); % Vector Normal (nw=-1 cte para caras west) 29 ne = ones(N,1); % Vector Normal (ne=1 cte para caras east) 30 rho = 1*ones(N,1); % Densidad del Fluido (kg/m^3) 31 u = 3*ones(N,1); % Velocidad Unidireccional (m/s) 32 Fw = rho.*Aw.*nw.*u; % Flujo de Fluido West 33 Fe = rho.*Ae.*ne.*u; % Flujo de Fluido East 34 35 %%% DIFERENCIAS CENTRADAS %%% 36 37 % COEFICIENTES 38 39 for I = 2:N-1 40 gw(I) = (caras(I)-nodos(I-1))/(nodos(I)-nodos(I-1)); % Factor Interpolación 41 ge(I) = (caras(I+1)-nodos(I+1))/(nodos(I)-nodos(I+1)); % Factor Interpolación 42 Dw(I) = Gamma_w(I)*nw(I)*Aw(I)/(nodos(I)-nodos(I-1)); 43 De(I) = Gamma_e(I)*ne(I)*Ae(I)/(nodos(I+1)-nodos(I)); 44 aW(I) = Fw(I)*(1-gw(I))+Dw(I); % Coeficiente phi_W 45 aP(I) = Fw(I)*gw(I)+Fe(I)*ge(I)-Dw(I)+De(I); % Coeficiente phi_P 46 aE(I) = Fe(I)*(1-ge(I))-De(I); % Coeficiente phi_E 47 aSP(I) = 0; % Término Asociado a la Generación 48 A(I,I-1) = aW(I); 49 A(I,I) = aP(I); 50 A(I,I+1) = aE(I); 51 b(I) = aSP(I); 52 end 53 54 % CONDICIONES DE CONTORNO 55 56 % LEFT 57 58 phi_0 = 1; % CdC Izquierda: phi impuesto (K); 59 gw(1) = 0; 60 ge(1) = (caras(2)-nodos(2))/(nodos(1)-nodos(2)); 61 Dw(1) = Gamma_w(1)*nw(1)*Aw(1)/(nodos(1)-caras(1)); 62 De(1) = Gamma_e(1)*ne(1)*Ae(1)/(nodos(2)-nodos(1)); 63 aW(1) = Fw(1)*(1-gw(1))+Dw(1); 64 aP(1) = Fw(1)*gw(1)+Fe(1)*ge(1)-Dw(1)+De(1); 65 aE(1) = Fe(1)*(1-ge(1))-De(1); 66 aSP(1) = 0; 67 A(1,1) = aP(1); 68 A(1,2) = aE(1); 69 b(1) = aSP(1)-aW(1)*phi_0; 70 71 % RIGHT 72 73 phi_L = 0; % CdC Izquierda: phi impuesto (K); 74 gw(N) = (caras(N)-nodos(N-1))/(nodos(N)-nodos(N-1)); 75 ge(N) = 0; 76 Dw(N) = Gamma_w(N)*nw(N)*Aw(N)/(nodos(N)-nodos(N-1)); 77 De(N) = Gamma_e(N)*ne(N)*Ae(N)/(caras(end)-nodos(N)); 78 aW(N) = Fw(N)*(1-gw(N))+Dw(N); 79 aP(N) = Fw(N)*gw(N)+Fe(N)*ge(N)-Dw(N)+De(N); 80 aE(N) = Fe(N)*(1-ge(N))-De(N); 81 aSP(N) = 0; 82 A(N,N-1) = aW(N); 83 A(N,N) = aP(N); 84 b(N) = aSP(N)-aE(N)*phi_L; 85 86 % SOLUCION
A.1 Convección-difusión estacionaria 1D 105 87 88 phi = A\b’; 89 plot([0,nodos,L],[phi_0;phi;phi_L],’o-’,’MarkerIndices’,2:length(nodos)+1,’linewidth’ ,1.2) 90 hold on 91 92 clearvars -except L N caras nodos Gamma_w Gamma_e Aw Ae nw ne rho u Fw Fe 93 94 %%% FIRST ORDER UPWIND %%% 95 96 % COEFICIENTES 97 98 for I = 2:N-1 99 gw(I) = 1-heaviside(u(I)); % Método Upwind 100 ge(I) = heaviside(u(I)); % Método Upwind 101 Dw(I) = Gamma_w(I)*nw(I)*Aw(I)/(nodos(I)-nodos(I-1)); 102 De(I) = Gamma_e(I)*ne(I)*Ae(I)/(nodos(I+1)-nodos(I)); 103 aW(I) = Fw(I)*(1-gw(I))+Dw(I); % Coeficiente phi_W 104 aP(I) = Fw(I)*gw(I)+Fe(I)*ge(I)-Dw(I)+De(I); % Coeficiente phi_P 105 aE(I) = Fe(I)*(1-ge(I))-De(I); % Coeficiente phi_E 106 aSP(I) = 0; % Término Asociado a la Generación 107 A(I,I-1) = aW(I); 108 A(I,I) = aP(I); 109 A(I,I+1) = aE(I); 110 b(I) = aSP(I); 111 end 112 113 % CONDICIONES DE CONTORNO 114 115 % LEFT 116 117 phi_0 = 1; % CdC Izquierda: phi impuesto (K) 118 gw(1) = 1-heaviside(u(1)); 119 ge(1) = heaviside(u(1)); 120 Dw(1) = Gamma_w(1)*nw(1)*Aw(1)/(nodos(1)-caras(1)); 121 De(1) = Gamma_e(1)*ne(1)*Ae(1)/(nodos(2)-nodos(1)); 122 aW(1) = Fw(1)*(1-gw(1))+Dw(1); 123 aP(1) = Fw(1)*gw(1)+Fe(1)*ge(1)-Dw(1)+De(1); 124 aE(1) = Fe(1)*(1-ge(1))-De(1); 125 aSP(1) = 0; 126 A(1,1) = aP(1); 127 A(1,2) = aE(1); 128 b(1) = aSP(1)-aW(1)*phi_0; 129 130 % RIGHT 131 132 phi_L = 0; % CdC Derecha: phi impuesto (K) 133 gw(N) = 1-heaviside(u(N)); 134 ge(N) = heaviside(u(N)); 135 Dw(N) = Gamma_w(N)*nw(N)*Aw(N)/(nodos(N)-nodos(N-1)); 136 De(N) = Gamma_e(N)*ne(N)*Ae(N)/(caras(end)-nodos(N)); 137 aW(N) = Fw(N)*(1-gw(N))+Dw(N); 138 aP(N) = Fw(N)*gw(N)+Fe(N)*ge(N)-Dw(N)+De(N); 139 aE(N) = Fe(N)*(1-ge(N))-De(N); 140 aSP(N) = 0; 141 A(N,N-1) = aW(N); 142 A(N,N) = aP(N); 143 b(N) = aSP(N)-aE(N)*phi_L; 144 145 % SOLUCION 146 147 phi = A\b’;
106 Capítulo A. Códigos de MATLAB 148 plot([0,nodos,L],[phi_0;phi;phi_L],’o-’,’MarkerIndices’,2:length(nodos)+1,’linewidth’ ,1.2) 149 hold on 150 151 clearvars -except L N caras nodos Gamma_w Gamma_e Aw Ae nw ne rho u Fw Fe 152 153 %%% SECOND ORDER UPWIND %%% 154 155 % COEFICIENTES 156 157 for I = 3:N-2 158 gwww(I) = -1/2*heaviside(u(I)); 159 gww(I) = 3/2*heaviside(u(I)); 160 gew(I) = -1/2*heaviside(u(I)); 161 gep(I) = 3/2*heaviside(u(I)); 162 gwp(I) = 3/2*(1-heaviside(u(I))); 163 gwe(I) = -1/2*(1-heaviside(u(I))); 164 gee(I) = 3/2*(1-heaviside(u(I))); 165 geee(I) = -1/2*(1-heaviside(u(I))); 166 Dw(I) = Gamma_w(I)*nw(I)*Aw(I)/(nodos(I)-nodos(I-1)); 167 De(I) = Gamma_e(I)*ne(I)*Ae(I)/(nodos(I+1)-nodos(I)); 168 aWW(I) = Fw(I)*gwww(I); % Coeficiente phi_WW 169 aW(I) = Fw(I)*gww(I)+Fe(I)*gew(I)+Dw(I); % Coeficiente phi_W 170 aP(I) = Fw(I)*gwp(I)+Fe(I)*gep(I)-Dw(I)+De(I); % Coeficiente phi_P 171 aE(I) = Fe(I)*gee(I)+Fw(I)*gwe(I)-De(I); % Coeficiente phi_E 172 aEE(I) = Fe(I)*geee(I); % Coeficiente phi_EE 173 aSP(I) = 0; % Término Asociado a la Generación 174 A(I,I-2) = aWW(I); 175 A(I,I-1) = aW(I); 176 A(I,I) = aP(I); 177 A(I,I+1) = aE(I); 178 A(I,I+2) = aEE(I); 179 b(I) = aSP(I); 180 end 181 182 % CONDICIONES DE CONTORNO 183 184 % LEFT 185 186 phi_0 = 1; % CdC Izquierda: phi impuesto (K) 187 188 gwww(1) = 0; 189 gww(1) = 1; 190 gew(1) = -1; 191 gep(1) = 2; 192 gwp(1) = 0; 193 gwe(1) = 0; 194 gee(1) = 0; 195 geee(1) = 0; 196 Dw(1) = Gamma_w(1)*nw(1)*Aw(1)/(nodos(1)-caras(1)); 197 De(1) = Gamma_e(1)*ne(1)*Ae(1)/(nodos(2)-nodos(1)); 198 aWW(1) = Fw(1)*gwww(1); % Coeficiente phi_WW 199 aW(1) = Fw(1)*gww(1)+Fe(1)*gew(1)+Dw(1); % Coeficiente phi_W 200 aP(1) = Fw(1)*gwp(1)+Fe(1)*gep(1)-Dw(1)+De(1); % Coeficiente phi_P 201 aE(1) = Fe(1)*gee(1)+Fw(1)*gwe(1)-De(1); % Coeficiente phi_E 202 aEE(1) = Fe(1)*geee(1); % Coeficiente phi_EE 203 aSP(1) = 0; % Término Asociado a la Generación 204 A(1,1) = aP(1); 205 A(1,2) = aE(1); 206 A(1,3) = aEE(1); 207 b(1) = aSP(1)-aW(1)*phi_0; 208
A.1 Convección-difusión estacionaria 1D 107 209 gwww(2) = -1; 210 gww(2) = 2; 211 gew(2) = -1/2; 212 gep(2) = 3/2; 213 gwp(2) = 0; 214 gwe(2) = 0; 215 gee(2) = 0; 216 geee(2) = 0; 217 Dw(2) = Gamma_w(2)*nw(2)*Aw(2)/(nodos(2)-nodos(1)); 218 De(2) = Gamma_e(2)*ne(2)*Ae(2)/(nodos(3)-nodos(2)); 219 aWW(2) = Fw(2)*gwww(2); % Coeficiente phi_WW 220 aW(2) = Fw(2)*gww(2)+Fe(2)*gew(2)+Dw(2); % Coeficiente phi_W 221 aP(2) = Fw(2)*gwp(2)+Fe(2)*gep(2)-Dw(2)+De(2); % Coeficiente phi_P 222 aE(2) = Fe(2)*gee(2)+Fw(2)*gwe(2)-De(2); % Coeficiente phi_E 223 aEE(2) = Fe(2)*geee(2); % Coeficiente phi_EE 224 aSP(2) = 0; % Término Asociado a la Generación 225 A(2,1) = aW(2); 226 A(2,2) = aP(2); 227 A(2,3) = aE(2); 228 A(2,4) = aEE(2); 229 b(2) = aSP(2)-aWW(2)*phi_0; 230 231 % RIGHT 232 233 phi_L = 0; % CdC Derecha: phi impuesto (K) 234 235 gwww(N-1) = -1/2*heaviside(u(N-1)); 236 gww(N-1) = 3/2*heaviside(u(N-1)); 237 gew(N-1) = -1/2*heaviside(u(N-1)); 238 gep(N-1) = 3/2*heaviside(u(N-1)); 239 gwp(N-1) = 3/2*(1-heaviside(u(N-1))); 240 gwe(N-1) = -1/2*(1-heaviside(u(N-1))); 241 gee(N-1) = 3/2*(1-heaviside(u(N-1))); 242 geee(N-1) = -1/2*(1-heaviside(u(N-1))); 243 Dw(N-1) = Gamma_w(N-1)*nw(N-1)*Aw(N-1)/(nodos(N-1)-nodos(N-2)); 244 De(N-1) = Gamma_e(N-1)*ne(N-1)*Ae(N-1)/(nodos(N)-nodos(N-1)); 245 aWW(N-1) = Fw(N-1)*gwww(N-1); % Coeficiente phi_WW 246 aW(N-1) = Fw(N-1)*gww(N-1)+Fe(N-1)*gew(N-1)+Dw(N-1); % Coeficiente phi_W 247 aP(N-1) = Fw(N-1)*gwp(N-1)+Fe(N-1)*gep(N-1)-Dw(N-1)+De(N-1); % Coeficiente phi_P 248 aE(N-1) = Fe(N-1)*gee(N-1)+Fw(N-1)*gwe(N-1)-De(N-1); % Coeficiente phi_E 249 aEE(N-1) = Fe(N-1)*geee(N-1); % Coeficiente phi_EE 250 aSP(N-1) = 0; % Término Asociado a la Generación 251 A(N-1,N-3) = aWW(N-1); 252 A(N-1,N-2) = aW(N-1); 253 A(N-1,N-1) = aP(N-1); 254 A(N-1,N) = aE(N-1); 255 b(N-1) = aSP(N-1); 256 257 gwww(N) = -1/2*heaviside(u(N)); 258 gww(N) = 3/2*heaviside(u(N)); 259 gew(N) = -1/2*heaviside(u(N)); 260 gep(N) = 3/2*heaviside(u(N)); 261 gwp(N) = 0; 262 gwe(N) = 0; 263 gee(N) = 0; 264 geee(N) = 0; 265 Dw(N) = Gamma_w(N)*nw(N)*Aw(N)/(nodos(N)-nodos(N-1)); 266 De(N) = Gamma_e(N)*ne(N)*Ae(N)/(caras(end)-nodos(N)); 267 aWW(N) = Fw(N)*gwww(N); % Coeficiente phi_WW 268 aW(N) = Fw(N)*gww(N)+Fe(N)*gew(N)+Dw(N); % Coeficiente phi_W 269 aP(N) = Fw(N)*gwp(N)+Fe(N)*gep(N)-Dw(N)+De(N); % Coeficiente phi_P 270 aE(N) = Fe(N)*gee(N)+Fw(N)*gwe(N)-De(N); % Coeficiente phi_E
108 Capítulo A. Códigos de MATLAB 271 aEE(N) = Fe(N)*geee(N); % Coeficiente phi_EE 272 aSP(N) = 0; % Término Asociado a la Generación 273 A(N,N-2) = aWW(N); 274 A(N,N-1) = aW(N); 275 A(N,N) = aP(N); 276 b(N) = aSP(N)-aE(N)*phi_L; 277 278 % SOLUCION 279 280 phi = A\b’; 281 plot([0,nodos,L],[phi_0;phi;phi_L],’o-’,’MarkerIndices’,2:length(nodos)+1,’linewidth’ ,1.2) 282 hold on 283 284 clearvars -except L N caras nodos Gamma_w Gamma_e Aw Ae nw ne rho u Fw Fe 285 286 %%% POWER-LAW SCHEME %%% 287 288 % COEFICIENTES 289 290 for I = 2:N-1 291 Dw(I) = Gamma_w(I)*nw(I)*Aw(I)/(nodos(I)-nodos(I-1)); 292 De(I) = Gamma_e(I)*ne(I)*Ae(I)/(nodos(I+1)-nodos(I)); 293 Pe_w(I) = Fw(I)/Dw(I); 294 Pe_e(I) = Fe(I)/De(I); 295 beta_w(I) = (1-0.1*Pe_w(I))^5/Pe_w(I); 296 beta_e(I) = (1-0.1*Pe_e(I))^5/Pe_e(I); 297 298 if abs(Pe_w(I)) < 10 299 gwp(I) = -beta_w(I); 300 gww(I) = 1+beta_w(I); 301 elseif Pe_w(I) >= 10 302 gwp(I) = 0-1/Pe_w(I); 303 gww(I) = 1+1/Pe_w(I); 304 else 305 gwp(I) = 1-1/Pe_w(I); 306 gww(I) = 0+1/Pe_w(I); 307 end 308 309 if abs(Pe_e(I)) < 10 310 gep(I) = 1+beta_e(I); 311 gee(I) = -beta_e(I); 312 elseif Pe_e(I) >= 10 313 gep(I) = 1+1/Pe_e(I); 314 gee(I) = 0-1/Pe_e(I); 315 else 316 gep(I) = 0+1/Pe_e(I); 317 gee(I) = 1-1/Pe_e(I); 318 end 319 320 aW(I) = Fw(I)*gww(I); % Coeficiente phi_W 321 aP(I) = Fw(I)*gwp(I)+Fe(I)*gep(I); % Coeficiente phi_P 322 aE(I) = Fe(I)*gee(I); % Coeficiente phi_E 323 aSP(I) = 0; % Término Asociado a la Generación 324 A(I,I-1) = aW(I); 325 A(I,I) = aP(I); 326 A(I,I+1) = aE(I); 327 b(I) = aSP(I); 328 end 329 330 % CONDICIONES DE CONTORNO 331
A.1 Convección-difusión estacionaria 1D 109 332 % LEFT 333 334 phi_0 = 1; % CdC Izquierda: phi impuesto (K) 335 Dw(1) = Gamma_w(1)*nw(1)*Aw(1)/(nodos(1)-caras(1)); 336 De(1) = Gamma_e(1)*ne(1)*Ae(1)/(nodos(2)-nodos(1)); 337 Pe_w(1) = Fw(1)/Dw(1); 338 Pe_e(1) = Fe(1)/De(1); 339 beta_w(1) = (1-0.1*Pe_w(1))^5/Pe_w(1); 340 beta_e(1) = (1-0.1*Pe_e(1))^5/Pe_e(1); 341 342 if abs(Pe_w(1)) < 10 343 gwp(1) = -beta_w(1); 344 gww(1) = 1+beta_w(1); 345 elseif Pe_w(1) >= 10 346 gwp(1) = 0-1/Pe_w(1); 347 gww(1) = 1+1/Pe_w(1); 348 else 349 gwp(1) = 1-1/Pe_w(1); 350 gww(1) = 0+1/Pe_w(1); 351 end 352 353 if abs(Pe_e(1)) < 10 354 gep(1) = 1+beta_e(1); 355 gee(1) = -beta_e(1); 356 elseif Pe_e(1) >= 10 357 gep(1) = 1+1/Pe_e(1); 358 gee(1) = 0-1/Pe_e(1); 359 else 360 gep(1) = 0+1/Pe_e(1); 361 gee(1) = 1-1/Pe_e(1); 362 end 363 364 aW(1) = Fw(1)*gww(1); 365 aP(1) = Fw(1)*gwp(1)+Fe(1)*gep(1); 366 aE(1) = Fe(1)*gee(1); 367 aSP(1) = 0; 368 A(1,1) = aP(1); 369 A(1,2) = aE(1); 370 b(1) = aSP(1)-aW(1)*phi_0; 371 372 % RIGHT 373 374 phi_L = 0; % CdC Derecha: phi impuesto (K) 375 Dw(N) = Gamma_w(N)*nw(N)*Aw(N)/(nodos(N)-nodos(N-1)); 376 De(N) = Gamma_e(N)*ne(N)*Ae(N)/(caras(end)-nodos(N)); 377 Pe_w(N) = Fw(N)/Dw(N); 378 Pe_e(N) = Fe(N)/De(N); 379 beta_w(N) = (1-0.1*Pe_w(N))^5/Pe_w(N); 380 beta_e(N) = (1-0.1*Pe_e(N))^5/Pe_e(N); 381 382 if abs(Pe_w(N)) < 10 383 gwp(N) = -beta_w(N); 384 gww(N) = 1+beta_w(N); 385 elseif Pe_w(N) >= 10 386 gwp(N) = 0-1/Pe_w(N); 387 gww(N) = 1+1/Pe_w(N); 388 else 389 gwp(N) = 1-1/Pe_w(N); 390 gww(N) = 0+1/Pe_w(N); 391 end 392 393 if abs(Pe_e(N)) < 10
110 Capítulo A. Códigos de MATLAB 394 gep(N) = 1+beta_e(N); 395 gee(N) = -beta_e(N); 396 elseif Pe_e(N) >= 10 397 gep(N) = 1+1/Pe_e(N); 398 gee(N) = 0-1/Pe_e(N); 399 else 400 gep(N) = 0+1/Pe_e(N); 401 gee(N) = 1-1/Pe_e(N); 402 end 403 404 aW(N) = Fw(N)*gww(N); 405 aP(N) = Fw(N)*gwp(N)+Fe(N)*gep(N); 406 aE(N) = Fe(N)*gee(N); 407 aSP(N) = 0; 408 A(N,N-1) = aW(N); 409 A(N,N) = aP(N); 410 b(N) = aSP(N)-aE(N)*phi_L; 411 412 % SOLUCION 413 414 phi = A\b’; 415 plot([0,nodos,L],[phi_0;phi;phi_L],’o-’,’MarkerIndices’,2:length(nodos)+1,’linewidth’ ,1.2) 416 hold on 417 418 phiexact = @(x)(phi_0+(phi_L-phi_0)*(exp(rho(1)*u(1)*x/Gamma_e(1))-1)/... 419 ((exp(rho(1)*u(1)*L/Gamma_e(1))-1))); 420 x = [0:0.002:L]; 421 plot(x,phiexact(x),’k’,’linewidth’,1.5) 422 hold on 423 xlabel(’$m$’,’Interpreter’,’latex’) 424 ylabel(’$\phi$’,’Interpreter’,’latex’) 425 title(’Distribuci\’’on de $\phi$’,’Interpreter’,’latex’) 426 legend({’Central Differencing Scheme’,’First Order Upwind’,’Second Order Upwind’,’Power- Law Scheme’,’Soluci\’’on Exacta’},’Interpreter’,’latex’, ... 427 ’Location’,’SouthEast’) A.2 Difusión no estacionaria 1D 1% PROBLEMA DE DIFUSION (NO ESTACIONARIO, EXPLICIT SCHEME) 2 3clear all;close all;clc 4 5TablaErrores = []; 6 7% MALLADO 8 9L = 0.02; % Longitud Barra (m) 10 N = 20; % Número de Nodos 11 phi_ini = 200*ones(N,1); % Valores Iniciales 12 phiexact_value = []; 13 14 caras = linspace(0,L,N+1); % Nodos Equiespaciados 15 16 % caras_theta = linspace(0,pi,N+1); % Nodos Concentrados en Extremos 17 % caras = L/2*(1 - cos(caras_theta)); 18 19 % rng = [0,L*rand(1,N-1),L]; % Nodos Aleatorios 20 % caras = sort(rng);
A.2 Difusión no estacionaria 1D 111 21 22 for i = 2:length(caras) 23 nodos(i-1) = (caras(i)-caras(i-1))/2 + caras(i-1); 24 end 25 26 % DATOS 27 28 Deltat = 0.25*ones(N,1); % Time Step (s) 29 Gamma_w = 10*ones(N,1); % Conductividad Térmica (W/(m*K)) -> CTE 30 Gamma_e = 10*ones(N,1); % Conductividad Térmica (W/(m*K)) -> CTE 31 Aw = 1*ones(N,1); % Ãrea Transversal (m^2) -> CTE 32 Ae = 1*ones(N,1); % Ãrea Transversal (m^2) -> CTE 33 nw = -ones(N,1); % Vector Normal (nw=-1 cte para caras west) 34 ne = ones(N,1); % Vector Normal (ne=1 cte para caras east) 35 rho = 1*ones(N,1); % Densidad del Fluido (kg/m^3) 36 rho_c = 10*10^6*ones(N,1); 37 u = 0*ones(N,1); % Velocidad Unidireccional (m/s) 38 Fw = rho.*Aw.*nw.*u; % Flujo de Fluido West 39 Fe = rho.*Ae.*ne.*u; % Flujo de Fluido East 40 41 % COEFICIENTES 42 43 for t = 0.25:0.25:120 44 45 phiexact_value = []; 46 47 for I = 2:N-1 48 Dw(I) = Gamma_w(I)*nw(I)*Aw(I)/(nodos(I)-nodos(I-1)); 49 De(I) = Gamma_e(I)*ne(I)*Ae(I)/(nodos(I+1)-nodos(I)); 50 DeltaV(I) = Ae(I)*(caras(I+1)-caras(I)); 51 aW(I) = 0; % Coeficiente phi_W 52 aP(I) = rho_c(I)*DeltaV(I)/Deltat(I); % Coeficiente phi_P 53 aE(I) = 0; % Coeficiente phi_E 54 aSP(I) = De(I)*phi_ini(I+1)-Dw(I)*phi_ini(I-1)+(Dw(I)-De(I)+rho_c(I)*DeltaV(I)/ Deltat(I))*phi_ini(I); 55 A(I,I-1) = aW(I); 56 A(I,I) = aP(I); 57 A(I,I+1) = aE(I); 58 b(I) = aSP(I); 59 end 60 61 % CONDICIONES DE CONTORNO 62 63 % LEFT 64 65 % qI = 0; % CdC Izquierda: flujo impuesto (K); 66 Dw(1) = Gamma_w(1)*nw(1)*Aw(1)/(nodos(1)-caras(1)); 67 De(1) = Gamma_e(1)*ne(1)*Ae(1)/(nodos(2)-nodos(1)); 68 DeltaV(1) = Ae(1)*(caras(2)-caras(1)); 69 aW(1) = 0; 70 aP(1) = rho_c(1)*DeltaV(1)/Deltat(1); 71 aE(1) = 0; 72 aSP(1) = De(1)*phi_ini(2)+(-De(1)+rho_c(1)*DeltaV(1)/Deltat(1))*phi_ini(1); 73 A(1,1) = aP(1); 74 A(1,2) = aE(1); 75 b(1) = aSP(1); 76 77 % RIGHT 78 79 phi_R = 0; % CdC Derecha: phi impuesto (K); 80 Dw(N) = Gamma_w(N)*nw(N)*Aw(N)/(nodos(N)-nodos(N-1)); 81 De(N) = Gamma_e(N)*ne(N)*Ae(N)/(caras(end)-nodos(N));
118 Capítulo A. Códigos de MATLAB 147 aSP(k) = 0; 148 A(k,k) = aP(k); 149 A(k,k+1) = aE(k); 150 A(k,k+M) = aN(k); 151 A(k,k-M) = aS(k); 152 b(k) = aSP(k); 153 end 154 end 155 156 % EAST 157 158 for j = 2:N-1 159 for i = M 160 k = (j-1)*M+i; 161 aW(k) = Dw; 162 aP(k) = Dn-Dw-Ds; 163 aE(k) = 0; 164 aS(k) = Ds; 165 aN(k) = -Dn; 166 aSP(k) = qE*Ae*ne; 167 A(k,k-1) = aW(k); 168 A(k,k) = aP(k); 169 A(k,k+M) = aN(k); 170 A(k,k-M) = aS(k); 171 b(k) = aSP(k); 172 end 173 end 174 175 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 176 177 % SOUTH-WEST CORNER 178 179 for j=1 180 for i = 1 181 k = (j-1)*M+i; 182 aW(k) = 0; 183 aP(k) = De+Dn-2*Ds; 184 aE(k) = -De; 185 aS(k) = 0; 186 aN(k) = -Dn; 187 aSP(k) = -2*Ds*TS; 188 A(k,k) = aP(k); 189 A(k,k+1) = aE(k); 190 A(k,k+M) = aN(k); 191 b(k) = aSP(k); 192 end 193 end 194 195 % SOUTH-EAST CORNER 196 197 for j = 1 198 for i = M 199 k = (j-1)*M+i; 200 aW(k) = Dw; 201 aP(k) = Dn-Dw-2*Ds; 202 aE(k) = 0; 203 aS(k) = 0; 204 aN(k) = -Dn; 205 aSP(k) = -2*Ds*TS+qE*Ae*ne; 206 A(k,k-1) = aW(k); 207 A(k,k) = aP(k); 208 A(k,k+M) = aN(k);
A.3 Difusión estacionaria 2D con mallado estructurado 119 209 b(k) = aSP(k); 210 end 211 end 212 213 % NORTH-WEST CORNER 214 215 for j = N 216 for i = 1 217 k = (j-1)*M+i; 218 aW(k) = 0; 219 aP(k) = De-Ds; 220 aE(k) = -De; 221 aS(k) = Ds; 222 aN(k) = 0; 223 aSP(k) = +qN*An*nn; 224 A(k,k) = aP(k); 225 A(k,k+1) = aE(k); 226 A(k,k-M) = aS(k); 227 b(k) = aSP(k); 228 end 229 end 230 231 % NORTH-EAST CORNER 232 233 for j = N 234 for i = M 235 k = (j-1)*M+i; 236 aW(k) = Dw; 237 aP(k) = -Dw-Ds; 238 aE(k) = 0; 239 aS(k) = Ds; 240 aN(k) = 0; 241 aSP(k) = qE*Ae*ne+qN*An*nn; 242 A(k,k-1) = aW(k); 243 A(k,k) = aP(k); 244 A(k,k-M) = aS(k); 245 b(k) = aSP(k); 246 end 247 end 248 249 % SOLUCION 250 251 phi1D = A\b’; 252 253 for j = 1:N 254 for i = 1:M 255 k = (j-1)*M+i; 256 phi2D(i,j) = phi1D(k); 257 xiplot(i,j) = nodos(k,1); 258 etaplot(i,j) = nodos(k,2); 259 end 260 end 261 phi2D 262 263 figure(3) 264 h=surf(xiplot,etaplot,phi2D); 265 colormap(turbo) 266 cb=colorbar; 267 shading interp 268 set(h,’EdgeColor’,’k’) 269 axis([0 L 0 H]) 270 xlabel(’$x\,(m)$’,’Interpreter’,’latex’)
120 Capítulo A. Códigos de MATLAB 271 ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 272 zlabel(’$T\,(^{\circ}C)$’,’Interpreter’,’latex’) 273 title(’Campo de Temperaturas’,’Interpreter’,’latex’) 274 275 figure(4) 276 contour(xiplot,etaplot,phi2D,150:10:270,’Linewidth’,2,’ShowText’,’on’); 277 colormap(turbo) 278 xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 279 ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 280 zlabel(’$T\,(^{\circ}C)$’,’Interpreter’,’latex’) 281 title(’Curvas de nivel del Campo de Temperaturas’,’Interpreter’,’latex’) A.4 Convección-difusión estacionaria 2D con mallado estructurado 1%% CONVECCION-DIFUSION 2D 2 3clear all;clc; 4 5% MALLADO EN X-Y 6 7M = 50; 8N = 10; 9L = 0.5; 10 H = 0.1; 11 Deltax = L/M; 12 Deltay = H/N; 13 14 for j = 1:(N+1) 15 ximat(1:(M+1),j) = linspace(0,L,(M+1)); 16 etamat(1:(M+1),j) = H*(j-1)/N; 17 end 18 19 % figure (1) 20 % 21 % for j = 1:(N+1) 22 % plot(ximat(1:(M+1),j), etamat(1:(M+1),j),’r’,’linewidth’,1.5) 23 % hold on 24 % end 25 % for i = 1:(M+1) 26 % plot(ximat(i,1:(N+1)), etamat(i,1:(N+1)),’r’,’linewidth’,1.5) 27 % hold on 28 % end 29 % xlabel(’$x\,(m)$’,’Interpreter’,’Latex’) 30 % ylabel(’$y\,(m)$’,’Interpreter’,’Latex’) 31 % title(’Mallado del dominio’,’Interpreter’,’Latex’) 32 % axis([-0.1 0.6 -0.1 0.2]) 33 34 % COORDENADAS DE LOS VERTICES Y NODOS 35 36 for j = 1:N 37 for i = 1:M 38 k = (j-1)*M+i; 39 xi1(k) = ximat(i,j); eta1(k) = etamat(i,j); 40 xi2(k) = ximat(i,j+1); eta2(k) = etamat(i,j+1); 41 xi3(k) = ximat(i+1,j+1); eta3(k) = etamat(i+1,j+1); 42 xi4(k) = ximat(i+1,j); eta4(k) = etamat(i+1,j); 43 xinodo(k) = xi1(k)+(xi3(k)-xi1(k))/2; 44 etanodo(k) = eta1(k)+(eta2(k)-eta1(k))/2; 45
A.4 Convección-difusión estacionaria 2D con mallado estructurado 121 46 end 47 end 48 49 % figure(1) 50 nodos = [xinodo’,etanodo’]; 51 % hold on 52 % plot(nodos(:,1),nodos(:,2),’.’,’MarkerSize’,25) 53 54 % DATOS 55 56 Np = M*N; 57 58 Gamma = 0.613; % Conductividad Térmica (W/(m*K)) -> CTE 59 nw = -1; % Vector Normal (nw=-1 cte para caras west) 60 ne = 1; % Vector Normal (ne=1 cte para caras east) 61 ns = -1; % Vector Normal (ns=-1 cte para caras south) 62 nn = 1; % Vector Normal (nn=1 cte para caras north) 63 rho = 1000; % Densidad del Fluido (kg/m^3) 64 cp = 4200; % Calor Específico del Agua (J/(kg*K)) 65 t = 1; 66 um = 2.92*10^-7; % Velocidad media (m/s) 67 Dw = Gamma*nw*Deltay*t/Deltax; 68 De = Gamma*ne*Deltay*t/Deltax; 69 Ds = Gamma*ns*Deltax*t/Deltay; 70 Dn = Gamma*nn*Deltax*t/Deltay; 71 qN = 15; 72 qE = -20; 73 TW = 303; 74 k = [1:1:Np]’; 75 76 % NODOS INTERMEDIOS 77 78 79 for j = 2:N-1 80 for i = 2:M-1 81 k = (j-1)*M+i; 82 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 83 Fe(k) = rho*cp*u(k)*ne*Deltay; 84 Fw(k) = rho*cp*u(k)*nw*Deltay; 85 aW(k) = Dw+Fw(k); 86 aP(k) = Fe(k)+De+Dn-Dw-Ds; 87 aE(k) = -De; 88 aS(k) = Ds; 89 aN(k) = -Dn; 90 aSP(k) = 0; 91 A(k,k-1) = aW(k); 92 A(k,k) = aP(k); 93 A(k,k+1) = aE(k); 94 A(k,k+M) = aN(k); 95 A(k,k-M) = aS(k); 96 b(k) = aSP(k); 97 end 98 end 99 100 %%% CONDICIONES DE CONTORNO %%% 101 102 % SOUTH 103 104 for j = 1 105 for i = 2:M-1 106 k = (j-1)*M+i; 107 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2);
122 Capítulo A. Códigos de MATLAB 108 Fe(k) = rho*cp*u(k)*ne*Deltay; 109 Fw(k) = rho*cp*u(k)*nw*Deltay; 110 aW(k) = Fw(k)+Dw; 111 aP(k) = Fe(k)+De+Dn-Dw; 112 aE(k) = -De; 113 aS(k) = 0; 114 aN(k) = -Dn; 115 aSP(k) = 0; 116 A(k,k-1) = aW(k); 117 A(k,k) = aP(k); 118 A(k,k+1) = aE(k); 119 A(k,k+M) = aN(k); 120 b(k) = aSP(k); 121 end 122 end 123 124 % NORTH 125 126 for j = N 127 for i = 2:M-1 128 k = (j-1)*M+i; 129 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 130 Fe(k) = rho*cp*u(k)*ne*Deltay; 131 Fw(k) = rho*cp*u(k)*nw*Deltay; 132 aW(k) = Fw(k)+Dw; 133 aP(k) = Fe(k)+De-Dw-Ds; 134 aE(k) = -De; 135 aN(k) = 0; 136 aS(k) = Ds; 137 aSP(k) = qN*Deltax*nn; 138 A(k,k-1) = aW(k); 139 A(k,k) = aP(k); 140 A(k,k+1) = aE(k); 141 A(k,k-M) = aS(k); 142 b(k) = aSP(k); 143 end 144 end 145 146 % WEST 147 148 for j = 2:N-1 149 for i = 1 150 k = (j-1)*M+i; 151 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 152 Fe(k) = rho*cp*u(k)*ne*Deltay; 153 Fw(k) = rho*cp*u(k)*nw*Deltay; 154 aW(k) = 0; 155 aP(k) = Fe(k)+De+Dn-2*Dw-Ds; 156 aE(k) = -De; 157 aS(k) = Ds; 158 aN(k) = -Dn; 159 aSP(k) = -(2*Dw+Fw(k))*TW; 160 A(k,k) = aP(k); 161 A(k,k+1) = aE(k); 162 A(k,k+M) = aN(k); 163 A(k,k-M) = aS(k); 164 b(k) = aSP(k); 165 end 166 end 167 168 % EAST 169
A.4 Convección-difusión estacionaria 2D con mallado estructurado 123 170 for j = 2:N-1 171 for i = M 172 k = (j-1)*M+i; 173 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 174 Fe(k) = rho*cp*u(k)*ne*Deltay; 175 Fw(k) = rho*cp*u(k)*nw*Deltay; 176 aW(k) = Fw(k)+Dw; 177 aP(k) = Fe(k)+Dn-Dw-Ds; 178 aE(k) = 0; 179 aS(k) = Ds; 180 aN(k) = -Dn; 181 aSP(k) = qE*Deltay*ne; 182 A(k,k-1) = aW(k); 183 A(k,k) = aP(k); 184 A(k,k+M) = aN(k); 185 A(k,k-M) = aS(k); 186 b(k) = aSP(k); 187 end 188 end 189 190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 191 192 % SOUTH-WEST CORNER 193 194 for j = 1 195 for i = 1 196 k = (j-1)*M+i; 197 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 198 Fe(k) = rho*cp*u(k)*ne*Deltay; 199 Fw(k) = rho*cp*u(k)*nw*Deltay; 200 aW(k) = 0; 201 aP(k) = Fe(k)+De+Dn-2*Dw; 202 aE(k) = -De; 203 aS(k) = 0; 204 aN(k) = -Dn; 205 aSP(k) = -(2*Dw+Fw(k))*TW; 206 A(k,k) = aP(k); 207 A(k,k+1) = aE(k); 208 A(k,k+M) = aN(k); 209 b(k) = aSP(k); 210 end 211 end 212 213 % SOUTH-EAST CORNER 214 215 for j = 1 216 for i = M 217 k = (j-1)*M+i; 218 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 219 Fe(k) = rho*cp*u(k)*ne*Deltay; 220 Fw(k) = rho*cp*u(k)*nw*Deltay; 221 aW(k) = Fw(k)+Dw; 222 aP(k) = Fe(k)+Dn-Dw; 223 aE(k) = 0; 224 aS(k) = 0; 225 aN(k) = -Dn; 226 aSP(k) = qE*Deltay*ne; 227 A(k,k-1) = aW(k); 228 A(k,k) = aP(k); 229 A(k,k+M) = aN(k); 230 b(k) = aSP(k); 231 end
124 Capítulo A. Códigos de MATLAB 232 end 233 234 % NORTH-WEST CORNER 235 236 for j=N 237 for i = 1 238 k = (j-1)*M+i; 239 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 240 Fe(k) = rho*cp*u(k)*ne*Deltay; 241 Fw(k) = rho*cp*u(k)*nw*Deltay; 242 aW(k) = 0; 243 aP(k) = Fe(k)+De-2*Dw-Ds; 244 aE(k) = -De; 245 aS(k) = Ds; 246 aN(k) = 0; 247 aSP(k) = -(2*Dw+Fw(k))*TW+qN*Deltax*nn; 248 A(k,k) = aP(k); 249 A(k,k+1) = aE(k); 250 A(k,k-M) = aS(k); 251 b(k) = aSP(k); 252 end 253 end 254 255 % NORTH-EAST CORNER 256 257 for j=N 258 for i = M 259 k = (j-1)*M+i; 260 u(k) = 6*um*(etanodo(k)/H-(etanodo(k)/H)^2); 261 Fe(k) = rho*cp*u(k)*ne*Deltay; 262 Fw(k) = rho*cp*u(k)*nw*Deltay; 263 aW(k) = Fw(k)+Dw; 264 aP(k) = Fe(k)-Dw-Ds; 265 aE(k) = 0; 266 aS(k) = Ds; 267 aN(k) = 0; 268 aSP(k) = qN*Deltax*nn+qE*Deltay*ne; 269 A(k,k-1) = aW(k); 270 A(k,k) = aP(k); 271 A(k,k-M) = aS(k); 272 b(k) = aSP(k); 273 end 274 end 275 276 % SOLUCION 277 278 phi1D = A\b’-273; 279 280 for j = 1:N 281 for i = 1:M 282 k = (j-1)*M+i; 283 phi2D(i,j) = phi1D(k); 284 xiplot(i,j) = nodos(k,1); 285 etaplot(i,j) = nodos(k,2); 286 end 287 end 288 phi2D 289 290 figure(4) 291 h=surf(xiplot,etaplot,phi2D); 292 colormap(turbo) 293 cb=colorbar;
A.4 Convección-difusión estacionaria 2D con mallado estructurado 125 294 shading interp 295 set(h,’EdgeColor’,’k’) 296 xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 297 ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 298 zlabel(’$T\,(^{\circ}C)$’,’Interpreter’,’latex’) 299 title(’Campo de Temperaturas’,’Interpreter’,’latex’) 300 301 figure(5) 302 contour(xiplot,etaplot,phi2D,[30:1:45],’Linewidth’,2,’ShowText’,’on’); 303 colormap(turbo) 304 xlabel(’$x\,(m)$’,’Interpreter’,’latex’) 305 ylabel(’$y\,(m)$’,’Interpreter’,’latex’) 306 zlabel(’$T\,(^{\circ}C)$’,’Interpreter’,’latex’) 307 title(’Curvas de nivel del Campo de Temperaturas’,’Interpreter’,’latex’)
Bibliografía [1] S. V. Patankar (1980): Numerical Heat Transfer and Fluid Flow, 1ªEdición, Boca Raton (Florida), CRC Press. [2] A. Barrero y M. Pérez-Saborid (2005): Fundamentos y Aplicaciones de la Mecánica De Fluidos, 1ªEdición, España, McGraw Hill. [3] A. W. Date (2005): Introduction to Computational Fluid Dynamics, Instituto Tecnológico de India (Bombay), Cambridge University Press. [4] H. K. Versteeg and W. Malalasekera (2007): An Introduction to Computational Fluid Dynamics, 2ªEdición, Inglaterra, Pearson Education. [5] V. R. Voller (2009): Basic Control Volume Finite Element Methods for Fluids and Solids, 1ª Edición, Singapore, IISc Press. [6] S. Mazumder (2016): Numerical Methods for Partial Differential Equations, Universidad Estatal de Ohio, Elsevier Inc. [7] ANSYS (2017): ANSYS Fluent Theory Guide, Release 18.2, ANSYS, Inc. and ANSYS Europe. [8] A. G. Galindo (2017): Modelado CFD de la atomización del combustible en quemadores de turbinas de gas, Trabajo Fin de Máster, Universitat Politécnica de València, Escuela Técnica Superior de Ingeniería del Diseño. Disponible en: https:// riunet.upv.es/ handle/ 10251/ 88196 [Consulta: 25-03-2022]. [9] A. Bakker (2008): Applied Computational Fluid Dynamics. CFD Lectures, Lecture 5-Solution Methods. Disponible en: http:// www.bakker.org [Consulta: 31-03-2022]. [10] A. Kaushik (2016): Critical Evaluation of Four Differencing Schemes for A Steady Convection- Diffusion Problem. Research Journal of Mathematical and Statistical Sciences, Volume 4, No. 2, pp 6-12 Disponible en: http:// www.isca.in/ MATH_SCI/ v4i2.php [Consulta: 12-04-2022]. [11] Y. Niño (2008): Método de Volúmenes Finitos. Dpto de Ing. Civil, Universidad de Chile. Disponible en: https:// www.u-cursos.cl/ ingenieria/ 2008/ 2/ CI71D/ 1/ material_docente/ detalle? id=200975 [Consulta: 20-04-2022]. [12] A. Ateş, O. Altun and A. Kilicman (2017): On A Comparison of Numerical Solution Methods for General Transport Equation on Cylindrical Coordinates. Applied Mathematics & Information Sciences, Volume 11, No. 2, pp 433-439. Disponible en: https:// www.naturalspublishing. com/ Article.asp?ArtcID=12846 [Consulta: 20-04-2022]. 127