Análisis de la influencia de la radiación térmica en los procesos de convección de Rayleigh-Benard

Equation Chapter 1 Section 1 Trabajo de Fin de Grado Ingeniería Aeroespacial Análisis de la influencia de la radiación térmica en los procesos de convección de Rayleigh-Benard Autor: Ángel Enrique Boyer Varela Tutor: Miguel Pérez-Saborid Sánchez-Pastor Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 iii Trabajo de Fin de Grado Ingeniería Aeroespacial Análisis de la influencia de la radiación térmica en los procesos de convección de Rayleigh-Benard Autor: Ángel Enrique Boyer Varela Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor titular Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Índice general 1. Introducción 7 1.1. Objetivos y descripción del proyecto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2. Concepto de convección. Convección forzada y natural . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3. Convecciónenlaatmósfera ....................................... 9 1.4. Transferencia de calor en el sol. Zonas de radiación y convectiva . . . . . . . . . . . . . . . . . . . 9 1.5. Estructuradeltrabajo.......................................... 11 2. Formulación del problema 13 2.1. Mecanismos físicos de la radiación y medio participante . . . . . . . . . . . . . . . . . . . . . . . 13 2.2. Ecuaciones del problema de convección con ujo de radiación. Aproximaciones . . . . . . . . . . 14 2.2.1. Ecuaciones del problema de convección con radiación. Camino libre medio . . . . . . . . . 14 2.2.2. EcuacionesdeKourgano.................................... 16 2.2.3. Aproximación ópticamente na o transparente . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.4. Aproximación ópticamente gruesa u opaca . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3. La convección de Rayleigh-Bénard. Aproximación de Boussinesq . . . . . . . . . . . . . . . . . . 22 2.4. Resultados principales del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3. Método numérico 25 3.1. Método de colocación en problemas 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2. Resoluciónnumérica........................................... 28 3.2.1. Condiciones de contorno e inicial. Problema con paredes no participativas . . . . . . . . . 29 3.2.2. Condiciones de contorno e inicial en el problema con paredes participativas en la radiación 35 4. Medio participativo dentro de una cavidad con paredes con absorción de radiación nula 45 4.1. Paredes verticales adiabáticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 1 2 Índice general 4.1.1. Número de Rayleigh crítico: resultados fundamentales . . . . . . . . . . . . . . . . . . . . 45 4.1.2. Variación del no de Rayleigh crítico con la relación de aspectos. Inuencia del medio . . . 47 4.1.3. Número de Nusselt y perles de temperatura . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2. Parbaroclínico .............................................. 50 4.3. Paredes horizontales adiabáticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3.1. Isocontornos de temperatura y función de corriente . . . . . . . . . . . . . . . . . . . . . . 51 4.4. Resolucióndelcasoopaco ........................................ 52 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula 55 5.1. Paredes con absorción de radiación no nula. Resolución numérica del problema para una cavidad rectangular con paredes horizontales adiabáticas . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.1.1. Perles de velocidades e isocontornos de temperatura . . . . . . . . . . . . . . . . . . . . 56 5.1.2. Número de Nusselt convectivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.1.3. Temperatura en las paredes aisladas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2. Corrección de los resultados: condiciones de contorno sin linealizar . . . . . . . . . . . . . . . . . 62 5.2.1. Códigomejorado......................................... 64 5.2.2. Temperatura de equilibrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2.3. Resultados corregidos para emisividades elevadas . . . . . . . . . . . . . . . . . . . . . . . 71 5.3. Paredes con absorción de radiación no nula. Resolución numérica del problema para una cavidad rectangular con paredes verticales adiabáticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.3.1. Isocontornos del campo de temperaturas y de la función de corriente . . . . . . . . . . . . 73 5.3.2. Número de Nusselt en función de Rayleigh . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.3.3. Perl de temperaturas en las paredes aisladas . . . . . . . . . . . . . . . . . . . . . . . . . 75 6. Conclusiones y líneas de desarrollo 77 2 Índice de guras 1.1. Esquema de un sistema concentrador de energía solar. . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2. Capassolares............................................... 10 1.3. Celdas convectivas en la supercie del sol, imagen cortesía de NASA. . . . . . . . . . . . . . . . . 10 2.1. Balance de calor en una pared lateral de la cavidad (aislada). . . . . . . . . . . . . . . . . . . . . 14 2.2. Esquema del emisor, medio y supercie detectora. . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3. Solución de la temperatura de equilibrio θe con parámetro de radiación-conducción Rc nulo para 3 valores del parámetro λ ......................................... 19 2.4. Representación del paralelepípedo que contiene el uido y de las condiciones de contorno del problema.................................................. 23 3.1. Ejemplo de función aproximada por interpolación usando nodos equiespaciados (a la izquierda) y nodos de Chebyshev (a la derecha). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2. Representación esquemática explicativa sobre la división de los contornos horizontales y verticales enelementosdepared........................................... 36 3.3. Factor de forma entre supercies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4. Determinación del factor de forma entre dos supercies: longitudes de cuerdas cruzadas ( L5y L6 ), no cruzadas ( L3y L4 ) y de las supercies ( L1y L2 ). ......................... 37 4.1. Gráca tomada como referencia para elegir los valores de λG=λ2 , donde λG es el parámetro queapareceenelejehorizontal...................................... 46 4.2. Número de Nusselt medio medido en la pared caliente, z=0, para el caso de convección natural, línea azul, ya resuelto en [1] y para los valores de λG= 1 , λG= 5 y λG= 10 . ........... 49 4.3. Perles de la temperatura y la de equilibrio, adimensionalizadas con la temperatura media, para varios valores de λ , jando el parámetro Ra/Rac= 2 para cada caso. . . . . . . . . . . . . . . . . 50 4.4. Par baroclínico, gura cortesía de [6] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.5. Isocontornos de los campos de temperatura y función de corriente. . . . . . . . . . . . . . . . . . 52 4.6. Ψ∗ e isocontornos de temperatura para paredes libres en el caso de aproximación opaca. . . . . . 53 5.1. Campo de velocidades e isocontornos de temperatura. . . . . . . . . . . . . . . . . . . . . . . . . 57 3 4 Índice de guras 5.2. Campo de velocidades e isocontornos de temperatura para números de Rayleigh bajos. . . . . . . 58 5.3. Tabla calculada frente a gráca dada por [24] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.4. Isocontorno de temperaturas para valores altos de la emisividad. . . . . . . . . . . . . . . . . . . 60 5.5. Perles de temperatura a lo largo de las paredes aisladas. . . . . . . . . . . . . . . . . . . . . . . 61 5.6. Temperatura de equilibrio con paredes emisoras. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.7. Comparativa del número de Nusselt medido en ambas paredes verticales con la gráca obtenido por Akiyama, para = 1 . ........................................ 71 5.8. Isocontornos de temperaturas para emisividad unidad. . . . . . . . . . . . . . . . . . . . . . . . . 72 5.9. Perl de temperaturas de las paredes aisladas para emisividad unidad. . . . . . . . . . . . . . . . 72 5.10.Comparativacon[23]........................................... 73 5.11. Isocontornos de la función de corriente (izquierda) y del campo de temperaturas (derecha). . . . 74 5.12. Nusselt medio en pared fría y caliente en función de Rayeligh. . . . . . . . . . . . . . . . . . . . . 75 5.13. Perl de temperaturas para paredes con emisividad nula y unidad. . . . . . . . . . . . . . . . . . 76 4 Índice de cuadros 4.1. Tabla comparativa con los resultados de Goody de la gura 4.1. . . . . . . . . . . . . . . . . . . . 46 4.2. Tabla comparativa del número de Rayleigh crítico con el caso de convección natural si se tiene en cuenta aproximación de medio transparente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3. Número de Rayleigh crítico en función de A para varios valores del parámetro λG . ........ 47 5 12 Capítulo 1. Introducción 12 Capítulo 2 Formulación del problema 2.1. Mecanismos físicos de la radiación y medio participante Todos los cuerpos emiten radiación en todas las direcciones, con distinta intensidad, a su alrededor a través de ondas electromagnéticas (fotones) debido a la conversión de la energía interna del cuerpo en radiación, por la agitación molecular y atómica que tiene asociadas. La radiación térmica es una forma de emisión electromagnética y, por tanto, puede propagarse por el vacío. Ejemplos usuales de radiación son la que llega a la Tierra proveniente del sol o la disipación de calor desde un objeto incandescente. Así, el calor es transmitido entre objetos distanciados con diferente temperatura. El intercambio de radiación depende de la temperatura y del estado de la supercie emisora. En el caso de líquidos y sólidos solo una na capa participa en la radiación, sobre todo en los sólidos donde el fenómeno puede considerarse totalmente supercial. Aunque, realmente, también depende del grosor de la capa y de la presión en el caso de cuerpos semitransparentes como el vidrio. En el caso de gases, la emisión y absorción de radiación pueden considerarse efectos volumétricos. Cuando la radiación llega a un cuerpo, este puede absorber parte, reejar parte y dejar pasar el resto. La parte que se absorbe es, lógicamente, la que se transforma en calor. La fracción (es decir la parte entre el total) de cada una de las partes se denominan, respectivamente, absortividad, reectividad y transmisividad. Si un cuerpo tiene una transmisividad igual a la unidad, su reectividad y absortividad serían nulas y la radiación incidente atravesaría por completo al cuerpo sin absorberse ni reejarse nada. Esto es lo que se conoce como cuerpo transparente . Por ejemplo, el aire es un gas que se puede considerar casi transparente, lo que no sucede con gases poliatómicos como el dióxido de carbono CO2 o el metano CH4 , que tienen capacidad de absorber parte de la radiación incidente. En contraste, un cuerpo cuya transmisividad sea nula, es un cuerpo opaco y por tanto el total de la radiación que le llega puede ser parcialmente reejada y parcialmente absorbida. En denitiva, un cuerpo de este tipo, si tiene buena reectividad, tendrá baja absortividad y viceversa. Además, se dene la dispersión ( scattering en inglés) como el cambio de dirección con respecto a la de propagación de la radiación sin pérdida de energía. Por otro lado, se dene de forma necesaria una magnitud llamada radiosidad J que representa a a la cantidad de energía de radiación que abandona un cuerpo por unidad de área y tiempo teniendo en cuenta todas las direcciones; es decir, esta magnitud tiene en cuenta la suma de la intensidad emitida y reejada. Una vez denidos estos conceptos, se procede a explicar qué se entiende por medios participativos y cuáles son sus propiedades. Considérese una intensidad Iν que se propaga a través de un medio en una dirección particular. Cuando atraviesa un elemento de espesor dS , se atenúa por la absorción y dispersión debidas a la inuencia del medio. La pérdida de intensidad en dicho elemento viene dada por dIν, atenuacion =−βνIν, atenuacion(S)dS, (2.1) donde βν es el coeciente de atenuación, que se puede demostrar que coincide con el inverso del camino libre medio y que es la suma de dos coecientes: βν=κν+σsν , con unidades de longitud inversa m−1 . El primero κν (véase [8]) corresponde al coeciente de absorción y el segundo σsν al coeciente de dispersión. Un medio que es capaz de inuir en la transferencia de radiación de esta manera, dispersándola o absorbiéndola, es lo que se conoce como un medio participativo (o participating medium en inglés) [13]. 13 14 Capítulo 2. Formulación del problema Se quiere hacer hincapié en que el estudio de la radiación es altamente complejo, y para poder obtener modelos que reproduzcan resultados, siempre es necesario realizar una serie de hipótesis simplicativas que permitan obtener soluciones, tanto para reducir la dicultad de las ecuaciones involucradas como las condiciones de contorno del problema. En el caso de este trabajo, al tenerse como condición de contorno paredes aisladas en la cavidad en el caso de que dicha pared no tenga temperatura impuesta (o lo que es lo mismo, es desconocida a priori), habría que tener en cuenta que en dichas paredes, la suma del calor que se reciba o emita por convección más el calor que se reciba o emita por radiación, sea nula (ver gura 2.1). Esta condición de contorno, denominada condición de contorno de convección-radiación, surge de la suposición de que la cavidad en la que está encerrado el uido está en contacto (por fuera de la cavidad) con un material aislante. Así, en esta situación no aparecen términos de convección en la condición de contorno, con la consecuencia de que habría que modelar de alguna manera ese calor intercambiado en forma de radiación que aparece en la gura mencionada como qr y está representado con echas no rectas (el término de calor por convección entre el udio y la pared sería simplemente −K∂T/∂z o −K∂T/∂x , dependiendo de si la pared aislada es horizontal o vertical respectivamente). En otros trabajos, como The inuence of radiative transfer on cellular convection o en trabajos donde se ha estudiado la inuencia de la radiación en procesos de convección en un uido encerrado entre dos esferas concéntricas, este término ni si quiera se plantea, pues la cavidad es innita en el primer caso y en el segundo no hay paredes laterales (por ejemplo eso sucedería en una atmósfera estelar). De este modo, el objetivo es al menos poder ilustrar de forma cualitativa el efecto que la radiación tiene en los procesos de convección de Rayleigh-Bénard y, en cualquier caso, debido a la enorme dicultad que supone siempre el estudio de la radiación, en cualquier documento en el que se estudien procesos en los que tenga relevancia, nunca se encuentran exentos de hipótesis que permitan obtener resultados que no dejan de ser aproximaciones, con menor o mayor grado de exactitud. Figura 2.1: Balance de calor en una pared lateral de la cavidad (aislada). El color que recubre a la cavidad está representando al material aislante mencionado. 2.2. Ecuaciones del problema de convección con ujo de radiación. Aproximaciones 2.2.1. Ecuaciones del problema de convección con radiación. Camino libre medio Supóngase que el paralelepípedo de la gura 2.4 recibe un ujo de calor por radiación RH , del inglés radiative heating , en su pared horizontal inferior z= 0 y que el uido dentro de él es capaz de absorber y emitir radiación térmica. Una solución completa para este problema sería extremadamente compleja, por lo que se van a utilizar las aproximaciones, una apropiada para medios opacos y otra para medios transparentes. En esta sección, se va a exponer la ecuación de la energía completa, teniendo en cuenta el calor que se transmite por conducción y radiación. De este modo, se tiene: ∂T ∂t +vx ∂T ∂x +vy ∂T ∂y +vz ∂T ∂z =α∇2T+RH ρ0cp , (2.2) donde α=K/(ρ0cp) , es la difusividad térmica. 14 2.2. Ecuaciones del problema de convección con ujo de radiación. Aproximaciones 15 La complejidad del asunto reside en que es necesario modelar el término RH para poder trabajar con el nuevo sistema. Para ello, se recurre en la sección próxima al uso de dos aproximaciones relacionadas con el valor del coeciente de absorción κν , o mejor dicho, con el valor inverso de dicho coeciente, denominado camino libre medio y comparado con la dimensión H de la cavidad en la que se mueve el uido. Entiéndase dicha magnitud κ−1 como la distancia que recorren las partículas de una onda electromagnética, en media, entre colisiones con las partículas del medio que atraviesa. Dicha magnitud puede ser muy pequeña (las partículas colisionarían con una alta frecuencia al atravesar un medio) o muy grande (las partículas tienen prácticamente vía libre e incluso una alta fracción de ellas podría llegar a atravesar el medio sin colisionar). Ambas situaciones, que se desarrollarán en las siguientes secciones, se pueden encontrar en medios como la atmósfera del Sol . Por ejemplo, en la zona de radiación y en la base de la zona de convección en el sol, los fotones rebotan en incontables ocasiones antes de atravesar dicha zona, por lo que se puede entender que en esta zona el medio es ópticamente grueso. Sin embargo, si se avanza hacia la fotosfera del sol, tan solo unos cientos de kilómetros dentro de ella, el medio pasa de opaco a transparente [14]. Una vez se ha denido el coeciente de absorción κν , se hace necesario explicar qué se entiende por profundidad óptica ( optical depth/thickness en inglés). Este parámetro describe cuánta absorción tiene lugar cuando la luz atraviesa un medio absorbente, por ejemplo, la atmósfera solar. Considérese como un haz de fotones la luz emitida por un emisor hacia una supercie, algunos de los cuales pueden ser absorbidos por el medio que atraviesan para llegar hasta ella. La probabilidad innitesimal dpν de un fotón emitido a una frecuencia ν en una rebanada de espesor ds (ver gura 2.2), es directamente proporcional a dicho espesor: dpν=κνds ⇒κν=dpν ds . (2.3) Figura 2.2: Esquema del emisor, medio y supercie detectora. Imagen extraída de: https://www.cv.nrao.edu/course/astr534/Radxfer.html Entonces, esto signica que la probabilidad de absorción no es constante. Para ilustrar esto, téngase en mente que a medida que la luz viaja, algunos fotones son absorbidos, de tal modo que a mayor distancia del emisor, menor será la cantidad de fotones restantes, con lo que la probabilidad de absorberlos aumentará de forma no lineal con el grosor de la rebanada. La fracción de la intensidad perdida debida solo a la absorción durante un desplazamiento innitesimal (ver ecuación 2.1) viene dada por: dIν Iν =−dpν=−κν. (2.4) Integrando a ambos lados de la ecuación a lo largo de todo el camino, se obtiene la cantidad de intensidad al nal de este, en función de la intensidad inicial. Iνf =Iν0e−Rsf 0κν(s0)ds, (2.5) 15 16 Capítulo 2. Formulación del problema donde la cantidad del exponente τν , τν=−Zsf 0 κν(s0)ds0, (2.6) se denomina profundidad óptica . Si τν<< 1 , la intensidad a la frecuencia ν permanece prácticamente constante, y se puede decir que el medio es ópticamente no o transparente , mientras que si τν>> 1 , la intensidad Iν cae rápidamente, absorbiéndose los fotones por el medio de forma inmediata y se considera que este medio es ópticamente grueso u opaco . En las siguientes secciones se desarrollarán estas ideas. 2.2.2. Ecuaciones de Kourgano La ecuación que gobierna en las variaciones de intensidad de radiación en una dirección ~s fue expuesta por Kourgano (1952) [15]: dI(~s) ds =κ[B−I(~s0)], (2.7) donde I(~s) es la intensidad de radiación in la dirección de ~s , B es la intensidad de Planck de un cuerpo negro 1 , ds es un desplazamiento innitesimal en la dirección ~s y en κ se ha omitido el subíndice. El ujo de calor por radiación RH es: RH =−Z4π dI(~s) ds dω, (2.8) donde ω representa a un elemento de ángulo sólido. Al ser isotrópica la intensidad de radiación de un cuerpo negro, se tiene, de la combinación de la ecuaciones anteriores: RH =−4πκB +κZI(~s)dω. (2.9) El primer término de 2.9 representa pérdidas de calor por emisión de radiación térmica en el punto del uido a la temperatura local, ya que B depende de la temperatura según la ley de Stefan, ecuación 2.10: B=σ πT4, (2.10) donde σ= 5.67 ·10−8W m2·K4 es la constante de Stefan-Boltzmann. El segundo término de 2.9 representa el calor absorbido en el punto del uido y emitido desde otros puntos del uido y su contorno. Se puede demostrar que las dimensiones de las celdas convectivas formadas en el seno de la cavidad, son del orden de H/a , siendo a un parámetro adimensional de orden unidad cuyos valores típicos se pueden consultar en Pellew & Southwell [16] y siendo H una dimensión característica de la cavidad. Por tanto, comparando los valores de el camino libre medio κ−1 y H/a se pueden obtener las aproximaciones tratadas en las próximas dos subsecciones. 2.2.3. Aproximación ópticamente na o transparente El primero de los casos abordados es aquel en el que se cumple que el orden del camino libre medio es mucho mayor que el orden del tamaño de la celda convectiva, es decir, si se cumple que κ−1>> H/a o 1 Aproximación física ideal que hace referencia a un cuerpo capaz de absorber toda la radiación electromagnética que le llega, independientemente de su frecuencia o ángulo de incidencia. 16 2.2. Ecuaciones del problema de convección con ujo de radiación. Aproximaciones 17 alternativamente κ−1H << a . En este caso, la ecuación 2.9, con el uso de la ecuación de Stefan-Boltzmann 2.10 se reduce a: RH =−4κσT4, (2.11) donde T es la temperatura. Seguidamente, es necesario introducir este término en la ecuación de la energía 2.2. Antes de ello, conviene aclarar que la diferencia de temperaturas entre las paredes inferior y superior, T1 y T2 respectivamente, varían relativamente poco, es decir, se cumple que: T1−T2 T1 << 1; T1−T2 T2 << 1. (2.12) Por tanto, el campo de temperaturas se puede modelar como una temperatura de equilibrio Te dependiente de las coordenadas x y z más una perturbación T0 debida al movimiento dependiente de las coordenadas espaciales y el tiempo. Además, esta temperatura de equilibro se puede descomponer a su vez como la suma de la temperatura media Tm=T1+T2 2 más una temperatura de perturbación de equilibrio dependiente de x y z , Te0(x, z) . En resumen: T(x, z, t) = Te(x, z) + T0(x, z, t); donde T0<< T,Te (2.13) Te(z) = Tm+T0 e(x, z); donde T0 e<< Tm,Te (2.14) El sistema de ecuaciones a resolver sería el conformado por las ecuaciones: ∂vx ∂x +∂vy ∂y +∂vz ∂z = 0; (2.15) ∂vx ∂t +vx ∂vx ∂x +vy ∂vx ∂y +vz ∂vx ∂z =−∂(pm/ρ0) ∂x +ν∇2vx; (2.16) ∂vy ∂t +vx ∂vy ∂x +vy ∂vy ∂y +vz ∂vy ∂z =−∂(pm/ρ0) ∂y +ν∇2vy; (2.17) ∂vz ∂t +vx ∂vz ∂x +vy ∂vz ∂y +vz ∂vz ∂z =−∂(pm/ρ0) ∂z +gβ(T−Tm) + ν∇2vz; (2.18) ρ0cp∂T ∂t +vx ∂T ∂x +vy ∂T ∂y +vz ∂T ∂z =K∇2T+RH, (2.19) donde RH es el término correspondiente al calor intercambiado por radiación, como se dijo en el apartado anterior. Dicho sistema 2.15-2.19 habrá que resolverlo bajo las condiciones de contorno siguientes, que se corresponden con paredes verticales aisladas y horizontales con temperaturas impuestas T1 y T2 : K∂Te ∂x x=0,L ±σ(T4 m−T4 e)=0, (2.20) T(z= 0) = T1⇒T0 e=T1−Tm=T1−T2 2, (2.21) 17 18 Capítulo 2. Formulación del problema T(z=L) = T1⇒T0 e=T2−Tm=−T1−T2 2, (2.22) donde  es la emisividad del uido. Si en la última ecuación se introducen las relaciones 2.13 y 2.14, se separa el problema de equilibrio y de movimiento, respectivamente, se obtienen las ecuaciones: K∇2T0 e−16κσT3 mT0 e= 4κσT4 m (2.23) ∂T0 ∂t +vx ∂(T0+Te) ∂x +vz ∂(T0+Te) ∂z =1 ρcpK∇2T0−4κσT4 m−4κσT3 mT0 (2.24) Se introducen las variables adimensionales siguientes: λ=16σκT3 mH2 K;Rc =4HσT3 m K (2.25) x=x∗H;z=z∗H;t=H2 α;vx=α Hv∗ x;vz=α Hv∗ z; ∆ ¯ T=T1−T2 Tm θe=T0 e Tm ;T0=∆¯ TTmT∗ Ra ;T0 e=∆¯ TTmT∗ e Ra ; (2.26) y la llamada función de corriente Ψ∗ , que se dene: v∗ x=∂Ψ∗ ∂z∗;v∗ z=−∂Ψ∗ ∂x∗, (2.27) que automáticamente satisface la ecuación 2.15 en el caso 2D en el plano xz . Introduciendo las variables adimensionales anteriores, tomando el caso bidimensional en el plano xz , derivado la ecuación 2.16 con respecto a z∗ y la ecuación 2.18 con respecto a x∗ y restándolas, se obtiene un sistema de dos ecuaciones con dos incógnitas: ∂∇2Ψ ∂t +∂Ψ ∂z ∇2∂Ψ ∂x −∂Ψ ∂x ∇2∂Ψ ∂z = Pr∇4Ψ−P r∂T ∂x −PrRa ∆¯ T ∂θe ∂x , (2.28) ∂T ∂t +∂Ψ ∂z ∂T ∂x −∂Ψ ∂x ∂T ∂z = −∂Ψ ∂z ∂θe ∂x −Ψ ∂x ∂θe ∂z Ra ∆¯ T+∇2T−λT, (2.29) donde se han omitido los superíndices ∗ por comodidad. Por otro lado, la ecuación que satisface la ecuación de equilibrio 2.23se reduce a: ∂2θe ∂x∗2+∂2θe ∂z∗2−λθe=λ 4, (2.30) 18 2.2. Ecuaciones del problema de convección con ujo de radiación. Aproximaciones 19 cuyas condiciones de contorno a satisfacer, en variables adimensionales se expresan: ∂θe ∂x∗ x∗=0 x∗=A = 0, (2.31) θe|x∗=0 x∗=A=±∆¯ T 2 (2.32) habiéndose asumido el caso en que no inuyen las paredes (Rc nulo), pues la inuencia de esta se considera a partir del capítulo 3. La ecuación 2.30 es una ecuación diferencial que se resolverá por el método de colocación , explicado en el capítulo 3 Método numérico , cuyos resultados para distintos valores del parámetros λ se representan en las guras de 2.3. Para ello, será necesario aplicar las condiciones de contorno de temperatura de equilibrio 2.31 y 2.32. (a) (a) λ= 0 , Rc=0 (b) (b) λ= 3 , Rc=0 (c) (c) λ= 10 , Rc=0 Figura 2.3: Solución de la temperatura de equilibrio θe con parámetro de radiación-conducción Rc nulo para 3 valores del parámetro λ . 19 20 Capítulo 2. Formulación del problema Por su parte, el problema de movimiento 2.28-2.29 queda determinado si se añaden las condiciones de contorno apropiadas: Para t∗= 0; 0 ≤x∗≤H L,0≤z∗≤1→Ψ∗= Ψ∗ 0(x∗, z∗), T∗= 0, (2.33) Para x∗= 0 y x∗= L/H=A→Ψ∗= 0,∂Ψ∗ ∂x∗= 0,∂T∗ ∂x∗ x∗=0 x∗=A =±RcT∗, (2.34) Para z∗= 0 y z∗= 1 →Ψ∗= 0,∂Ψ∗ ∂z∗= 0, T∗= 0, (2.35) donde las condiciones de contorno 2.34 y 2.35 se corresponden al caso en el que las paredes inferior y superior son rígidas , es decir, la componente horizontal del campo de velocidades es nula vx=−∂Ψ∗/∂z = 0 . En el apartado de Condiciones de contorno dentro del capítulo de Método numérico se verán el caso de ambas paredes libres y el de pared inferior rígida y superior libre. Las ecuaciones 2.28 y 2.29 se resolverán numéricamente en los siguientes capítulos, teniendo en cuentas las condiciones de contorno e inicial que se acaban de exponer. 2.2.4. Aproximación ópticamente gruesa u opaca El segundo caso abordado es aquel en el que el camino libre medio κ−1 es de un orden muy inferior al orden de magnitud de la altura de la cavidad. Dicha condición puede ser expresada matemáticamente como κ∗h >> a , donde h representa la altura de la cavidad en la que se mueve el uido y a es un número característico al que ya se ha hecho referencia en el apartado anterior. Recuérdese que dicho número es de orden unidad [16]. Siguiendo la misma estructura que con anterioridad, en primer lugar se tiene, a partir de una ecuación 2.9, para valores altos de κ , se puede desarrollar RH en series de potencias en términos de κ−2 una solución formal de dicha ecuación viene dada por Goody (ec. 13 del documento) [8]: I(~s) = e−κs Zs q κeκσB(σ)dσ, (2.36) donde la contribución del límite inferior se debe determinar a partir de las condiciones de contorno, pero dichas condiciones de contorno solo contribuyen de manera apreciable a distancias menores que κ−1 del contorno, por lo que se puede despreciar su contribución fuera de dichas regiones. Si se integra de forma reiterada la ecuación anterior, se obtiene que: I(~s) = B−κ−1dB ds +κ−2d2B ds2−κ−3d3B ds3... (2.37) Entonces, usando la ecuación 2.8 y sabiendo que dB ds =µ1 ∂B ∂x +µ2 ∂B ∂y +µ3 ∂B ∂z , (2.38) donde µ1, µ2y µ3 son los cosenos direccionales de ~s . Al realizar las integrales de 2.37 y utilizar la ecuación de Stefan-Boltzmann 2.10, se obtiene que el término de transferencia de energía por radiación RH para la aproximación opaca es: H(s) = κ−14σ 3∇2T4. (2.39) 20 2.2. Ecuaciones del problema de convección con ujo de radiación. Aproximaciones 21 Para realizar desarrollar la ecuación 2.39 es necesario desarrollar cómo va a ser el campo de temperaturas de manera cualitativa mediante el uso de una serie de hipótesis sobre dicho campo , análogas a las del apartado anterior: T(~x, t) = Te(~x, t) + T0(z); donde T0<< T,Te (2.40) Te(z) = Tm+T0 e(z); donde T0 e<< Tm,Te (2.41) T1−T2 T1 << 1; T1−T2 T2 << 1; (2.42) De esta manera, introduciendo 2.40 y 2.41 en la ecuación de la energía 2.2 (2D) se obtiene: ∂T ∂t +vx ∂T ∂x +vz ∂T ∂z =α∇2T+16σ 3κ∇ · (Te3∇Te +Te3∇T0+ 3Te2T0∇T e), (2.43) donde se ha hecho T3∼ =Te3+ 3Te2T0 y se han despreciado los productos entre variaciones. Para hallar Te(z) , se impone el estado estático en la ecuación anterior: ~v =~ 0; T0= 0; (2.44) se determina que se cumple: d dz αdTe dz +16σ 3κTe3dTe dz = 0, (2.45) y en virtud de la hipótesis 2.41, si se sustituye Te3 por Tm3 , se tiene que: dTe dz =C. (2.46) De esta forma, aplicando las condiciones de contorno T(0) = T1 y T(H) = T2 , se obtiene denitivamente Te(z) : Te(z) = T2−T1 Hz+T1. (2.47) Entonces, la ecuación de la energía se convierte en: ∂T0 ∂t +vx ∂T0 ∂x +vz ∂T0 ∂z =   α∇2Te +   16σ 3κ∇2Te4 4+∇2T0 +16σ 3κ(2∇Te3· ∇T0+Te3∇2T0+ 3T0(2Te∇Te · ∇T e +  Te2∇2Te), (2.48) donde se han utilizado las ecuaciones 2.45 y 2.47 para la cancelación de términos. ∂T0 ∂t +~v · ∇T0+vz T2−T1 H= α+16σ 3κρ0cp Tm3∇2T0+32σ 3κρ0cp Tm2T2−T1 H ∂T0 ∂z +T0T2−T1 H232σ 3κρ0cp Tm (2.49) 21 28 Capítulo 3. Método numérico En la siguiente sección se detalla cómo se va a utilizar el método explicado para resolver las ecuaciones Saltzmann para la convección y el sistema equivalente obtenido en el capítulo 2 para los casos de radiación con aproximación transparente y opaco. 3.2. Resolución numérica En esta sección se va a explicar de qué manera se han discretizado las ecuaciones 2.28 y 2.29 para aproximar las derivadas espaciales y las derivadas con respecto al tiempo, que se aproximarán mediante diferencias regresivas . Dichas ecuaciones son las que se obtuvieron en el apartado 2.2.3. Posteriormente, se verá que solo hay que introducir ligeros cambios en los programas para tener en cuenta distintas condiciones de contorno y según si se quiere estudiar el caso de temperatura impuesta en las paredes horizontales o en las paredes verticales y paredes rígidas o libres. Para discretizar las ecuaciones anteriores a la hora de desarrollar códigos numéricos y obtener soluciones aproximadas, se van a introducir las siguientes aproximaciones para las derivadas espaciales: ∂ ∂x∗F∼ =Dx∗F;∂ ∂z∗F∼ =Dz∗F;∇2F∼ =(D2 x+D2 z)∗F=DL∗F;∇4F∼ =D2 L∗F, (3.9) donde F representa a Ψ o a T∗ mientras que para las derivadas con respecto al tiempo se utiliza el Método de Diferencias regresivas ([18], capítulo 1), por lo que los términos en los que aparece ∂F ∂t∗ , se tiene: ∂F ∂t∗=F|tn −Ftn−1 ht ; dondeht<< 1. (3.10) Por último, téngase en cuenta que en las ecuaciones aparecen productos entre las derivadas de las funciones Ψ∗ y T∗ , siendo estos productos términos de carácter no lineal si se sustituyen de forma directa por Ψ∗|tn y T∗|tn , respectivamente, por lo que en dichos términos se realizarán las aproximaciones: Ψ∗|tn ≈Ψ∗|tn−1;T∗|tn ≈T∗|tn−1, (3.11) que serán válidas siempre que ht sea lo sucientemente pequeño . Teniendo en mente todo lo expuesto con anterioridad, se está en condiciones de escribir qué ecuaciones se obtienen nalmente para poder ser implementadas en un lenguaje de programación que permita resolver sistemas de ecuaciones lineales de manera práctica (MATLAB en este caso). Así, se tiene: DL∗Psi −Psinm1 ht +Dz∗Psinm1∗Dx∗DL∗Psinm1−Dx∗Psi ∗Dz∗DL∗Psi | {z } NLP SI = PrD2 LPsi −PrDxTp −Pr ∗Ra/DeltaTb ∗DxThetae; (3.12) Tp −Tpnm1 ht +Dz∗Psinm1∗Dx∗Tpnm1−Dx∗Psinm1∗Dz∗Tpnm1 | {z } NLT = −(Dz∗Psi ∗Dx∗Thetae −Dx∗Psi ∗Dz∗Thetae)∗Ra/DeltaT b +DL∗Tp −lambda ∗Tp, (3.13) donde se han utilizado los nombres de las variables que se les ha dado en el código y donde los términos NLPSI y NLT reciben esos nombres pues corresponden con términos que serían no lineales si no se hubiera realizado en ellos las aproximaciones Ψ∗|tn ≈Ψ∗|tn−1yT∗|tn ≈T∗|tn−1 . De esta manera, se tiene un sistema de ecuaciones, donde los vectores: 28 3.2. Resolución numérica 29 Psi =  Psi1 ... PsiNt   Tp =  Tp1 ... TpNt   son las 2×Nt incógnitas de un sistema de 2×Nt ecuaciones, cuya forma compacta escrita se obtiene a partir de las ecuaciones 3.12 y 3.13 y se verá implementada en la matriz del sistema en el código que se adjunta más adelante. Convección: programa para la convección. Si se quieren obtener resultados para el caso de convección pura sin tener en cuenta efectos de radiación, tan solo hay que hacer en el programa las variables lambda y Rc nulas. En los resultados que se expongan en los próximos apartados, se reproducen los obtenidos en [1] y se comparan con resultados que se obtengan con valores no nulos de las variables anteriores. 3.2.1. Condiciones de contorno e inicial. Problema con paredes no participativas En el apartado anterior se ha explicado cómo se discretizan las variables de las ecuaciones que aparecen en el problema de radiación-convección. En este, se expone cómo se han impuesto numéricamente las condiciones de contorno e inicial en el problema de radiación en el que el medio participa, λ > 0 , pero las paredes no absorben radiación, Rc = 0 . El cumplimiento de dichas condiciones no es trivial. Esto es debido a que se tienen condiciones de contorno sobre las derivadas de Ψ∗ y T∗ a la vez que sobre las propias funciones, por ejemplo, para el caso de bordes rígidos/rígidos. Al no ser posible imponer dos valores distintos a la vez sobre los vectores solución de Ψ∗ y T∗ , la forma de resolver este conicto es imponer las condiciones de contorno de las funciones sobre los puntos que se encuentran sobre los contornos, mientras que para las condiciones sobre las derivadas se usan los puntos de los subcontornos. En denitiva, los pasos a seguir serían: Antes que nada es esencial aclarar cuáles son los valores de los índices I para cada pared del contorno. A saber: I=       (i−1) ∗Nz + 1 si z∗= 0 (i−1) ∗Nz +Nz =i∗Nz si z∗= 1 jsi x∗= 0 (Nx −1) ∗Nz +jsi x∗=L/H donde i = 1, ..., Nx y j = 1, ..., Nz (3.14) Ahora, en primer lugar, debido a que Matlab trabaja de forma mucho más eciente "llenandoçolumnas que las, se transponen todos los vectores y matrices que forman parte del sistema general. En segundo lugar, para imponer las condiciones de Ψ∗ , en el caso de que se esté estudiando el caso de paredes horizontales rígidas se procede de la siguiente forma: se anula en las matrices (transpuestas) APSI_tr, ATp_tr, BPSI_tr y BTp_tr las columnas cuyo índice I corresponde con el de los puntos de los contornos horizontales. Posteriormente, se impone un valor unidad en los elementos APSI_(I,I) , así se impone que Ψ∗ I= 0 una vez que se transponga de nuevo la matriz original del sistema. Análogamente se procede para imponer que T∗= 0 en los contornos horizontales: se impone que los elementos BT_tr(I,I)=1. Además, para asegurar que para todo instante t se cumplen las condiciones de contorno, se crean dos factores de condiciones de contorno o boundary condition factors , BCF_PSI_tr y BCF_Tp_tr cuyo valor es nulo para todo índice I perteneciente a cualquier punto de los contornos horizontales. Ahora hay que imponer que en las paredes verticales la derivada de T∗ es nula y que Ψ∗ es también nula. Entonces, tal y como se hizo en el punto anterior, se impone en los puntos de los contornos verticales que las columnas de las matrices (transpuestas) APSI_tr, ATp_tr, BPSI_tr y BT_tr cuyo índice sea I , sean nulas. Posteriormente, se hacen igual a la unidad los elementos APSI_tr(I,I) para imponer Ψ∗ I= 0 29 30 Capítulo 3. Método numérico tras la transposición de la matriz del sistema. Para imponer que ∂T∗ ∂x∗= 0 se hace BT_tr(:,I)=Dx_tr(:,I). Se hacen nulos también los términos BCF_PSI_tr(I) y BCF_Tp_tr(I). Por último, quedaría imponer la condición de que las paredes horizontales y verticales son rígidas, por lo que habría que imponer en los puntos del contorno que las velocidades v∗ x y vz son nulas, respectivamente. Para ello, teniendo en cuenta que v∗ x=−∂Ψ∗/∂z y v∗ z=∂Ψ∗/∂x , en principio habría que proceder como en los dos puntos anteriores para imponer condiciones de contorno. Sin embargo, no queda más remedio que hacer uso de los puntos del subcontorno, pues los del contorno ya están utilizados. Para ello, se anulan los las columnas de todas las matrices que componen la matriz completa del sistema con índices I pertenecientes a dichos subcontornos, donde es evidente que las esquinas son comunes a los horizontales y verticales y en este caso se han incluido en los primeros. En el caso de subcontorno inferior se tiene I= (i−1) ∗Nz + 2 y en el superior I= (i−1) ∗Nz +Nz −1 con i=2:Nx −1 , mientras que para los subcontornos izquierdo y derecho, se tiene, respectivamente, I=Nz +j y I= (Nx −2) ∗Nz +j , donde j= 3 : Nz −2 . Después, en la matriz APSI_tr, se impone en las columnas del subcontorno vertical es igual a la columna con el mismo índice de la matriz Dz_tr y en las horizontales se hacen igual a las columnas de Dx_tr. Con todo esto ya se tienen impuestas las condiciones de contorno. En el siguiente punto se discute qué se ha de hacer si se tiene alguna pared libre (sin esfuerzos tangenciales) en vez de paredes rígidas. En el caso de tener una pared libre, por ejemplo la horizontal inferior, los esfuerzos tangenciales han de ser nulos en la pared: τ=µ∂vx ∂z z=0 = 0 ⇒∂2Ψ∗ ∂z∗z∗=0 = 0 (3.15) Por tanto, si los esfuerzos son nulos, habrá que imponer, en vez de lo dicho en el punto anterior, en el que se habló del caso rígido-rígido en los contornos horizontales, que la columna con índice I perteneciente al subcontorno horizontal de la matriz APSI_tr es igual a la columna del mismo índice de la matriz que proporciona la segunda derivada con respecto de z , Dz2_tr. Análogamente se haría con la pared horizontal superior o, en el caso de paredes libres verticales, en cuyo caso se necesitaría igualar las columnas correspondientes a las de la matriz que proporciona la segunda derivada con respecto x , Dx2_tr. Para las condiciones iniciales, se ha usado una función de corriente que cumpla con las condiciones de contorno y lo mismo se ha hecho con el campo de temperaturas (ver [1], capítulo 4): Ψ∗(x, z, t = 0) = 0.0005x∗2(x∗ max −x∗)2z∗2(z∗ max −z∗)2cos πx∗z∗ 5x∗ maxz∗ max  (3.16) T∗(x, z, t = 0) = 0 (3.17) A continuación se muestran los códigos principales que se han utilizado para resolver el problema de convección sin radiación con paredes horizontales rígidas. A partir de él, con las modicaciones explicadas en las condiciones de contorno y cambiando las ecuaciones según se expuso en la sección de resolución numérica en función de la aproximación de radiación que se quiera estudiar, es suciente para obtener nuevos códigos que permitan efectivamente resolver los problemas. Programa básico Se utiliza el siguiente programa como base y ciertas modicaciones para resolver los distintos apartados. 1 2 % _____PARAMETROS_____ % 3 Pr=0.73; 4 5 %Ra=1565+1000−500−300+200+50+500−1000; 6 DeltaTb=0.1; Rc=0; 7 Ra=8800+15; 30 3.2. Resolución numérica 31 8 lambda_G=5; lambda=lambda_G^2; 9 10 dt=0.0001; 11 % ______GEOMETRIA_____ % 12 13 zmin=0; zmax=1; 14 15 Nz=30 16 17 xmin=0; xmax=.5; 18 19 Nx=30 20 21 tic 22 23 [Lpx,Dx,Nx,xch]=matricesx(Nx,Nz,xmin,xmax); 24 [Lpz,Dz,Nz,zch]=matricesz(Nz,Nx,zmin,zmax); 25 % _____MATRICES____ % 26 Nt=Nz*Nx; 27 28 DL=Dx*Dx+Dz*Dz; 29 DL2=DL*DL; 30 DxDL=Dx*DL; 31 DzDL=Dz*DL; 32 Dz2=Dz*Dz; 33 34 %%SPARSE 35 DL=sparse(DL); 36 DL2=sparse(DL2); 37 DxDL=sparse(DxDL); 38 DzDL=sparse(DzDL); 39 Dz2=sparse(Dz2); 40 41 %%TRANSPUESTAS 42 DL_tr=DL'; DL2_tr=DL2'; DxDL_tr=DxDL'; 43 DzDL_tr=DzDL'; Dz2_tr=Dz2'; Dx_tr=Dx'; 44 Dz_tr=Dz'; 45 46 %% % %SOLUCION EQUILIBRIO % % % % 47 Ae=DL−lambda*eye(Nt); 48 be=lambda*ones(Nt,1)/4; 49 50 % Modificaciones de Ae y be debidas a condiciones de contorno: 51 % Contornos horizontales (insulated) 52 53 Ae_tr=Ae'; be_tr=be'; 54 for i=2:Nx−1, 55 I1=(i−1)*Nz+1; 56 I2=(i−1)*Nz+Nz; 57 % Radiacion conveccion horiz. 58 % Ae(I1,:)=Dz_tr(I1,:); Ae(I1,I1)=Ae(I1,I1)−Rc; be(I1,1)=0; 59 % Ae(I2,:)=Dz_tr(I2,:); Ae(I2,I2)=Ae(I2,I2)+Rc; be(I2,1)=0; 60 % Temperatura impuesta horiz. 61 Ae_tr(:,I1)=0; Ae_tr(I1,I1)=1; be_tr(1,I1)=DeltaTb/2; 62 Ae_tr(:,I2)=0; Ae_tr(I2,I2)=1; be_tr(1,I2)=−DeltaTb/2; 63 end 64 % Contornos verticales (imposed temperature): 65 % 66 for j=1:Nz, 67 I1=j; 68 I2=(Nx−1)*Nz+j; 31 32 Capítulo 3. Método numérico 69 Ae_tr(:,I1)=Dx_tr(:,I1); Ae_tr(I1,I1)=Ae_tr(I1,I1)−Rc; be_tr(1,I1)=0; 70 Ae_tr(:,I2)=Dx_tr(:,I2); Ae_tr(I2,I2)=Ae_tr(I2,I2)+Rc; be_tr(1,I2)=0; 71 end 72 % 73 Ae=Ae_tr'; be=be_tr'; 74 Thetae=Ae\be; 75 76 DxThetae=Dx*Thetae; DzThetae=Dz*Thetae; 77 for j=1:Nz, 78 Thetaemat(1:Nx,j)=Thetae(((1:Nx)−1)*Nz+j,1); 79 end 80 for I=1:Nt, 81 aux1(I,:)=Dz(I,:)*DxThetae(I,1); 82 aux2(I,:)=Dx(I,:)*DzThetae(I,1); 83 end 84 for j=1:Nz, 85 xmat(1:Nx,j)=xch(1:Nx)'; zmat(1:Nx,j)=zch(j); 86 end 87 mesh(xmat,zmat,Thetaemat) 88 89 % _____CONDICIONES DE CONTORNO_____ % 90 %Definimos las matrices necesarias: 91 APSI=DL−dt*Pr*DL2; 92 AT=dt*Pr*Dx; 93 BPSI=Ra/DeltaTb*dt*(aux1−aux2); 94 BT=eye(Nt)*(1+dt*lambda)−dt*DL; 95 %Factor de condiciones de contorno 96 FBC_PSI=ones(Nt,1); 97 FBC_T=ones(Nt,1); 98 99 APSI_tr=APSI'; AT_tr=AT'; BPSI_tr=BPSI'; BT_tr=BT'; 100 FBC_PSI_tr=FBC_PSI'; FBC_T_tr=FBC_T'; 101 102 %Ponemos las condiciones que hacen Psi y T igual a cero. Esto se da en 103 %contornos horizontales: 104 for i=2:(Nx−1) 105 %Contornos de abajo 106 I=(i−1)*Nz+1; 107 %Imponemos Psi y T: 108 %psi=0 109 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI_tr(1,I)=0; AT_tr(:,I)=0; 110 %T=0 111 BPSI_tr(:,I)=0; BT_tr(:,I)=0; BT_tr(I,I)=1; FBC_T_tr(1,I)=0; 112 % % ConRad 113 % BPSI_tr(:,I)=0; BT_tr(:,I)=Dz_tr(:,I); BT_tr(I,I)=BT_tr(I,I)−Rc; FBC_T_tr(1,I) =0; 114 %Contorno superior 115 I=(i−1)*Nz+Nz; 116 %psi=0 117 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI_tr(1,I)=0; AT_tr(:,I)=0; 118 %T=0 119 BPSI_tr(:,I)=0; BT_tr(:,I)=0; BT_tr(I,I)=1; FBC_T_tr(1,I)=0; 120 % % ConvRad 121 % BPSI_tr(:,I)=0; BT_tr(:,I)=Dz_tr(:,I); BT_tr(I,I)=BT_tr(I,I)+Rc; FBC_T_tr(1,I) =0; 122 end 123 % 124 %Ahora vamos a contornos verticales 125 for j=1:Nz 126 %Contornos de izq 127 I=j; 32 3.2. Resolución numérica 33 128 %Imponemos Psi y T: 129 %psi=0 130 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI_tr(1,I)=0; AT_tr(:,I)=0; 131 %dT/dx=0 132 BPSI_tr(:,I)=0; BT_tr(:,I)=Dx_tr(:,I); BT_tr(I,I)=BT_tr(I,I)−Rc; FBC_T_tr(1,I)=0; 133 %T=0 134 % BPSI_tr(:,I)=0; BT_tr(:,I)=0; BT_tr(I,I)=1; FBC_T_tr(1,I)=0; 135 %Contornos de dcha 136 I=(Nx−1)*Nz+j; 137 %Imponemos Psi y T: 138 %psi=0 139 APSI_tr(:,I)=0; APSI_tr(I,I)=1; FBC_PSI_tr(1,I)=0; AT_tr(:,I)=0; 140 %dT/dx=0 141 BPSI_tr(:,I)=0; BT_tr(:,I)=Dx_tr(:,I); BT_tr(I,I)=BT_tr(I,I)+Rc; FBC_T_tr(1,I)=0; 142 % %T=0 143 % BPSI_tr(:,I)=0; BT_tr(:,I)=0; BT_tr(I,I)=1; FBC_T_tr(1,I)=0; 144 end 145 %Falta por aplicar velocidades nulas en las paredes. Para ello, se usan las 146 %filas correspondientes a los subcontornos de abajo. 147 %Contornos horizontales d^2psi/dz^2=0 148 for i=3:(Nx−2) 149 K=(i−1)*Nz+2; 150 APSI_tr(:,K)=Dz_tr(:,(K−1)); FBC_PSI_tr(1,K)=0; AT_tr(:,K)=0; 151 K=(i−1)*Nz+Nz−1; 152 APSI_tr(:,K)=Dz_tr(:,(K+1)); FBC_PSI_tr(1,K)=0; AT_tr(:,K)=0; 153 end 154 %Contornos verticales dpsi/dx=0 155 for j=2:(Nz−1) 156 K=(2−1)*Nz+j; 157 APSI_tr(:,K)=Dx_tr(:,(K−Nz)); FBC_PSI_tr(1,K)=0; AT_tr(:,K)=0; 158 K=(Nx−2)*Nz+j; 159 APSI_tr(:,K)=Dx_tr(:,(K+Nz)); FBC_PSI_tr(1,K)=0; AT_tr(:,K)=0; 160 end 161 162 APSI=APSI_tr'; AT=AT_tr'; BPSI=BPSI_tr'; BT=BT_tr'; 163 FBC_PSI=FBC_PSI_tr'; FBC_T=FBC_T_tr'; 164 165 %Formamos la matriz del sistema: 166 Asyst=[APSI AT; BPSI BT]; 167 Asystm1=Asyst\eye(size(Asyst)); 168 toc 169 for i=1:Nx, 170 for j=1:Nz, 171 I=(i−1)*Nz+j; 172 Psi(I,1)=0.05*xch(i)^2*(xmax−xch(i))^2*zch(j)^2*(zmax−zch(j))^2; 173 T(I,1)=0.05*sin(3*pi*zch(j)); 174 end 175 end 176 AbPsi=sparse([DL, zeros(Nt,Nt)]); 177 AbT=sparse([zeros(Nt,Nt), speye(Nt)]); 178 [X, Z]=meshgrid(xch,zch); 179 Ntime=5e5; 180 for nt=1:Ntime, 181 nt; 182 t=nt*dt; 183 %tic 184 bNLPsi(1:Nt,1)=(−(Dz*Psi).*(DxDL*Psi)+(Dx*Psi).*(DzDL*Psi))*dt−Pr*Ra/DeltaTb*dt* DxThetae; 185 %tiempo1=toc 186 187 bNLT(1:Nt,1)=(−(Dz*Psi).*(Dx*T)+(Dx*Psi).*(Dz*T))*dt; 33 34 Capítulo 3. Método numérico 188 189 bLinPsi(1:Nt,1)=AbPsi*[Psi ; T]; 190 bLinT(1:Nt,1)=AbT*[Psi ; T]; 191 bsyst(1:(2*Nt),1)=[FBC_PSI.*(bNLPsi+bLinPsi); FBC_T.*(bNLT+bLinT)]; 192 bn=Asystm1*bsyst; Psi=bn(1:Nt,1); T=bn((Nt+1):(2*Nt),1); 193 194 if nt==10 || nt==Ntime*0.2 || nt==Ntime*0.4 || nt==Ntime/2 || nt==Ntime*0.6 || nt== Ntime*0.8 || nt==Ntime 195 nt 196 for j=1:Nz, 197 Psimat(1:Nx,j)=Psi(((1:Nx)−1)*Nz+j,1); Tmat(1:Nx,j)=T(((1:Nx)−1)*Nz+j,1); 198 Tphysmat(1:Nx,j)=−1/2+Thetaemat(1:Nx,j)/DeltaTb+Tmat(1:Nx,j)/Ra; 199 end 200 subplot (2,2,1) 201 contour (X, Z, Psimat') 202 subplot (2,2,2) 203 contour (X, Z, Tphysmat') 204 Tmat_c(nt)=Tmat(ceil(Nx/2),ceil(Nz/2)); 205 Tmat_max=max(abs(Tmat_c(1:nt))); 206 subplot (2,2,3) 207 plot((1:nt)*dt,Tmat_c(1:nt)/Tmat_max) 208 axis([0 Nt*dt −2 2]) 209 subplot (2,2,4) 210 plot (Tphysmat(ceil(Nx/2),:),zch,'r') 211 axis([−2 2 0 1]) 212 pause(0.00001) 213 hold off 214 end 215 end A continuación se exponen también los códigos de los que hace uso el programa principal para obtener las matrices que permiten calcular las derivadas de forma numérica, matricesx (línea 6 y 7 del código). Es análogo para la variable z . 1 %% % %SE OMITEN TILDES % % % % 2 %cambiar las condiciones como se desee 3 function [Lpx,Dx,Nx,x]=matricesx(Nx,Nz,x1,xn) 4 x=zeros(1,Nx); x(1)=x1; x(Nx)=xn; 5 %nodos de chebyshev 6 for cl=1:Nx 7 x(cl)=(x(Nx)+x(1))/2−(x(Nx)−x(1))/2*cos((cl−1)/(Nx−1)*pi); 8 end 9 % 10 %derivada interpolantes 11 num=0; 12 for j=1:Nx 13 for i=1:Nx 14 den=x(i)−x; den(i)=1; den=prod(den); flag=1; 15 if i==j %la expresion es distinta cuando coinciden los dos indices, pues aparecen mas terminos en el numerador 16 sumnum=0; 17 for k1=1:Nx 18 if k1~=i 19 num=x(j)−x; num(j)=1; num(k1)=1 ; sumnum=prod(num)+sumnum; 20 end 21 end 22 num=sumnum; 23 Lpx(j,i)=num/den; 24 else 25 num=(x(j)−x); 26 num(j)=1; num(i)=1; 34 3.2. Resolución numérica 35 27 num=prod(num); 28 Lpx(j,i)=num/den; 29 30 end 31 Lpx(j,i)=num/den; %matriz con las derivadas de las funciones de Lagrange 32 end 33 end 34 Lpx; 35 Dx=kron(Lpx,eye(Nz)); 36 end Donde en la línea 35 aparece la función de Matlab kron(A,B) , que calcula el producto de Kronecker 2 de dos matrices [19], lo cual resulta muy eciente para calcular una matriz de tamaño (Nx×Nz)×(Nx×Nz) = Nt×Nt , en el caso de Dx , y, análogamente, la matriz Dz , de tamaño (Nz×Nx)×(Nz×Nx) = Nt×Nt . 3.2.2. Condiciones de contorno e inicial en el problema con paredes participativas en la radiación En esta sección se trata el problema en el que las paredes son capaces de emitir y absorber radiación ( Rc >0 ). Es destacable sin duda la manera en la que se ha modelado dicho intercambio de calor entre elementos de pared , habiéndose apoyado en la referencia [21] (veánse los capítulos 12 y 13) y [22] (especialmente sección 5.3). Las claves aquí proporcionadas relevantes para este estudio son: En primer lugar, dividir la supercie de las paredes en elementos de pared, que se han elegido de tal modo que están formados por el tramo comprendido entre el punto medio entre dos nodos consecutivos, excepto en el caso de ser el primer o último trozo de pared dentro de uno de los contornos verticales (u horizontales), es decir, si un trozo de pared está comprendido entre el nodo de una esquina y el anterior (o posterior en su caso), en cuyo caso el trozo de pared abarca desde el nodo en la esquina hasta el punto medio entre los dos nodos anteriores (o posteriores en su caso). Para una mejor explicación y entendimiento, ver gura 3.2. Además, se han considerado supercies difusas , entonces sus propiedades son independientes de la dirección, y grises , siendo entonces también independientes de las longitudes de onda, además de ser uniformes en cada trozo de supercie las radiaciones entrantes y salientes. Asimismo, como consecuencia directa de las dos hipótesis anteriores, en un material así descrito se tiene que su absortividad αi es igual a su emisividad i . 2 El producto de Kronecker de dos matrices A ( p×q )yB( m×n ) es una matriz de tamaño (p×m)×(q×n) : A⊗B=          a11   b11 ... b1n ... ... ... bm1... bmn  ... a1q  b11 ... b1n ... ... ... bm1... bmn   ... ... ... ap1  b11 ... b1n ... ... ... bm1... bmn  ... apq   b11 ... b1n ... ... ... bm1... bmn            35 36 Capítulo 3. Método numérico Figura 3.2: Representación esquemática explicativa sobre la división de los contornos horizontales y verticales en elementos de pared. Los índices de elementos en la pared son los representados dentro de un círculo, correspondiendo los otros números al índice de nodos de contornos, no incluyendo las esquinas. Obsérvese que el sentido seguido es horario y que habrá 2Nz+ 2Nx−8 elementos de pared, por Nw= 2Nz+ 2Nx−4 nodos en los contornos. En segundo lugar, para tener en cuenta la orientación de las distintas divisiones que conforman las paredes, se tiene un parámetro que solo depende de la geometría de la cavidad, conocido como factor de de visión (ver gura 3.3) o de forma ( view or shape factor ): Fij : fracción de radiación que sale de la supercie i y llega directamente a la supercie j . Figura 3.3: Factor de forma entre supercies. Estos factores de forma cuentan con dos propiedades que se han usado, demostradas en la referencia citada (Cengel, 2006) [21]: FijAi=FjiAj, (3.18) conocida como relación de reprocidad y Ne X j Fij = 1, (3.19) donde Ne es el número total de elementos de pared en la cavidad. 36 3.2. Resolución numérica 37 Además, se aporta la regla necesaria para el cálculo de los factores de forma para cavidades bidimensionales, que fue dada por Hottel , denida como: Fij =PLc−PLnc 2×Li , (3.20) donde Lc representa la longitud de cuerdas cruzadas, Lnc , la longitud de cuerdas no cruzadas y Li , la longitud de cuerda de la supercie i , que se pueden ver en la gura 3.4. Figura 3.4: Determinación del factor de forma entre dos supercies: longitudes de cuerdas cruzadas ( L5y L6 ), no cruzadas ( L3y L4 ) y de las supercies ( L1y L2 ). En tercer lugar, teniendo en cuenta que la cantidad de radiación que sale de un cuerpo es la suma de su calor reejado y del emitido, teniendo en cuenta las hipótesis del primer punto, la ecuación de Stefan-Boltzmann 2.10 y denición de radiosidad J , se tiene que la radiosidad es igual a: Ji=iσT4+ (1 −i)Gi, (3.21) donde el primer término representa la potencia radiativa de la supercie i y Gi=PN j=1 JjFij es la irradiancia , es decir, la potencia incidente en la supercie i por unidad de área. El balance de energía en un elemento de pared puede obtenerse por un lado a partir de las radiosidades del resto, su propia radiosidad y la ganancia neta Qij de calor por otros medios distintos a la radiación, como por ejemplo la convección. Este último, en el caso de una pared aislada, será igual a la pérdida neta de calor debida a la radiación que emite, dada por la ley de Stefan-Boltzman. El balance queda por tanto: Qi Ai −Ji+ N X j=1 JjFij = 0 ⇒Qi Ai = N X j=1 (Ji−Jj)Fij.. (3.22) A su vez, el balance de energía neta en un elemento de pared también se puede calcular a partir de la energía que sale de ella y la energía que absorbe: Qi Ai =iσT4 i−iGiAi. (3.23) Tras exponer los puntos anteriores, se trata ahora de encontrar una ecuación que sirva como condición de contorno a cumplir por parte los elementos de pared, imponiéndola en los nodos. En los nodos pertenecientes a las paredes aisladas se conoce el calor que intercambian con el uido por convección Qi/Ai y se les denotará por ic (o jc ), mientras que los nodos que pertenezcan a paredes con temperatura impuesta, se denotarán iT (o jT ). Consecuentemente, teniendo en cuenta que en los nodos cuya temperatura es conocida, si se sustituye la irradiancia en función de la radiosidad, se pueden calcular estas últimas mediante: Ji =Arad−1 ij bi, (3.24) donde Aradij =δij −Fij(1 −i) , y bi=σiT4 i . Así se pueden determinar los valores del calor por unidad de área en dichos nodos: Qi Ai =AQijbi, (3.25) 37 44 Capítulo 3. Método numérico 44 Capítulo 4 Medio participativo dentro de una cavidad con paredes con absorción de radiación nula 4.1. Paredes verticales adiabáticas En esta sección se exponen resultados sobre el problema de convección con aproximación de medio transparente cuando las paredes que están aisladas no tienen capacidad de absorber o emitir radiación, es decir, cuando la emisividad de dichas paredes  es nula (materiales como el oro pulido, el cobre pulido y el aluminio tienen emisividades casi nulas), condición que según el modelo desarrollado se corresponde con un valor parámetro de radiación-convección Rc nulo , para cavidades con temperatura impuesta en las paredes horizontales y verticales adiabáticas y viceversa. En particular, se estudia cómo afecta aumentar la importancia del parámetro de radiación, λ↑ , al número de Rayleigh crítico para diferentes valores de la relación de aspectos A=L/H , donde xmax =L y zmax =H para una cavidad rectangular con paredes libres o rígidas, siendo x la coordenada horizontal y z la vertical (véase gura 2.4). Además, se estudia también cómo dicho parámetro inuyen en la λc , el número de Nusselt, etc. Este problema, como se expuso en el capítulo 2, está gobernado por las ecuaciones 2.28 y 2.29, bajo las condiciones de contorno 2.31 y 2.32. 4.1.1. Número de Rayleigh crítico: resultados fundamentales Este apartado está dedicado al caso de paredes verticales adiabáticas, pues no existen una valor crítico del número de Rayleigh si las paredes con temperatura impuesta son las verticales, pues ya se ha comentado que siempre es inestable dicho problema. Se tiene como objetivo estudiar cuál es el número de Rayleigh a partir del cual las perturbaciones iniciales se propagan de forma inestable, para el caso de cavidad con paredes horizontales libres en cavidades de gran relación de aspectos A= 10 , ver 4.1, y compararlos con resultados obtenidos por Goody [8], como punto de partida para empezar nuevos estudios. De esta misma gráca se toman los valores del parámetro λ , que se relaciona con el de la gura de tal forma que λ=λ2 G . 45 46 Capítulo 4. Medio participativo dentro de una cavidad con paredes con absorción de radiación nula Figura 4.1: Gráca tomada como referencia para elegir los valores de λG=λ2 , donde λG es el parámetro que aparece en el eje horizontal. Los valores del número de Rayleigh crítico aquí mostrados, corresponden a los valores de la teoría lineal con paredes horizontales libres. Obsérvese también que se ha usado en esta gráca una doble escala logarítmica. Es también notorio que los valores del parámetro χ inuyen solo En el caso de radiación opaca, como es lógico después de lo visto en la teoría del capítulo 2 para el caso opaco. Se observa que los valores calculados siguen un comportamiento similar al que se produce en el lado izquierdo de las curvas, es decir, en el caso de radiación transparente, que es el de aumentar considerablemente el valor de Rac cuando λG≥10 . Imagen tomada de [8]. Así, los resultados de número de Rayleigh crítico en comparativa con los de Goody se muestran en la tabla 4.1: Cuadro 4.1: Tabla comparativa con los resultados de Goody de la gura 4.1. Si se tiene en cuenta la escala logarítmica se observa que para los valores del parámetro de λG tomados tienen el mismo comportamiento que los recogidos en la tabla, poniéndose de maniesto el aumento de la estabilidad del uido, ya que es capaz de emitir radiación a las paredes, estabilizándose con mayor rapidez. Asimismo, para una cavidad con paredes horizontales rígidas con temperatura impuesta y paredes verticales adiabáticas, para una relación de aspectos A= 2 , se muestra la inuencia del parámetro de radiación λ sobre el número de Rayleigh crítico para la misma cavidad en el problema de convección, cuyo resultado se 46 4.1. Paredes verticales adiabáticas 47 puede consultar en [1], que obtuvo un valor de Rac= 2014 , frente a los que se han obtenido aquí, a saber: Cuadro 4.2: Tabla comparativa del número de Rayleigh crítico con el caso de convección natural si se tiene en cuenta aproximación de medio transparente. Obsérvese de nuevo como al aumentar el parámetro λG , aumenta la estabilidad del uido y por tanto aumenta el número de Rayleigh crítico, siendo muy importante el aumento para números altos como λG= 10 ⇒λ= 100 . De esta manera, se reproduce el comportamiento que Goody obtuvo. 4.1.2. Variación del no de Rayleigh crítico con la relación de aspectos. Inuencia del medio En este apartado se analiza cómo varía el número de Rayleigh crítico cuando las dimensiones de la cavidad pasan desde una relación de aspectos muy pequeña (A=0.5) hasta una relación de aspectos alta. Intuitivamente, antes de mostrar los resultados se puede comentar qué se espera obtener. A priori, una menor relación de aspectos, debido a la condición de paredes verticales rígidas, diculta que se produzcan rollos de convección, por lo que el número de Rayleigh crítico a partir de el cual la convección sea posible, tendrá que aumentar, tanto más cuanto mayor sea la inuencia de la presencia de las paredes verticales, esto es, a medida que la relación de aspectos disminuya. Por otro lado, se espera que comparando para los distintos valores de λ , aumenten los números de Rayleigh críticos, como se ha visto en el apartado anterior. A posteriori, se incluyen comentarios sobre dichos resultados y se justican a partir de conclusiones y comparaciones con la teoría disponible. A continuación se muestran los resultados obtenidos para varios valores del parámetro λ denido en la aproximación ópticamente na, desde valores bajos en los que la radiación es inapreciable hasta valores altos , habiéndose usado la gráca de Goody (gura 4.1) como referencia. De [8] (página 434), teniendo en cuenta que χ >> 1 , se obtiene la relación entre el parámetro de la gura, λG y el del capítulo 2, λ : λ=λ2 G; (4.1) Según los resultados de Goody, el número de Rayleigh crítico aumenta a partir de λG∼1 por ser el valor a partir del cual la radiación comienza a tener peso cuando la relación de aspectos es grande (teóricamente innita). A continuación se muestra el número de Rayleigh crítico calculado en función de la relación de aspectos para varios valores de λG : Cuadro 4.3: Número de Rayleigh crítico en función de A para varios valores del parámetro λG . Cavidad con paredes rígidas. Las celdas de la tabla que aparecen en color azul signican que no se ha calculado el valor correspondiente, por considerarse innecesario e irrelevante, pues ya fueron calculados en [1]. 47 48 Capítulo 4. Medio participativo dentro de una cavidad con paredes con absorción de radiación nula Se observa en la tabla anterior que para un valor de λ dado, a medida que aumenta la relación de aspectos, el número de Rayleigh crítico disminuye, pues la inuencia de la condición de contorno de rigidez o velocidad nula en las paredes, se atenúa, debido a que la separación entre ellas va en aumento. Además, si se ja una relación de aspectos, a mayor valor de λ , mayor valor de Rac , a excepción del caso en el que la relación de aspectos es muy pequeña, A=0.5, donde dicho parámetro comienza a disminuir a partir del caso en el que se considera la participación del medio transparente, λ > 0 , alcanza un mínimo para λG= 5 aproximadamente, y luego vuelve a aumentar. Probablemente, este último comportamiento se pueda explicar como que en primer lugar el efecto del rozamiento de las paredes rígidas disminuye al poder transmitirse energía por otro medio, por lo que el comportamiento inicialmente con λ↑ es una disminución del número de Rayleigh crítico. En cambio, cuando el parámetro aumenta demasiado, la inuencia sobre el medio provoca su rápida homogeinización, de tal forma que sea necesario un número de Rayleigh mayor para provocar inestabilidades en el uido. En cualquier caso, para pequeñas relaciones de aspectos posiblemente se debería de tener en cuenta la participación de las paredes en la radiación, que es precisamente en lo que se centra el capítulo siguiente, el más importante del trabajo. 4.1.3. Número de Nusselt y perles de temperatura El número de Nusselt Nu sirve como medida de la relación entre el calor que se trasmite por convección (movimiento) y el calor que se transferiría si no existiese movimiento (equilibrio), por tanto sirve para generar una idea de cuánta importancia tiene, bajo determinados valores de los parámetros del problema, la convección a la hora de disipar calor. Por tanto, el número de Nusselt local NuL medido en la pared caliente ( x∗= 0 o z∗= 0 dependiendo del problema) responde a la siguiente denición matemática: NuL(x) = −K∂T ∂z z=0 −K∂Te ∂z z=0 ; (4.2) Entonces, usando las variables adimensionales expuestas en la sección Aproximación ópticamente na o transparente , se tiene:: NuL(x) = ∂(θe+∆¯ T Ra T∗) ∂z∗z∗=0 ∂θe ∂z∗z∗=0 . (4.3) Es conveniente también denir un número de Nusselt local medio, ¯ NuL , denido como el calor intercambiado, teniendo en cuenta la convección, por el uido a lo largo de la placa (por unidad de tiempo y de longitud y ) y el calor intercambiado en equilibrio, también por unidad de tiempo y de longitud de la placa según y , que da como resultado: ¯ NuL=RA 0 ∂(θe+∆¯ T Ra T∗) ∂z∗z∗=0 dx∗ RA 0 ∂θe ∂z∗z∗=0 dx∗, (4.4) Esta última denición es la que se expondrá en los resultados pare evitar un gran número de grácos, a los que se recurre solo en caso de que se estime oportuno. De esta manera, el número de Nusselt ha de variar con la inuencia de la radiación. Para estudiarlo, se ha calculado, con los valores de la gura 4.1 de Rac , el número de Nusselt medio en función de Ra/Rac , comparándolo con el caso en el que no hay radiación, habiéndose recreado los resultados del número de Nusselt local obtenido por [1]: 48 4.1. Paredes verticales adiabáticas 49 Figura 4.2: Número de Nusselt medio medido en la pared caliente, z=0, para el caso de convección natural, línea azul, ya resuelto en [1] y para los valores de λG= 1 , λG= 5 y λG= 10 . En la gura anterior 4.2, se observa que el número de Nusselt local disminuye cuando λG (o λ ) aumenta. Téngase en cuenta que se tienen dos efectos que son contrapuestos en cuanto a su inuencia en este número adimensional: por un lado, al aumentar el parámetro λ , el número de Rayleigh crítico aumenta, según se ha visto en la sección anterior, por lo que según esto, el calor disipado por movimientos convectivos, podría ser mayor (a mayor número de Rayleigh mayor intensidad de la convección); por otro lado, al aumentar el valor de λ , el medio cada vez tiene mayor capacidad de evacuar calor por radiación, por lo que el número de Nusselt debería disminuir. Consecuentemente, ha quedado claro, tras estos resultados, que el efecto homogeneizador de la radiación es superior al efecto de tener mayores valores de Rayleigh para cada punto calculado de las grácas (ya que lo que se representa en el eje horizontal es el valor Ra/Rac , no Ra). Se muestran a continuación los distintos perles de temperatura que se obtienen cuando cambia el parámetro de radiación λ . 49 50 Capítulo 4. Medio participativo dentro de una cavidad con paredes con absorción de radiación nula (a) (a) λG= 0 (b) (b) λG= 1 (c) (c) λG= 5 (d) (d) λG= 10 Figura 4.3: Perles de la temperatura y la de equilibrio, adimensionalizadas con la temperatura media, para varios valores de λ , jando el parámetro Ra/Rac= 2 para cada caso. En las guras 4.3, donde se ha jado la relación de aspectos en A=2, se observa en primer lugar que para z∗= 0 y z∗= 1 se tiene que Te/Tm= 0.95 y T/Tm= 1.05 en cualquiera de ellas, debido a la imposición de las condiciones de temperatura ( ∆¯ T se ha jado en 0.1). Asimismo, es notorio que medida que λ aumenta, la importancia del movimiento en la transmisión de calor disminuye, pues la temperatura total se hace cada vez más similar a la temperatura que existe solo en equilibrio, si no existiese movimiento. 4.2. Par baroclínico De forma introductoria previa a la exposición de resultados se va a explicar la diferencia fundamental existente en cuanto a la física, cuando las temperaturas están impuestas en paredes verticales y cuando están impuestas en las paredes horizontales. En ausencia de absorción de radiación por parte de las paredes, Rc=0, el gradiente de temperaturas en el segundo de los casos no depende de la coordenada x por lo que el gradiente de temperaturas es paralelo al de presión y por tanto el de densidad también sigue la misma dirección. En cambio, en el primero de los casos, las condiciones de contorno imponen valores constantes de temperatura T1 y T2 en las paredes verticales izquierda y derecha, respectivamente, cumpliéndose que T1>T2 . Esto provoca un gradiente de temperaturas horizontal, que al dejar de ser paralelo al gradiente de 50 4.3. Paredes horizontales adiabáticas 51 presiones, provoca que el gradiente de densidad tampoco sea paralelo al de presiones generándose un par baroclínico . Esto se explica con mayor grado de detalles a continuación. Debido al no paralelismo de los gradientes de presión y densidad, deducción que se obtiene de tomar el gradiente en la ecuación de estado: ρ=p RgT ⇒ ∇ρ=1 RgT ∇p−1 RgT2∇T, (4.5) pues ∇T∦∇p y al desplazarse el centro de masas de las partículas en la dirección del gradiente de densidad, por ser esta no uniforme, al tener en cuenta que las fuerzas de presiones, −∇pdΩ , sí están aplicadas en el centro de masas, se crea un momento con respecto a él que se conoce como par baroclínico (véase la gura 4.4). Figura 4.4: Par baroclínico, gura cortesía de [6] El lector interesado puede consultar [6] (sección 4.4) para más detalles. La parte esencial y relevante de este fenómeno es que en el problema que a continuación se presenta, no existe un número de Rayleigh crítico a partir del cuál se inicie el movimiento de convección, pues existe la inestabilidad expuesta en cualquier caso. 4.3. Paredes horizontales adiabáticas Como se ha explicado con anterioridad, en este caso no existe un número de Rayleigh crítico, por lo que la obtención de resultados se ha enfocado a grácas que representen los isocontornos del campo de temperaturas, al campo de velocidades y también se obtienen resultados para el número de Nusselt medio y en función de la coordenada z . . 4.3.1. Isocontornos de temperatura y función de corriente En este apartado se ha estudiado cómo afecta el parámetro λ a los campos de temperatura y de función de corriente (aquí se observa su efecto sobre los rollos de convección). 51 52 Capítulo 4. Medio participativo dentro de una cavidad con paredes con absorción de radiación nula (a) λG= 1 (b) λG= 4 (c) λG= 25 (d) λG= 100 Figura 4.5: Isocontornos de los campos de temperatura y función de corriente. 4.4. Resolución del caso opaco Mención aparte merece el caso en el que se considera la aproximación opaca . En el momento en el que se obtuvieron las ecuaciones se observó que las ecuaciones de Saltzmann válidas para la convección natural eran análogas a las que se obtenían para este caso si se reescribían como se hizo en 2.52. Por tanto se dijo entonces que el único cambio en el que se podrían diferenciar los valores de este caso con los del caso de convección natural, es que los valores de Rac estarían multiplicados por un factor de 1 + χ , donde χ=16σT3 m 3κk , por tanto este problema se presupone trivial. En cualquier caso, se ha corroborado que así sucede a través del siguiente resultado: en la gura 4.6 se ha muestra que el valor numérico calculado como número de Rayleigh crítico para esta aproximación, si se tienen paredes libres es de 671.3(1 + 103) con A=10, si χ= 103 . Este valor puede ser comparado con el valor teórico para una cavidad innita con paredes libres, que es 657.25 o con el valor numérico calculado en apartados posteriores. 52 4.4. Resolución del caso opaco 53 Figura 4.6: Ψ∗ e isocontornos de temperatura para paredes libres en el caso de aproximación opaca. En conclusión, se ha comprobado que el caso de aproximación ópticamente gruesa es extrapolable de los resultados del problema de convección natural, por lo de aquí en adelante se pondrá la atención en el caso de aproximación ópticamente na y su comparativa con el caso de radiación natural. De hecho, la longtidud crítica entre dos rollos convectivos o período espacial , λc resulta ser, a la vista de la gura 4.6: λc=7.143 −2.857 3/2= 2.857, que, comparándolo con el valor de la Teoría Lineal de Rayleigh de λc para una cavidad innita con paredes horizontales libres, λc= 2.828 , asumiendo que pueden existir pequeños errores numéricos, se refuerza la idea de que este caso tan solo diere en el factor (1 + χ) del número de Rayleigh. 53 60 Capítulo 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula (a) = 0.7 (b) = 0.85 (c) = 1 Figura 5.4: Isocontorno de temperaturas para valores altos de la emisividad: obsérvense los errores numéricos que se cometen en las paredes horizontales en la zona cercana a las esquinas cuando la emisividad es igual a la unidad. En la gura anterior se observa cómo la emisividad no altera en exceso el campo de temperaturas, excepto para la gura 5.4 donde = 1 en las zonas cercanas a las paredes verticales, donde en las esquinas en la pared izquierda se forman zonas con temperatura superior a la temperatura de la pared caliente, y zonas con temperatura inferior a la temperatura de la pared fría, en las esquinas de la derecha de la cavidad. Naturalmente, estos resultados no corresponden con la realidad, ni tampoco con los obtenidos por [24], por lo que se tratan de limitaciones del método numérico desarrollado. No obstante, se ha intentado a posteriori desarrollar un método eliminando la linealización que antes se realizaba, para comprobar que efectivamente era la causa raíz de los errores obtenidos numéricamente. En la siguiente sección ( ?? ) se tratará este problema, como mejora de los resultados hasta ahora reproducidos. 5.1.3. Temperatura en las paredes aisladas En esta sección se ha puesto interés en calcular la evolución de la temperatura en las paredes aisladas, para observar cómo afectan la emisividad y el número de Rayleigh, es decir la inuencia de la radiación y de la convección, a la distribución de temperatura que alcanzan en las paredes horizontales. Como 60 5.1. Paredes con absorción de radiación no nula. Resolución numérica del problema para una cavidad rectangular con paredes horizontales adiabáticas 61 referencia para estos resultados se pueden consultar las grácas de la página 193 de la referencia [3]. (a) Ra = 103 , = 0 (b) Ra = 7 ·104 , = 0 (c) Ra = 103 , = 0.5 (d) Ra = 7 ·104 , = 0.5 (e) Ra = 103 , = 1 (f) Ra = 7 ·104 , = 1 Figura 5.5: Perles de temperatura a lo largo de las paredes aisladas. 61 62 Capítulo 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula En primer lugar se observa que a mayor número de Rayleigh, mayor es la diferencia de temperatura entre las paredes aisladas, especialmente en el caso en el que no se considera radiación, lo cual quiere decir que el transporte de calor por convección calienta de forma notoria a la pared superior y enfría a la inferior cuando la recirculación aumenta (Ra aumenta), cuyas temperaturas solo se igualan en x∗= 0 y x∗= 1 por condición de contorno. De hecho, si no hubiese movimientos de convección (Ra=0), ambas paredes tendrían la misma temperatura. Por otro lado, al aumentar la emisividad para tener en cuenta las paredes, se halla que la consideración de la radiación facilita la equidad de las temperaturas entre las paredes aisladas, ya que pueden intercambiar calor entre sí directamente y no solo a través del uido. 5.2. Corrección de los resultados: condiciones de contorno sin linealizar En este apartado se expone el proceso seguido para el desarrollo de un nuevo código más preciso que el de la versión anterior, a n de comparar resultados para valores de emisividad altos, donde se encontraban fallos. La diferencia radica en que no se linealiza la temperatura, la cual se va a considerar como la suma de tres términos: T=Tm+T0+T+T0 e, como se ha hecho anteriormente. Cabe destacar que para el código se mantiene la notación de índices y el nombre de las matrices que se ha usado en la sección anterior. A continuación se explica cómo se han impuesto las condiciones de contorno en el programa para el cálculo de la temperatura de equilibrio (es decir, obviando la hidrodinámica), y de manera análoga se procede con la temperatura completa. De este modo, en el equilibrio, en primer lugar se tiene la condición que se ha de cumplir en las paredes con temperatura impuesta: ∂T ∂t −α∇2T= 0; (5.3) Para imponer la condición de contorno en una pared aislada hay que hacer un balance entre los calores intercambiados por convección con el uido y el calor neto intercambiado por radiación, de tal modo que: K∂T ∂z ∓Qr= 0; paraz=(0 H (5.4) Se ha optado por introducir las variables adimensionales x∗ , z∗ y t∗ ya presentadas previamente (se omitirán los superíndices por comodidad), además de una variable adimensional de temperatura de equilibrio θE : T=Tm+ (T1−T2)θE; donde Tm=T1+T2 2, (5.5) cuyo subíndice E se omite en adelante por comodidad, de tal modo que se denota θ≡θE . Así, la ecuación 5.3 quedaría: ∂θ ∂t − ∇2θ= 0 Discretiz. =====⇒Tn−Tn−1 dt −DL ∗Tn= 0 ⇒eye(Nt)−dt ∗DL)Tn=Tn−1, (5.6) donde el subíndice n representa la iteración actual en el tiempo y n−1 , la anterior, que es conocida y a la matriz eye(Nt)−dt ∗DL se le denota en adelante y en el código por BT E . Para la ecuación para pared aislada 5.4, usando el desarrollo completo de la expresión (1 + θ/θ0)4 , donde θ0= 1/∆¯ T y teniendo en cuenta que el calor que recibiría un elemento de pared sería nulo si todos los 62 5.2. Corrección de los resultados: condiciones de contorno sin linealizar 63 elementos estuviesen a la misma temperatura, condición que se expresa como: Nwall X j=1 AQicjjσ·1=0, (5.7) se tiene que la ecuación 5.4 se convierte en: 1∂θ ∂z ic −nic Rcθ0 4 Nwall X j=1 AQicjjgjθn,donde nic=(1si z = 0 −1si z =zmax , (5.8) donde gj es la aproximación numérica 4/θ0+ 6θn−1/θ2 0+ 4θ2 n−1/θ3 0+θ3 n−1/θ4 0 , de tal manera que para cada iteración en el tiempo represente un valor conocido, y que se ha obtenido de una aproximación del desarrollo de (1 + θ/θ0)4 y usar 5.7. Para implementar las condiciones de contorno de temperatura impuesta (obsérvese que corresponden con θ= 1/2 en x= 0 y θ=−1/2 en x=xmax ) es necesario imponer en la matriz del sistema de la temperatura de equilibrio BT E , como se ha hecho en repetidas ocasiones, los siguientes cambios: Teniendo el nodo de la pared identicado con su índice I , se anula dicha la I completamente. El elemento BT E(I, I) se hace igual a la unidad. Por último, en el elemento I del vector que almacena los términos independientes (llamado bTE ) se asigna el valor de ±1/2 según si corresponde a uno nodo de la pared caliente, x= 0 , o fría, x= 1 . Para implementar la condición de contorno dada por 5.8, que es la parte más compleja e interesante, se siguen los pasos siguientes: Siguiendo la notación usada en el capítulo de Método numérico, se denen varios índices para hacer referencias a los nodos y a los paneles de los contornos, a saber: un índice kw para paneles con condición de calor impuesta (paredes horizontales en este caso) tal que empieza numerando como primer panel al panel que se sitúa entre el primer nodo de la esquina superior izquierda y el punto medio entre los dos nodos siguientes. A partir de ahí, en la pared superior recorre de izquierda a derecha (es decir los primeros Nx-2 paneles) y continúa en la pared inferior de derecha a izquierda (numerando por tanto los Nx-2 paneles de esta pared). Adicionalmente, a cada panel kw le corresponde un índice ic asociado a la enumeración original de los paneles usada en el programa previo (empieza numerando por el contorno izquierdo en la zona inferior y en sentido horario desde ese primer panel). Además, para recorrer el resto de paneles habiendo jado uno ic , se utiliza un índice j que recorre todos los paneles de los contornos, ya sean paneles con temperatura impuesta o aislados, por lo que j=1:Nwall , siendo Nwall = 2 ∗Nx + 2 ∗Nz −8 . Por último, cada índice I denido al principio del capítulo de Método numérico está relacionado con un par (ic, j) . Una vez explicados los índices que se han usado, se exponen las matrices y vectores creados con el n de exponer la condición de contorno de pared aislada sin linealizar las temperaturas. En primer lugar se crea una matriz auxiliar BT EauxBC para almacenar los valores del segundo término de la ecuación 5.8, es decir, para cada nodo de las paredes aisladas se almacena en una la de Nt columnas los valores −normalH(kw)Rc/4∗AQ(ic, j)(j) , donde normalH(kw) tiene un valor de 1 si 1≤kw ≤Nx −2 o un valor de -1 si Nx −1≤kw ≤2∗Nx −4 y la matriz AQ es la misma que la del apartado 3.2. Después, se crea un vector auxiliar de Nt componentes, de tal modo que la componente J (correspondiente a la columna J de la matriz BTEauxBC) de dicho vector es 1 si J es un nodo de los contornos y cero si no pertenece a los contornos. Este vector sirve para tener en cuenta que sólo transmiten calor, según el término −nic Rcθ0 4PNwall j=1 AQicjjσgjθn , los paneles que pertenecen al contorno. Para nalizar, téngase en cuenta que debido a que la condición de contorno en los paneles aislados depende de la temperatura cambiante de los paneles de la otra pared aislada, por lo que la condición de contorno cambia cada vez que se produce una iteración. Esto lleva a recalcular lo que se ha llamado auxv en el programa, que contiene los valores de la expresión de gj en los nodos que pertenecen a paneles del contorno. Debido a esto, la matriz del sistema, BT E , cambia con cada iteración del tiempo, por lo que hay que invertirla para resolverlo en cada iteración. 63 64 Capítulo 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula Para el cálculo de la temperatura T∗ las líneas de código necesarias son prácticamente análogas a las que se han necesitado para el cálculo de la temperatura de equilibrio. Del mismo modo, se puede aprovechar el nuevo programa para la obtención de resultados en el caso en el que las paredes aisladas son las verticales. Simplemente Para obtener la temperatura con la hidrodinámica hay que tener en cuenta que en el desarrollo de la ecuación 5.8 tiene que tomarse su expresión completa: T=Tm1 + θE θ0 +T∗ Raθ0, (5.9) quedando así la ecuación de la condición de contorno denitivamente: ∂T∗ ∂z (n)−nic Rcθ0 4 Nwall X j=1 AQijjT∗(n) "41 + θE θ031 θ0 + 6 1 + θE θ02T∗(n−1) Raθ2 0 +1 + θE θ0T∗2(n−1) Ra2θ3 0 +T∗3 Ra3θ4 0#= 0, (5.10) donde se ha desarrollado la expresión de T4 y se ha usado la condición de equilibrio 5.8. 5.2.1. Código mejorado En este apartado se expone la versión del código sin uso de linealización para el caso en el que las temperaturas están impuestas en las paredes verticales. 1 2 3 clear all; 4 close all; 5 clc; 6 % _____PARAMETRES_____ % 7 Pr=0.73; 8 Ra=1e5; theta0=29.35; DeltaTb=1/theta0; 9 10 Tm=293.5; 11 DeltaT=10; %NOT DeltaTb 12 beta_nu_alpha=3.7836e+05; %let's keep this as a constant 13 d=(Ra/DeltaT/beta_nu_alpha)^(1/3); %as a function of the rest of the parameters. This value is needed to calculate RaAki 14 sigma=5.67e−8; %constant 15 K=0.0257; %conductivity at Tm 16 emissivity=1; 17 18 RcAki=Tm^3*sigma*d*theta0/K; 19 20 if emissivity>0 21 Rc=4*RcAki*DeltaTb/emissivity; 22 else 23 Rc=0; 24 end 25 26 27 dt=0.0001; 28 % ______GEOMETRY_____ % 29 30 zmin=0; zmax=1; 31 %Chebyshev nodes 64 5.2. Corrección de los resultados: condiciones de contorno sin linealizar 65 32 Nz=30; 33 zch(1:Nz)=(zmax+zmin)/2−(zmax−zmin)/2*cos(((1:Nz)−1)*pi/(Nz−1)); 34 35 xmin=0; xmax=1; 36 %Chebyshev nodes 37 Nx=30; 38 xch(1:Nx)=(xmax+xmin)/2−(xmax−xmin)/2*cos(((1:Nx)−1)*pi/(Nx−1)); 39 tic 40 % _____MATRICES____ % 41 [lpx]=dCheby(xch,Nx); 42 [lpz]=dCheby(zch,Nz); 43 lx=eye(Nx); 44 lz=eye(Nz); 45 Nt=Nz*Nx; 46 47 for i=1:Nx 48 for j=1:Nz 49 I=(i−1)*Nz+j; 50 X(I)=xch(i); Z(I)=zch(j); 51 for m=1:Nx, 52 Ki=(m−1)*Nz+1; 53 Kf=(m−1)*Nz+Nz; 54 Dx(I,Ki:Kf)=lpx(i,m)*lz(j,1:Nz); 55 Dz(I,Ki:Kf)=lx(i,m)*lpz(j,1:Nz); 56 end 57 end 58 end 59 % 60 % 61 Nwall=2*Nx+2*Nz−8; 62 kw1L=1; kw2L=Nz−2; 63 kw1U=Nz−1; kw2U=Nx+Nz−4; 64 kw1R=Nx+Nz−3; kw2R=Nx+2*Nz−6; 65 kw1B=Nx+2*Nz−5; kw2B=2*Nx+2*Nz−8; 66 for kw=1:Nwall, 67 if kw>=kw1L && kw<=kw2L, 68 Iwall(kw)=kw+1; Iw1=Iwall(kw)−1; Iw2=Iwall(kw)+1; 69 end 70 if kw>=kw1U && kw<=kw2U,, 71 Iwall(kw)=(kw−kw1U+2)*Nz; Iw1=Iwall(kw)−Nz; Iw2=Iwall(kw)+Nz; 72 end 73 if kw>=kw1R && kw<=kw2R, 74 Iwall(kw)=Nx*Nz−(kw−kw1R+1); Iw1=Iwall(kw)+1; Iw2=Iwall(kw)−1; 75 end 76 if kw>=kw1B && kw<=kw2B 77 Iwall(kw)=(Nx−1)*Nz+1−(kw−kw1B+1)*Nz; Iw1=Iwall(kw)+Nz; Iw2=Iwall(kw)−Nz; 78 79 end 80 Iw=Iwall(kw); 81 xwall(kw,1)=(X(Iw1)+X(Iw))/2; xwall(kw,2)=(X(Iw2)+X(Iw))/2; 82 zwall(kw,1)=(Z(Iw1)+Z(Iw))/2; zwall(kw,2)=(Z(Iw2)+Z(Iw))/2; 83 if (kw−kw1L)*(kw−kw1U)*(kw−kw1R)*(kw−kw1B)==0, 84 xwall(kw,1)=X(Iw1); 85 zwall(kw,1)=Z(Iw1); 86 end 87 if (kw−kw2L)*(kw−kw2U)*(kw−kw2R)*(kw−kw2B)==0, 88 xwall(kw,2)=X(Iw2); 89 zwall(kw,2)=Z(Iw2); 90 end 91 % [kw Iw Iw1 Iw2] 92 % [kw Iwall(kw) Iw1 Iw2 xwall(kw,1) xwall(kw,2)] 65 66 Capítulo 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula 93 %pause 94 end 95 % 96 % Emissivities and View factors by Hottel's rule 97 % 98 %emiss(1:Nwall)=0.7; 99 emiss(1:Nwall)=emissivity; 100 for kw=1:Nwall, 101 x1A=xwall(kw,1); x2A=xwall(kw,2); 102 z1A=zwall(kw,1); z2A=zwall(kw,2); 103 Long(kw)=sqrt((x2A−x1A)^2+(z2A−z1A)^2); 104 for mw=1:Nwall, 105 x1B=xwall(mw,1); x2B=xwall(mw,2); 106 z1B=zwall(mw,1); z2B=zwall(mw,2); 107 L1A1B=sqrt((x1A−x1B)^2+(z1A−z1B)^2); 108 L2A2B=sqrt((x2A−x2B)^2+(z2A−z2B)^2); 109 L2A1B=sqrt((x2A−x1B)^2+(z2A−z1B)^2); 110 L1A2B=sqrt((x1A−x2B)^2+(z1A−z2B)^2); 111 Fview(kw,mw)=(L1A1B+L2A2B−L1A2B−L2A1B)/2/Long(kw); 112 end 113 Fview(kw,kw)=0; 114 end 115 Ident=eye(Nwall); 116 for kw=1:Nwall, 117 Femiss(kw,1:Nwall)=emiss(kw)*Fview(kw,:); 118 Arad(kw,1:Nwall)=Ident(kw,:)−(1−emiss(kw))*Fview(kw,:); 119 end 120 Aradm1=inv(Arad); 121 AQ=Ident−Femiss*Aradm1; 122 % 123 DL=Dx*Dx+Dz*Dz; 124 DL2=DL*DL; 125 DxDL=Dx*DL; 126 DzDL=Dz*DL; 127 % 128 % Preparing BC's: 129 NwallH=2*(Nx−2); 130 NwallT=2*(Nz−2); 131 kwallT(1:NwallH)=[kw1L:kw2L, kw1R:kw2R]; 132 nwallT(1:NwallH)=[ones(1,kw2L−kw1L+1), −ones(1,kw2R−kw1R+1)]; 133 kwallH(1:NwallT)=[kw1U:kw2U, kw1B:kw2B]; 134 nwallH(1:NwallH)=−[ones(1,kw2U−kw1U+1), −ones(1,kw2B−kw1B+1)]; 135 BTauxBC=zeros(NwallH,Nt); TauxBC=zeros(1,Nt); normalH=zeros(Nt,Nt); 136 for kw=1:NwallH; 137 ic=kwallH(kw); Ic=Iwall(ic); IwH(kw)=Ic; 138 for mw=1:NwallH 139 jc=kwallH(mw); Jc=Iwall(jc); 140 BTauxBC(kw,Jc)=−nwallH(kw)*Rc*theta0/4*AQ(ic,jc)*emiss(jc); 141 TauxBC(1,Jc)=1; 142 end 143 for nw=1:NwallT, 144 jT=kwallT(nw); JT=Iwall(jT); 145 BTauxBC(kw,JT)=−nwallH(kw)*Rc*theta0/4*AQ(ic,jT)*emiss(jT); 146 TauxBC(1,JT)=1; 147 end 148 end 149 %EQUILIBRIUM 150 % 151 dte=10*dt; 152 BTE=eye(Nt)−dte*DL; 153 bTE=zeros(Nt,1); 66 5.2. Corrección de los resultados: condiciones de contorno sin linealizar 67 154 fBC(1:Nt,1)=0; fInt(1:Nt,1)=1; 155 % 156 157 for j=1:Nz, 158 I1=j; 159 I2=(Nx−1)*Nz+j; 160 161 BTE(I1,:)=0; BTE(I1,I1)=1; bTE(I1,1)=1/2; fBC(I1,1)=1; fInt(I1,1)=0; 162 BTE(I2,:)=0; BTE(I2,I2)=1; bTE(I2,1)=−1/2; fBC(I2,1)=1; fInt(I2,1)=0; 163 end 164 165 theta0=1/DeltaTb; theta02=theta0*theta0; theta03=theta02*theta0; theta04=theta03*theta0; 166 Thetanm1(1:Nt,1)=−1/2; 167 for it=1:10000, 168 [it] 169 vaux1=TauxBC.*Thetanm1'; 170 vaux2=vaux1.*vaux1; 171 vaux3=vaux2.*vaux1; 172 auxv=4/theta0+6*vaux1/theta02+4*vaux2/theta03+vaux3/theta04; 173 for kw=1:NwallH 174 Ikw=IwH(kw); 175 BTE(Ikw,:)=Dz(Ikw,:)+BTauxBC(kw,:).*auxv(1,:); 176 fBC(Ikw,1)=1; fInt(Ikw,1)=0; 177 end 178 % BTEinv=inv(BTE); 179 % Thetan=BTEinv*(fBC.*bTE+fInt.*Thetanm1); 180 AUX=fBC.*bTE+fInt.*Thetanm1; 181 Thetan=BTE\AUX; 182 for j=1:Nz, 183 xmat(1:Nx,j)=xch(1:Nx)'; zmat(1:Nx,j)=zch(j); 184 Thetanm1mat(1:Nx,j)=Thetanm1(((1:Nx)−1)*Nz+j,1); 185 end 186 mesh(xmat,zmat,Thetanm1mat) 187 pause(0.001) 188 hold off 189 if max(abs(Thetanm1−Thetan)) <= 10^−5, 190 break 191 end 192 Thetanm1=Thetan; 193 end 194 [max(abs(Thetanm1−Thetan))] 195 ThetaE=Thetan; 196 pause 197 198 Thetae=DeltaTb*ThetaE; 199 DxThetae=Dx*Thetae; DzThetae=Dz*Thetae; 200 for j=1:Nz, 201 Thetaemat(1:Nx,j)=Thetae(((1:Nx)−1)*Nz+j,1); 202 end 203 for I=1:Nt, 204 aux1(I,:)=Dz(I,:)*DxThetae(I,1); 205 aux2(I,:)=Dx(I,:)*DzThetae(I,1); 206 end 207 % 208 AuxThetae=1+Thetae'; 209 AuxThetae2=AuxThetae.*AuxThetae; 210 AuxThetae3=AuxThetae2.*AuxThetae; 211 212 % _____BOUNDARY CONDITIONS_____ % 213 214 APSI=DL−dt*Pr*DL2; 67 68 Capítulo 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula 215 AT=dt*Pr*Dx; 216 BPSI=Ra/DeltaTb*dt*(aux1−aux2); 217 BT=eye(Nt)−dt*DL; 218 %Boundary conditions factors 219 FBC_PSI=ones(Nt,1); 220 FBC_T=ones(Nt,1); 221 222 for i=2:(Nx−1) 223 I=(i−1)*Nz+1; 224 APSI(I,:)=0; APSI(I,I)=1; FBC_PSI(I,1)=0; AT(I,:)=0; 225 I=(i−1)*Nz+Nz; 226 APSI(I,:)=0; APSI(I,I)=1; FBC_PSI(I,1)=0; AT(I,:)=0; 227 end 228 229 for j=1:Nz 230 I=j; 231 APSI(I,:)=0; APSI(I,I)=1; FBC_PSI(I,1)=0; AT(I,:)=0; 232 I=(Nx−1)*Nz+j; 233 APSI(I,:)=0; APSI(I,I)=1; FBC_PSI(I,1)=0; AT(I,:)=0; 234 end 235 236 for i=3:(Nx−2) 237 K=(i−1)*Nz+2; 238 APSI(K,:)=Dz((K−1),:); FBC_PSI(K,1)=0; AT(K,:)=0; 239 K=(i−1)*Nz+Nz−1; 240 APSI(K,:)=Dz((K+1),:); FBC_PSI(K,1)=0; AT(K,:)=0; 241 end 242 243 for j=2:(Nz−1) 244 K=(2−1)*Nz+j; 245 APSI(K,:)=Dx((K−Nz),:); FBC_PSI(K,1)=0; AT(K,:)=0; 246 K=(Nx−2)*Nz+j; 247 APSI(K,:)=Dx((K+Nz),:); FBC_PSI(K,1)=0; AT(K,:)=0; 248 end 249 % 250 251 for i=1:Nx, 252 for j=1:Nz, 253 I=(i−1)*Nz+j; 254 Psi(I,1)=0.05*xch(i)^2*(xmax−xch(i))^2*zch(j)^2*(zmax−zch(j))^2; 255 T(I,1)=0.05*sin(3*pi*zch(j)); 256 end 257 end 258 AbPsi=sparse([DL, zeros(Nt,Nt)]); 259 AbT=sparse([zeros(Nt,Nt), speye(Nt)]); 260 261 [X, Z]=meshgrid(xch,zch); 262 % 263 figure 264 Ntime=5000; 265 for nt=1:Ntime, 266 t=nt*dt; 267 % 268 %BPSI y BT: 269 for j=1:Nz 270 I=(1−1)*Nz+j; 271 BPSI(I,:)=0; BT(I,:)=0; BT(I,I)=1; FBC_T(I,1)=0; 272 I=(Nx−1)*Nz+j; 273 BPSI(I,:)=0; BT(I,:)=0; BT(I,I)=1; FBC_T(I,1)=0; 274 end 275 vaux1=TauxBC.*T'; 68 5.2. Corrección de los resultados: condiciones de contorno sin linealizar 69 276 vaux2=vaux1.*vaux1; 277 vaux3=vaux2.*vaux1; 278 auxv=4*AuxThetae3/theta0+6*AuxThetae2.*vaux1/theta02/Ra+4*vaux2.*AuxThetae/Ra^2/theta03+ vaux3/Ra^3/theta04; 279 for kw=1:NwallH 280 Ikw=IwH(kw); 281 BT(Ikw,:)=Dz(Ikw,:)+BTauxBC(kw,:).*auxv(1,:); 282 FBCT(Ikw,1)=1; 283 end 284 Asyst=[APSI AT; BPSI BT]; 285 bNLPsi(1:Nt,1)=(−(Dz*Psi).*(DxDL*Psi)+(Dx*Psi).*(DzDL*Psi))*dt−Pr*Ra/DeltaTb*dt* DxThetae; 286 bNLT(1:Nt,1)=(−(Dz*Psi).*(Dx*T)+(Dx*Psi).*(Dz*T))*dt; 287 bLinPsi(1:Nt,1)=AbPsi*[Psi ; T]; 288 bLinT(1:Nt,1)=AbT*[Psi ; T]; 289 bsyst(1:(2*Nt),1)=[FBC_PSI.*(bNLPsi+bLinPsi); FBC_T.*(bNLT+bLinT)]; 290 bn=Asyst\bsyst; Psi=bn(1:Nt,1); T=bn((Nt+1):(2*Nt),1); 291 for j=1:Nz, 292 Psimat(1:Nx,j)=Psi(((1:Nx)−1)*Nz+j,1); Tmat(1:Nx,j)=T(((1:Nx)−1)*Nz+j,1); 293 Tphysmat(1:Nx,j)=−1/2+Thetaemat(1:Nx,j)/DeltaTb+Tmat(1:Nx,j)/Ra; 294 end 295 subplot (2,2,1) 296 contour (X, Z, Psimat') 297 subplot (2,2,2) 298 contour (X, Z, Tphysmat') 299 Tmat_c(nt)=Tmat(ceil(Nx/2),ceil(Nz/2)); 300 Tmat_max=max(abs(Tmat_c(1:nt))); 301 subplot (2,2,3) 302 plot((1:nt)*dt,Tmat_c(1:nt)/Tmat_max) 303 axis([0 Nt*dt −2 2]) 304 subplot (2,2,4) 305 plot (Tphysmat(ceil(Nx/2),:),zch,'r') 306 axis([−2 2 0 1]) 307 pause(0.01) 308 hold off 309 end 5.2.2. Temperatura de equilibrio A continuación se presenta la temperatura de equilibrio, ya calculada con el programa anterior, para varios valores de los parámetros Rc y  , a n de estudiar su inuencia: 69 76 Capítulo 5. Medio no participativo dentro de una cavidad con paredes con absorción de radiación no nula (a) = 0 (b) = 0 Figura 5.13: Perl de temperaturas para paredes con emisividad nula y unidad. 76 Capítulo 6 Conclusiones y líneas de desarrollo En este trabajo se ha estudiado la inuencia de la radiación en el movimiento de convección de Rayleigh-Bénard, con aproximación de medio transparente (principalmente) y opaco, lo cual supone un tratamiento simplicado del medio participativo, el gas encerrado en la cavidad. Se han considerado dos paredes adiabáticas, lo cual conlleva la necesidad de modelar las condiciones de contorno, que se han implementado mediante el uso de un parámetro de radiación-convección y considerando la capacidad de emitir radiación por parte de las paredes, incluso se puede asignar una emisividad diferente a cada uno de los paneles de las paredes (estudio que aquí no se ha hecho, tomando todos los paneles del contorno con una misma emisividad) y dos paredes con temperatura impuesta. La motivación principal de este trabajo, cuyos resultados han hecho hincapié en el problema (el capítulo 5 es sin duda el capítulo más importante del tema) en el que un medio no es participativo pero las paredes de la cavidad que lo encierran sí que intercambian calor por radiación, se debe a su aplicación a hornos industriales donde el medio (aire) al estar compuesto principalmente por moléculas no polares se puede considerar no participativo, mientras que sus paredes sí que intercambian calor. De hecho, una tesis doctoral reciente (2007) [26] está centrada en el modelado de la radiación en hornos y cámaras de combustión. Asimismo, en acondicionamiento de habitaciones de edicios también se encuentra una situación similar (ver la cita [4], donde se estudia este tema). De hecho, para estudiar un caso de este estilo, como futura mejora, se podría dividir la cavidad estudiada en varios compartimentos simulando salas. Por otro lado, el otro problema estudiado, el caso de medio participativo, se puede dar por ejemplo en el aire, que se puede considerar un medio participativo siempre y cuando las distancias que ha de recorrer en él son muy grandes, como es el caso de la atmósfera. En cuanto al método de resolución del problema se ha utilizado el método de colocación con polinomios de Lagrange usando nodos de Chebyshev, lo cual facilita la convergencia de los resultados al aumentar el número de nodos. De esta forma, con códigos que podrían desarrollar alumnos que tengan conocimientos de MATLAB se demuestra que es posible calcular resultados que hace tan solo unas décadas no había posibilidad de obtener, salvo con cálculos engorrosos y usando demasiadas simplicaciones si se quería llegar a alguna solución de un problema tan complejo. Gracias a esos códigos del capítulo 4, se han reproducido en primer lugar y con alta precisión los resultados de Goody de estabilidad lineal, calculándose los números de Rayleigh críticos y posibilitándose también el estudio de estabilidad no lineal con los códigos implementados. Además, se ha estudiado cómo afecta el parámetro de radiación λ al número de Nusselt denido en la pared caliente y también cómo varía este para un valor nulo de dicho parámetro (medio totalmente transparente) pero con paredes que pueden intercambiar calor por radiación entre sí, es decir, para valores no nulos del parámetro de radiación-convección y de la emisividad de las paredes. Es importante notar también que, a diferencia de las citas aportadas como bibliografía que se centran en los resultados del régimen estacionario, con estos códigos se puede estudiar la evolución temporal de la variable en la que se tenga interés. Como posibles líneas de ampliación de este proyecto se tiene un gran abanico de posibilidades. Quizás la extensión más evidente es la reproducción del mismo problema pero de manera tridimensional, pero hay otras muchas que se podrían aplicar. Por ejemplo, se puede ampliar el estudio para medios que no se puedan aproximar ni como transparentes ni como opacos. A su vez, hay casos que debido a las aproximaciones realizadas sobre las temperaturas en las paredes, cuya diferencia es pequeña, es una hipótesis que no puede aplicarse en todos los problemas: por ejemplo, en atmósfera de estrellas, el gas se puede considerar como un 77 78 Capítulo 6. Conclusiones y líneas de desarrollo medio opaco, pero la diferencias de temperaturas es muy destacable, por lo que habría que usar otro tipo de hipótesis menos simplicativas para estudiarlo. También sería posible ahondar más en algunos de los problemas resueltos y modicar las condiciones de contorno de pared rígida o libre en ellos o estudiar cavidades con distintas relaciones de aspectos aprovechando los mismos códigos que aquí se han desarrollado. Téngase en cuenta que el objetivo no era el de aportar un sinfín de resultados sin explicar, si no solo estudiar algunos parámetros de interés como el número de Nusselt o los campos de temperatura y velocidad o función de corriente, como se ha hecho. Otra propuesta sería la posibilidad de extender el estudio de los problema para valores muy altos del número de Rayleigh, de tal modo que se alcance un régimen turbulento (se necesitaría también desarrollar un modelo tridimensional) ya que en los resultados obtenidos para Ra ≥106 se pierde precisión, como se vio en el capítulo 5. Para nalizar, también se podría implementar la división de la cavidad en varios compartimentos para estudiar la transferencia de calor entre las paredes de un edicio, como se ha descrito en párrafos previos. 78 Bibliografía [1] El Método de Colocación para el problema de convección de Rayleigh-Bénard , autor: Pablo José Ruiz Contreras, tutor: Miguel Pérez-Saborid Sánchez-Pastor, 2013. [2] Combined radiation and natural convection in a two-dimensional participating square medium , ZHIQIANG TAN and JOHN R. HOWELL , 1990. [3] RADIATION HEAT TRANSFER: FUNDAMENTALS AND APPLICATIONS: ANALYSIS OF RADIATION-NATURAL CONVECTION INTERACTIONS IN 1-G AND LOW-G ENVIRONMENTS USING THE DISCRETE EXCHANGE FACTOR METHOD , 1990 [4] Heat conduction in two and three dimensions : computer modelling of building physics applications Byggnadsfysik LTH, Lunds Tekniska Högskola , Blomberg Thomas , 1996 [5] Investigation of heat loss from a solar cavity receiver , E. Abbasi-Shavazi & G.O. Hughes & J.D. Pye , 2014 [6] Fundamentos y Aplicaciones de la Mecánica de Fluidos , autor: Antonio Barrero Ripoll & Miguel Pérez-Saborid Sánchez-Pastor [7] Procesos de convección natural con hipótesis anelástica , autor: Eduardo M. García Juárez , tutor: Miguel Pérez-Saborid Sánchez-Pastor, 2014 [8] The inuence of radiative transfer on cellular convection , autor: R. M. Goody, 1956. [9] Convection Heat Transfer, 4th Edition , cap. 1, Adrian Bejan , 2014. [10] Heat Transfer , cap. 4-5, Gregory Nellis & Sanford Klein , 2008. [11] Atmospheric Convection , Kerry A. Emanuel , 1994. [12] https://www.nasa.gov/ [13] , https://www.thermaluidscentral.org/encyclopedia/index.php/Properties_of_participating_media, references: Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO. [14] https://imagine.gsfc.nasa.gov/science/objects/sun1.html [15] , V. Basic Methods in Transfer Problems . Oxford University Press, 1952. Kourgano [16] On maintained convective motion in a uid heated from below , Anne Pellew & R. V. Southwell , 1940 [17] Multi Variable Calculus and Linear Algebra, with Applications to Dierential Equations and Probability , vol. II, John Wiley & Sons , 1969 by Xerox Corporation. [18] Introducción al Método de Diferencias Finitas y su Implementación Computacional , autor: Antonio Carrillo Ledesma y Omar Mendoza Bernal, Facultad de Ciencias UNAM, 2015. [19] http://mathworld.wolfram.com/KroneckerProduct.html [20] Numerical Methods Using MATLAB, fourth edition , cap. 4, John H. Mathews & Kurtis D. Fink [21] Heat and Mass Transfer: A Practical Approach , Yunus A. Cengel 3rd edition, 2006. 79 80 Bibliografía [22] Thermal Radiation Heat Transfer, 5th Edition , cap. 5, John R. Howell, Robert Siegel & M.Pinar Mengüç , 2010. [23] E. H. Ridouane , M. Hasnaoui , A. Amahmid & A. Raji (2004) INTERACTION BETWEEN NATURAL CONVECTION AND RADIATION IN A SQUARE CAVITY HEATED FROM BELOW, Numerical Heat Transfer, Part A: Applications, 45:3, 289-311, DOI: 10.1080/10407780490250373 [24] M. Akiyama & Q. P. Chong (1997) NUMERICAL ANALYSIS OF NATURAL CONVECTION WITH SURFACE RADIATION IN A SQUARE ENCLOSURE, Numerical Heat Transfer, Part A Applications, 32:4, 419-433, DOI: 10.1080/10407789708913899 [25] http://www.aerospaceweb.org/design/scripts/atmosphere/ [26] http://lup.lub.lu.se/record/548795 80