Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre

Equation Chapter 1 Section 1 Proyecto Fin de Grado Ingeniería Aeroespacial Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Antonio Enríquez López Tutor: Miguel Pérez-Saborid Sánchez-Pastor Sevilla, 2020 ii Proyecto Fin de Carrera Ingeniería Aeroespacial Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Autor: Antonio Enríquez López Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor titular Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 iv Proyecto Fin de Carrera: Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Autor: Antonio Enríquez López Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2020 El Secretario del Tribunal vi Agradecimientos Quisiera agradecer todos los desafíos que la Escuela me ha planteado en mis años de grado, porque me han ayudado a enfocar mis esfuerzos de manera eficiente. Además, son varios los docentes en los que he hallado, no solo personas apasionadas en su materia, sino también buenos interlocutores con los alumnos, que, en beneficio de éstos, procuran mejorar sus asignaturas, destacando, entre otros, los responsables de Cálculo de Aeronaves, Ingeniería Electrónica o Aviónica. A ellos les doy las gracias. También le agradezco a mi tutor, Miguel Pérez Saborid, la ayuda que me prestó para la asimilación de los conceptos involucrados en este proyecto, y la empatía que me demostró con el paso de las semanas, tanto durante las citas presenciales que durante los primeros meses mantuvimos, como durante las videoconferencias y los mensajes que intercambiamos durante los meses de confinamiento y de restricciones a la movilidad debido a la pandemia por Covid-19. Para terminar, y siendo uno de los miles de estudiantes que han hecho y defendido su trabajo fin de grado en este año marcado por la alerta sanitaria, quiero que conste en acta mi agradecimiento sincero al esfuerzo titánico que los servicios esenciales han realizado, realizan, y seguirán realizando para protegernos del Covid-19, y asegurar que la sociedad no se paralice en esta crisis que en septiembre de 2020 seguimos afrontando. También le doy las gracias a Felipe, una nueva amistad que trasciende a lo espiritual, que me ha dado esperanza, así como a mi familia y amigos, que los siento cerca de mi, ahora y siempre, y, por último, pero no menos importantes, a mis padres, por creer en mi y motivarme para realizarme como persona y sacar lo mejor de mi. viii Resumen Este proyecto presenta, a grandes rasgos, dos partes bien diferenciadas. En primer lugar, se presenta una introducción a modo de fundamentación teórica, cuyo objetivo es, por un lado, desarrollar los conceptos que involucran a las capas límite, generadas por convección forzada (flujo laminar subsónico ideal y viscoso confinado entre dos placas horizontales que se comportan como cuerpos negros) y por convección libre (flujo inducido laminar viscoso en placa plana vertical), y a la radiación térmica, prestando atención a las hipótesis de medio con espesor óptico delgado, y medio con espesor óptico grueso; y por otro, poner en contexto las ecuaciones a resuelver, partiendo del sistema de ecuaciones diferenciales parciales de Navier Stokes, así como de las expresiones del flujo de calor por radiación. En segundo lugar, una vez puestas en contexto las ecuaciones involucradas, se desarrolla sucintamente los métodos numéricos empleados para resolver los perfiles de velocidad y de temperatura de las capas límite en estudio, y se procede al despliegue de una serie de parámetros adimensionales cuyo sentido físico se alumbra. A continuación, partiendo de unos códigos originales, se encara el estudio de una serie de resultados gráficos, la mayor parte de ellos sobre una placa plana, obtenidos con MATLAB R2020a, que, en ocasiones, están acompañados de resultados de la literatura para demostrar su aplicabilidad. Una muestra amplia de los códigos empleados para su obtención se añade en un anexo para su consulta. Para la elaboración de este trabajo han sido fundamentales la utilización de algunas referencias especializadas. Caben destacar las siguientes. En el primer apartado se emplean [1] y [2] para exponer aquellos aspectos más esenciales de las ecuaciones de Navier-Stokes y de la capa límite, así como las hipótesis básicas. Especial mención merece el trabajo fin de grado [3] con autoría de Adrián Carriba Merino, que me inspiró en una buena parte en lo concerniente a la comprensión del comportamiento que tienen las partículas de un flujo en presencia de radiación. Entre otros, [4] y [5] me ayudaron a aclarar los fundamentos de la radiación para poder estructurar de manera eficiente el apartado 2, mientras que [2] y el trabajo fin de grado [6] me sirvieron de referencia para las hipótesis que sustentan convección libre. El modelo teórico sin dispersión con el que se obtiene la expresión general del gradiente del flujo de calor por radiación proviene de [7]. Las referencias con las que se validan los resultados correspondientes a capa límite viscosa thick y thin generada por convección forzada son [8] y [4], mientras que [9] y [10] me ayudaron al análisis del comportamiento de la capa límite térmica según la longitud recorrida en la placa. Para terminar, los resultados concernientes a convección libre se basan en los de [6]. xvi Índice de Figuras Figura 1: distribución de velocidades de un flujo viscoso cerca de la pared, siempre que Re>>1 [2] 5 Figura 2: distribución de velocidades de un flujo viscoso cerca de la pared [2]. Nótese el desprendimiento de la capa límite cuando avanzamos más allá del punto S 8 Figura 3: capa límite por convección natural en una placa vertical [12] 8 Figura 4: aeronave. Fuente: El País Figura 5: reentrada de nave espacial. Fuente: Wikipedia 10 Figura 6: ejemplo de convección libre en naturaleza [6] 10 Figura 7: Representación de los mecanismos físicos que tiene lugar en el seno del fluido en convección libre [6] 11 Figura 8: espectro electromagnético [5] 13 Figura 9: propiedades radiativas de una superficie [15] 14 Figura 10: atenuación de radiación debido a absorción y dispersión [13] 15 Figura 11: reflejo de la luna en el casco de un astronauta que presenta un fino recubrimiento de oro [5] 16 Figura 12: equilibrio radiativo en la Tierra [3] 16 Figura 13: absortividad espectral de una mezcla isoterma de nitrógeno y dióxido de carbono [13] 17 Figura 14: ley de Kirchhoff [13] 18 Figura 15: energía específica del cuerpo negro en función de su temperatura y de la longitud de onda [13] 20 Figura 16: hemisferio mostrando las dimensiones de un diferencial de ángulo sólido [5] 21 Figura 17: intensidad de radiación [4] 22 Figura 18: modelo físico de radiación unidimensional [7] con aproximación al método de las diferencias finitas; se estudia la placa 1 como la condición de contorno “w” y la placa 2 como la “0” 23 Figura 19: atenuación debido a absorción y dispersión a lo largo de ds [4]; el razonamiento espectral es análogo al aquí expresado, en función de magnitudes totales 26 Figura 20: malla de discretización [2] 33 Figura 21: resolución numérica de la capa límite mecánica para una forma genérica; y condición de desprendimiento 35 Figura 22: tratamiento numérico a las componentes de velocidad de la capa límite mecánica 38 Figura 23: separación del flujo entre dos regiones con intercambios energéticos diferentes 43 Figura 24: valores de Ncr para diferentes sustancias según la temperatura [4] 45 Figura 25: sistema de ejes coordenados empleado en [8] 51 Figura 26: comparación con [8] de la distribución de temperaturas en función de η escalada, en la última estación de la placa plana. 𝜃𝑤=0.1; 𝜃0=1; 𝑁𝑐𝑟={0.1,1,10} 53 Figura 27: comparación con [8] de la distribución de temperaturas en función de η escalada, en la última estación de la placa plana. 𝜃𝑤=1; 𝜃0=0.1; 𝑁𝑐𝑟={0.1,1,10} 53 Figura 28: comparación con [8] del gradiente de temperaturas, para un caso de pared caliente y otro de pared fría, en función de η escalada, en la última estación de la placa plana. 𝑁𝑐𝑟={0.1,1,10} 54 Figura 29: comparación con [8] del gradiente de temperaturas en la pared fría (𝜃0=1) en función de la temperatura de la pared, en la última estación de la placa plana. 𝑁𝑐𝑟={0.1,1,10} 55 Figura 30: comparación con [8] del gradiente de temperaturas en la pared caliente (𝜃𝑤=1) en función de la xviii temperatura del infinito, en la última estación de la placa plana. 𝑁𝑐𝑟={0.1,1,10} 55 Figura 31: comparación de la distribución de temperaturas en función de η, para Ncr=0.1 (a, b y c) y para Ncr=1 (d,e y f), para tres estaciones longitudinales hasta ξ=1, Nxi=2000; y para placa plana, cilindro ideal (σ=1) y elipse ideal con σ=0.4 57 Figura 32: distribución de temperatura para flujo no viscoso que no conduce calor y que no dispersa radiación, calculada por el código original atribuido a este trabajo, y calculada por el código basado en [3], cuando ambos códigos han alcazado el régimen estacionario a lo largo de la placa plana 58 Figura 33: salto de temperatura en la pared en flujo no viscoso. Pared caliente (𝜃0=0.1;𝜃0=0.9) y pared fría (𝜃𝑤=0.1;𝜃𝑤=0.9). Comparación con resultados análogos para radiación lineal [4] 59 Figura 34: distribución de temperaturas para flujo no viscoso con Ncr=0 (a, b y c), y con Ncr=1 (d, e y f); para tres posiciones en la placa plana, hasta ξ=1, Nxi=2000, y para τ=20 60 Figura 35: distribución de temperaturas para flujo no viscoso con Ncr=0 (a, b y c), y con Ncr=1 (d); para tres posiciones en la placa plana, hasta ξ=1, Nxi=2000, y para τ=0.2 60 Figura 36: distribución de temperaturas para flujo no viscoso con Ncr=0 (a, b y c), y con Ncr=1 (d); para tres posiciones en la placa plana, hasta ξ=1, Nxi=2000, y para τ=0.002 61 Figura 37: comparación entre segunda componente del vector de temperaturas del flujo ideal designado, y condición de contorno exterior de Howell para el caso específico de hipótesis thin en capas límite viscosas 62 Figura 38: código que implementa la condición de contorno de Howell 62 Figura 39: evolución de 𝑁𝑢𝑅𝑒 con ξ para pared fría. Comparación con [4] 63 Figura 40: evolución que describe el gradiente de temperaturas Θ𝑤′ en función de 𝜃𝑤/𝜃0 para cuatro estaciones diferentes de la placa plana. Comparación con [4]. 64 Figura 41: comparación de la distribución de temperaturas en capa límite viscosa con espesor thin, con Pr=1 y con Nxi=2000 resuelta con códigos originales y con consideraciones expuestas en [4] 64 Figura 42: distribución de temperaturas en función de τ de las tres posibles soluciones originales, relajando hipótesis (𝑓(𝜏)), hipótesis thick e hipótesis thin para Ncr=2, Ntau=1000 y Nxi=2000. En este caso se produce en torno al 0.01% de error entre flujo general y thick en ξ=1, ópticamente cerca de la placa 65 Figura 43: distribución de temperaturas en función del espesor óptico de f(τ)), thick y thin, para Ncr=2⋅10− 6, Ntau=1000 y Nxi=2000. En este caso se produce menos del 0.01% de error entre flujo general y thin en ξ=0.0045023. 66 Figura 44: distribución de temperaturas en función del espesor óptico de f(τ), thick y thin para Ncr=0.1, Ntau=1000 y Nxi=2000; y superposición con la distribución de temperaturas de flujo ideal 66 Figura 45: distribución de temperaturas en función de η de f(τ) para los tres valores típicos de Ncr, Ntau=1000 y Nxi=2000 67 Figura 46: evolución de 𝑁𝑢𝑅𝑒 con ξ en el caso de pared fría y pared caliente, y los tres valores típicos de Ncr 68 Figura 47: gradiente de temperaturas frente a 𝜂𝑣 calculado a partir de f(τ) en las estaciones indicadas en la leyenda, para los valores típicos de Ncr 69 Figura 48: perfiles de temperatura adimensional θ, y perfiles de velocidad u y v, para tres estaciones diferentes de la placa plana, Ncr={0.1,1,10,∞} y Pr=1 72 Figura 49: perfiles de temperatura adimensional θ, y perfiles de velocidad u y v, para tres estaciones diferentes de la placa plana, Ncr={0.1,1,10,∞} y Pr=100 72 Figura 50: perfiles de Θ respecto a 𝜂𝑒𝑠𝑐 para Pr=1 y Pr=100, Ncr={0.1,1,10,∞}, y en la ultima estación de la placa plana, es decir, ξ=1 73 Figura 51: perfiles de 𝑢𝑒𝑠𝑐 respecto a 𝜂𝑒𝑠𝑐 para Pr=1 y Pr=100, Ncr={0.1,1,10,∞}, y en la ultima estación de la placa plana, es decir, ξ=1 74 Figura 52: evolución de 𝑁𝑢/𝐺𝑟1/4 respecto a ξ, para Pr=1 y Pr=100, y Ncr={0.1,1,10,∞} 75 xx 1 1 Capa límite: convección forzada y convección libre 1.1 Introducción a las capas límite. Tratamiento de las ecuaciones de Navier Stokes En este proyecto se trabaja con flujo ideal, o también denominado flujo no viscoso, y con flujo real, o viscoso. Interaccionan, en general, con una o dos placas planas que están orientadas horizontal o verticalmente. En todos los casos, el concepto de capa límite salta a la palestra. Cuando hablamos de capa límite nos referimos indistintamente al comportamiento, en términos de velocidad y de temperatura, que presenta un flujo muy cerca (más adelante se discute cuánto de cerca) de una superficie, o a esa zona cercana a la pared donde el flujo exhibe dicho comportamiento. Para que se genere una capa límite mecánica y térmica es condición necesaria, entre otros, que el fluido en estudio tenga viscosidad. Si el fluido en estudio no es viscoso, pero transmite calor por conducción, solo se genera una capa límite térmica. Por otro lado, es importante conocer si el campo de velocidades y la densidad del fluido dependen de su campo de temperaturas o no. En este trabajo se va a estudiar radiación en convección forzada, cuando no existe dicha dependencia; y radiación en convección libre, cuando sí existe dicha dependencia. Si en la convección forzada también existiese, se hablaría de convección mixta. Así pues, en este primer apartado, se realiza una somera descripción de la naturaleza física de los fluidos ideales y no ideales; se explica aquellos términos que se van a conservar de las ecuaciones de Navier-Stokes; se enfatiza la necesidad de trabajar con números adimensionales en la mecánica de fluidos, en presencia o no de radiación térmica; y se explican las hipótesis que nos permiten sostener el modelo de capa límite mecánica y térmica en convección forzada y en convección natural que se va a resolver a través de los métodos numéricos del apartado 3. En el apartado 2, se explican los fundamentos de la radiación térmica que son relevantes en este proyecto. Finalmente, en los apartados 4 y 5 se muestran los resultados buscados, validados tras contrastarse con referencias de interés. En los fluidos no viscosos, no existen fuerzas tangenciales entre capas adyacentes de fluido, ni entre el fluido y las condiciones de contorno, solo existen fuerzas normales (presión). Dicho de otro modo, [2] el flujo ideal no se opone a un cambio en su forma a través de cualquier resistencia interna. De esta manera, si bien esta teoría ha tenido una gran evolución desde el punto de vista matemático, y ofrece una descripción satisfactoria en el movimiento de ondas en un fluido, o en la formación de chorros líquidos y capilares, no es útil para calcular la resistencia que experimenta un cuerpo que se mueve de manera relativa respecto a un fluido. En este caso la hipótesis de flujo ideal, al no tener en cuenta las fuerzas tangenciales antes mencionadas, predice la paradoja de D’Alembert, por la cual, un cuerpo que se mueve de manera subsónica y uniforme a través de un fluido ideal infinito no experimenta resistencia al avance, cuando la experiencia dicta que sí la experimenta. Si bien en este trabajo, no se pretende calcular dicha resistencia al avance, sí se pretende estudiar las propiedades del fenómeno físico que la explica cuando se consideran efectos de radiación térmica. Del mismo modo, para conocer la razón de ser de dicho fenómeno, se debe considerar una propiedad física de los fluidos, la viscosidad. Históricamente, la viscosidad de los fluidos se había despreciado de manera generalizada, ya que fluidos tan usados como el agua o el aire, tienen una viscosidad muy baja. En muchos casos, estos fluidos de baja viscosidad se comportan como los flujos ideales, porque las fuerzas tangenciales en las capas del fluido son muy pequeñas. Sin embargo, en el flujo ideal hay una diferencia entre las velocidades tangenciales que presentan las moléculas de fluido que se encuentran en las proximidades de una condición de contorno y en la propia condición de contorno. En teoría, nada impide que el fluido deslice por la condición de contorno. En cambio, el fluido real no desliza, porque la viscosidad, por muy pequeña que sea, le obliga a adherirse a la condición de contorno, debido a las fuerzas tangenciales que aparecen entre ésta y el fluido. La presencia de dichas fuerzas tangenciales y la condición de no deslizamiento constituyen la diferencia esencial entre fluidos viscosos y no viscosos. Capa límite: convección forzada y convección libre 2 2 La mecánica de fluidos es una disciplina que presenta una fuerte dependencia con la experimentación para validar las hipótesis susceptibles de ser planteadas para modelar y comprender un fenómeno. Al mismo tiempo, debe apoyarse en todo momento en las ecuaciones simplificadas de Navier-Stokes. Se trata de un conjunto de ecuaciones en derivadas parciales no lineales que describen el movimiento general de un fluido respecto a un objeto, entendido este último como parte de las condiciones de contorno junto con la condición de no deslizamiento, la velocidad lejos del objeto, y la temperatura en él y lejos de él. Deducidas por Claude-Louis Naviere y por George Gabriel Stokes, aún no se dispone de una solución general, debido a su colosal complejidad. La manera que tiene la mecánica de fluidos de tratar con ellas para resolver problemas es a través de números adimensionales definidos convenientemente, a través de hipótesis y a través de métodos numéricos. El sistema de ecuaciones de Navier Stokes [1] proporciona los campos de densidad (ρ), velocidad (𝑣), presión (p) y temperatura (T) en cualquier punto del dominio fluido, incluida la superficie del objeto con el que el fluido interactúa. Se compone de seis ecuaciones y seis incónitas. Las ecuaciones son la conservación de la masa (1 ecuación), la conservación de la cantidad de movimiento (3 ecuaciones), la conservación de la energía (1 ecuación) y la ecuación de estado para los gases perfectos (1 ecuación). Las incógnitas son ρ, 𝑣 (tres componentes, u, v y w), p y T. Se muestran las ecuaciones en su forma más general para encarar los problemas que se resuelven en este trabajo: ∇⋅(ρv 󰇍  )=0 (1.1) ρv 󰇍  ⋅∇v 󰇍  =−∇p+ρfm 󰇍 󰇍 󰇍 󰇍  +∇⋅τ´ (1.2) ρCpv 󰇍  ⋅∇T=∇⋅(k∇T)−∇⋅qr (1.3) pρ=RgT (1.4) Se trata de las ecuaciones para un movimiento estacionario en el tiempo. La densidad es constante, a expensas de matizar esta información en convección natural. Sabiendo además que, desde un punto de vista matemático, la divergencia de un campo vectorial como es 𝑣 mide la diferencia entre el flujo saliente y entrante de dicho campo vectorial sobre la superficie que rodea a un volumen de control, entonces (1.1) expresa un movimiento incompresible. En el primer miembro de (1.2) tenemos el término convectivo de la derivada total, que mide la variación de la cantidad de movimiento de una partícula fluida a lo largo de su movimiento en el volumen de control en estudio. Dicho término constituye una fuerza que se iguala con los posibles gradientes de presión, fuerzas másicas y difusión del momento. Con 𝐶𝑝 el calor específico a presión constante, y “k” la conductividad térmica, en (1.3) se tiene el balance de energía que tiene lugar en el volumen de control, el cual, en nuestro caso será el flujo confinado entre las dos placas horizontales, o el flujo adyacente a una placa vertical. Según el primer principio de la termodinámica, el incremento de energía total que se produce en el volumen de control es igual al calor suministrado a cada elemento de masa. El término de la izquierda mide la variación de la energía interna de una particula fluida a lo largo de su movimiento en el volumen de control en estudio. Dado que el transporte de calor por convección se debe al movimiento del fluido, esté o no esté obligado a circular a cierta velocidad por la superficie en estudio; y dado que el término de la derecha se corresponde con el calor suministrado al volumen de control por radiación y por conducción, el balance de energía cobra sentido físico: el calor que se propaga por convección se iguala al calor por conducción y por radiación. Estos dos últimos medios de propagación del calor, como se verá, gobiernan los problemas a resolver en este trabajo. Como se verá en los apartados 4 y 5, dicho balance puede modificarse durante el avance del fluido. Por otro lado, la ecuación de estado se empleará cuando haya que justificar ciertas hipótesis en el marco de la convección natural. El fluido de trabajo es en todo momento un gas perfecto. En el apartado 2 se justifica su elección. 3 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 1.2 Teorema Pi. Parámetros adimensionales empleados Como se ha comentado previamente, uno de los recursos más versátiles para exprimir el potencial de las ecuaciones de Navier-Stokes es el análisis dimensional. Esta herramienta [11] se basa en el principio de que las relaciones funcionales entre las variables que intervienen en un fenómeno físico dado son independientes del sistema de unidades elegido. Dicho principio permite expresar toda la información acerca de las variables físicas contenida en las relaciones funcionales en términos de un número menor de variables o parámetros adimensionales obtenidos a partir de los originales. En otras referencias, se refieren a dicho principio como leyes de similaridad, según las cuales, [2] las variables a calcular pueden referenciarse respecto a valores de referencia, obteniéndose de la misma manera parámetros adimensionales. Este principio queda plasmado en el Teorema Pi o de Vaschy-Buckingham. Considerese una relación funcional [11] entre las cantidades físicas 𝑎0,𝑎1,...𝑎𝑛. Si se tienen “k” dimensiones físicas independientes, como pueden ser longitud, tiempo, masa y temperatura, este teorema establece que pueden definirse “n+1-k” variables o parámetros adimensionales. A continuación, se hace una relación de aquellos que son de nuestro interés, tanto para convección forzada como para convección libre, en presencia de radiación, y que irán apareciendo a lo largo del trabajo conforme se vaya desarrollando el proyecto. • Número de Eckert: mide la importancia relativa de la convección de energía cinética frente al efecto de convección de calor. Se ha supuesto mucho menor que la unidad tanto en convección forzada como en convección natural [1] ya que se desprecia la disipación de energía cinética en energía interna. Ec=u∞ 2 Cp∆T (1.5) • Número de Froude: mide la importancia relativa de la convección de energía cinética frente al término de fuerzas másicas. Fr=U∞ 2 Lg (1.6) En convección forzada es mucho menor que la unidad, mientras que en convección libre es de orden unidad. • Número de Reynolds: se define de la siguiente forma. Re=ρU∞L μ=U∞L ν (1.7) Siendo µ el coeficiente de viscosidad del fluido, ν su viscosidad cinemática y L la longitud característica del objeto en estudio, generalmente una placa plana (del orden de metros). Mide la importancia relativa entre las fuerzas de inercia convectiva, y las fuerzas de viscosidad. En este trabajo se va a trabajar con régimen laminar, y “Re” toma un valor de 𝑅𝑒=104. • Número de Prandtl: se define de la siguiente forma. Pr=να=μ⋅Cp k (1.8) Siendo α la difusividad térmica. Se empleará profusamente en el apartado 3 la primera definición. Se trata de una propiedad del fluido, al igual que la viscosidad, y compara la capacidad de un fluido para transferir, por difusión, cantidad de movimiento y calor. Toma un valor de Pr=1, salvo en convección libre, que se le dan dos valores, Pr=1 y Pr=100. Capa límite: convección forzada y convección libre 4 4 • Número de Peclet: se define de la siguiente forma. Pe=Re⋅Pr=U∞L α (1.9) Mide la importancia relativa de la convección de calor frente al calor recibido por conducción • Número de Boltzmann: se define de la siguiente forma. Bo=ρCpU∞ σT∗3 (1.10) Establece una medida cualitativa del orden de magnitud de la convección frente al flujo de calor por radiación. Así mismo, el cambio de variable que se le aplica a la variable longitudinal que se emplea en el análisis de las capas límite en función del espesor óptico puede expresarse, cuando la convección de calor entra en juego, en función de Bo: x=ρCpU∞ 2βσT∗3⋅ξ=Bo 2β⋅ξ (1.11) • Parámetro Ncr: se define de la siguiente forma. Ncr=κk 4σT∗3 (1.12) Siendo κ la absortividad del medio, con unidades recíprocas, este número adimensional aparece de manera natural en el apartado 3. Establece una relación entre el término conductivo de la capa límite y el término radiativo. De hecho, su valor puede dar pistas sobre si el número de Pr puede tomarse grande o pequeño. Si Ncr<<1, entonces Pr podría tomarse Pr>>1, momento en el que la capa límite térmica tendría un espesor mucho menor que el de la capa límite mecánica. Si Ncr>>1, entonces Pr podría tomarse Pr<<1, momento en el que la capa límite térmica tendría un espesor mucho mayor que el de la capa límite mecánica. En la página 213 de [2] puede hallarse más información sobre estos supuestos. • Espesor óptico: en el apartado 2.4 se aportará más información. De momento se tiene: τ=∫κ⋅dy y 0 (1.13) Siendo “y”, la coordenada vertical que se maneja en este trabajo. El espesor óptico y Ncr son los parámetros que modificaremos con mayor frecuencia con el objeto de conocer cómo se comporta un medio, en este caso flujo ideal y capa límite, frente a radiación térmica. • Número de Nusselt: se trata de un número adimensional que mide el aumento de la transmisión de calor por convección comparada con la transferencia de calor si ésta ocurriera solamente por conducción. Puede existir interés por un Nusselt local, correspondiente a un punto concreto de la superficie en estudio, o por un Nusselt medio, que puede obtenerse integrando el Nusselt local en toda la superficie. En este trabajo, se va a estudiar por defecto el número de Nusselt local, que se toma como el gradiente de temperaturas en la superficie, en tres formas diferentes: gradiente de temperaturas persé, 𝑁𝑢/√𝑅𝑒 (convección forzada) y Nu/(Gr)1/4 (convección libre). • Número de Rayleigh: mide la importancia relativa entre los efectos de las fuerzas de flotabilidad y los efectos de las fuerzas de viscosidad de la conducción térmica. • Número de Gr: mide la importancia relativa entre los efectos de las fuerzas de flotabilidad y los efectos de la viscosidad del fluido. Ambos parámetros se usan en el estudio de los fenómenos de convección libre: 5 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Ra=gβΔTL3 να ; Gr=gβΔT⋅L3 ν2 (1.14) Valga el siguiente comentario para contextualizar la relevancia del análisis dimensional. El choque de partículas entre sí en un flujo sea ideal o viscoso, permite la transmisión de calor por conducción o difusión. Sin embargo, la absorción de las partículas de la energía de los fotones no tiene porqué venir acompañada de estos choques. Por un lado, el parámetro Ncr da una idea de la existencia o no de estos choques de partículas, mientras que el espesor óptico da una idea de los obstáculos que los fotones encuentran en su recorrido a la velocidad de la luz hasta que finalmente son absorbidos por las partículas del medio. Efectivamente, a la hora de validar en la realidad los resultados que se muestran en este proyecto, los experimentos a realizar deberían considerar entre otros, los valores de la absortividad, la conductividad térmica, la velocidad del fluido o la temperatura de referencia, de los que dependen τ y Ncr, los parámetros que gobiernan los problemas en estudio de convección forzada. 1.3 Capas límite generadas por convección forzada. Hipótesis desplegadas En este apartado se alumbran las implicaciones físicas que tiene los órdenes de magnitud que toman algunos de los números adimensionales referidos anteriormente, en el caso de convección forzada. En el siguiente apartado se hace lo mismo, en el caso de convección libre. Cuando se obliga a un flujo viscoso a circular por una superficie horizontal a una velocidad de referencia constante 𝑈∞, [1] la cual, para las aplicaciones en las que se genera capa límite, toma valores del orden de 100m/s, el término de fuerzas másicas es varios órdenes menor que el de la aceleración convectiva. En consecuencia, Fr>>1, y se desprecia, de momento, la aceleración de la gravedad. Por otro lado, resultan de nuestro interés los fluidos de baja viscosidad. Recordando (1.7), se tiene que Re>>1, y, por tanto, Pe>>1, ya que Pr=1. Es decir, siempre y cuando nos situemos a una distancia L (contemplada en Re) lo suficientemente lejos del objeto que interacciona con el flujo, se pueden despreciar los términos viscosos y los de conducción de calor, siendo válida la hipótesis de flujo ideal cometiéndose errores relativos del orden de 𝑅𝑒−1 y de (RePr)−1. Sin embargo, estos términos han de ser retenidos en una región del espacio cercana a la superficie del objeto, ya que representan en las ecuaciones aquellos mecanismos físicos responsables de forzar tanto la condición de contorno de igualdad de velocidades tangenciales entre el fluido y el sólido, como la igualdad de temperaturas y de flujo de calor entre el fluido y el sólido: Figura 1: distribución de velocidades de un flujo viscoso cerca de la pared, siempre que Re>>1 [2] Ese espacio cercano a la superficie del sólido se conoce como capa límite. Fue identificado por Ludwig Prandtl en 1904 en su trabajo “On the motion of fluids in very little friction”. Como se ve en la Figura 1, [2] a medida que el flujo recorre una distancia mayor, la capa de partículas que es frenada por la fricción se hace cada vez más grande. Si bien la transición entre capa límite y flujo exterior se produce de manera contínua en régimen laminar, se introduce de manera artificial “𝛿” para designar al espesor de la capa límite. Capa límite: convección forzada y convección libre 12 12 −∂p ∂x+ρfm 󰇍 󰇍 󰇍 󰇍  =−(ρ−ρ∞)g 𝐢 (1.35) Con esto en mente, y [6] siempre que las variaciones relativas de densidad en el proceso sean pequeñas: • La densidad se considera constante e igual a la de referencia, en este caso, 𝜌∞, en todos los términos de las ecuaciones de Navier-Stokes excepto en el de la variación de la presión en función de la variable longitudinal de la placa, en el cual, se aproxima por: ρ=ρ∞[1−β(T−T∞))] (1.36) Siendo β el coeficiente de expansión térmica. • Todas las propiedades del fluido (ν, α y β) se consideran constantes y tomadas a la temperatura y densidad de referencia. • Los términos de disipación viscosa son despreciables. Es decir, como ya se hace en convección forzada, Ec<<1 (ver apartado 2.5). Por tanto, (1.33) queda: ρ(u𝑥∂u𝑥 ∂x+v𝑦∂u𝑥 ∂y)=gβ(T−T∞)+μ∂2u𝑥 ∂y2 (1.37) Como puede verse, la ecuación de la cantidad de movimiento depende del perfil de temperaturas, ya que, como se comenta en el primer apartado del trabajo, la convección libre se modela con una densidad dependiente de la temperatura. En el apartado 3.4 se le añade al sistema (1.32)-(1.34) el término de radiación pertinente, y se le aplican los métodos numéricos que se explican en los apartados precedentes. Para cerrar el sistema, las condiciones de contorno quedarían: ux(x,0)=vy(x,0)=ux(x,y→∞)=0;T(0,y)=T∞ ,T(x,y→∞)=T∞,T(x,0)=Tw (1.38) Aquí concluye esta introducción de aquellos aspectos de la mecánica de fluidos relevantes en este proyecto. Pasamos a explicar los fundamentos de la radiación térmica. 13 2 Introducción a la radiación térmica 2.1 La radiación como onda-corpúsculo. Aplicaciones en ciencia e ingeniería Cualquier tipo de materia cuya temperatura sea mayor que el cero absoluto emite energía electromagnética, que se propaga en todas las direcciones a la velocidad de la luz. Dicha energía, según la dualidad onda-corpúsculo, puede tener, o bien, naturaleza de onda electromagnética, o bien, naturaleza de fotones. Los dos puntos de vista se complementan. Por un lado, [3] [4] la teoría de las ondas electromagnéticas suele usarse para predecir las propiedades radiativas de ciertas superficies, mientras que la teoría de los cuantos de energía se emplea para determinar las propiedades de los medios que absorben y emiten. En este trabajo no se asume ninguno de los dos puntos de vista en especial, aunque a la hora de comentar los resultados desde el punto de vista físico, se hará uso del concepto de fotón, en vista de que se va a proceder al estudio de los efectos de la radiación en capas límite que, en mayor o menor medida, pueden absorber y emitirla. Dicho esto, la energía electromagnética viene dada por tres parámetros: la frecuencia ν, la longitud de onda λ, y la velocidad de la onda o de los fotones c, siendo independientes dos de ellos, a través de la siguiente relación: c=λν (2.1) La velocidad c y la longitud de onda en un cierto medio se relacionan con sus correspondientes valores en el vacío, 𝑐0 𝑦 𝜆0=𝑐0/𝜈, mediante: c=c0 n λ=λ0 n (2.2) donde n es el índice de refracción (en este trabajo se asume de valor unidad). De esta manera, [5] el interés del proyecto reside en la energía electromagnética con un rango de longitud de onda de entre 0.1 𝜇𝑚 𝑦 100 𝜇𝑚 (parte del espectro ultravioleta, todo el espectro visible y parte del espectro infrarrojo). Figura 8: espectro electromagnético [5] Dicho rango se corresponde con energía, denominada de aquí en adelante radiación térmica, emitida por un medio debido únicamente a su temperatura. Es una de las tres maneras en las que el calor puede transmitirse, junto con la convección y la conducción. En este trabajo se estudia cómo se relacionan estos tres modos de transmisión de calor. Como característica diferenciadora, la radiación térmica puede propagarse por el vacío. La radiación puede, en general, [3] reflejarse, transmitirse, ser emitida por un cuerpo o partícula (emitting), ser Introducción a la radiación térmica 14 14 absorbida por dicho cuerpo (absortion), a la vez que emite radiación propia, dispersarse (scattering) o atenuarse. La relación existente entre los niveles cuánticos de energía de átomos y moléculas, de origen electrónico, vibracional y rotacional gobierna las leyes de la emisión de energía, que se propaga por un medio entre el objeto emisor y el objeto que absorbe. Por otro lado, el efecto de absorción, sumado al de dipersión, se denomina atenuación (extinction). Si los fotones impactan con moléculas del medio en estudio ( [4] o incluso con inhomogeneidades), [14] éstas se excitan a un nivel de energía mayor y casi inmediatamente (en la escala de nanosegundos), vuelven a su estado energético previo y emiten fotones de igual o diferente frecuencia que aquellos que absorbieron en primer lugar. En este caso se ha producido dispersión de la radiación. Si en cambio, las moléculas excitadas colisionan con otras antes de que se produzca la reemisión, el salto de energía se traduce en un aumento de energía cinética de las moléculas que reciben el impacto de los otras. En este caso se produciría absorción de los fotones, y, en función del valor de Ncr, transmisión de calor por conducción. En este proyecto se supone que, si bien las moléculas se ven excitadas por el impacto de los fotones, no existe dispersión, por lo que el nivel de atenuación de la radiación coincide con la absorción producida. Si la radiación se atenúa completamente en su camino por el medio, éste es opaco. Si lo atraviesa sin atenuarse, la onda se transmite completamente y el medio es transparente. Si se atenúa parcialmente, entonces el medio es semitransparente. Del mismo modo, [7] la radiación puede interactuar con la convección y la conducción de dos maneras diferentes, ya comentadas anteriormente. Por un lado, puede darse fenómenos de radiación a través de las condiciones de contorno de los procesos que involucran conducción y convección. En este trabajo dicha condición de contorno será generalmente una placa plana (si bien también se harán algunos comentarios para otras formas). Así, es posible que la temperatura de dicha placa se vea modificada debido a la interacción de radiación térmica con los otros dos modos de transmisión de calor, si bien, en este trabajo se hace la hipótesis de placa isoterma. Dicho esto, según las propiedades radiativas de la superficie que contacte la radiación (más adelante se definen), una fracción de ésta se refleja, otra se convierte en energía interna a través de la absorción de energía, y otra se transmite a través de la placa (en nuestro caso, trabajaremos con placas opacas). Figura 9: propiedades radiativas de una superficie [15] Por otro lado, puede darse fenómenos de radiación volumétricos; es decir, la energía atraviesa un medio que puede absorberla o emitirla. Debido al hecho de que la expresión analítica de la potencia emisiva depende de la potencia a la cuarta de la temperatura, problemas de este tipo presentan un claro comportamiento no lineal. Así mismo, la radiación que finalmente absorbe un elemento, tras ser emitida por otros elementos situados en la vecindad o lejanos, da lugar a una expresión integral. Teniendo en cuenta además el tratamiento numérico de dicha expresión integral, así como el empleo del método de las diferencias finitas, que nos aportan expresiones algebraicas; así como la dependencia de las magnitudes, a priori, de las 3 dimensiones del espacio, o la frecuencia; y la necesidad de resolver la capa límite para introducirla, como dato conocido en cada estación, en la ecuación de conservación de la energía, consecuentemente el sistema de ecuaciones de Navier-Stokes que se resuelve de manera recursiva en este trabajo es, en los sucesivos pasos que se dan en su resolución, no linealintegro-algebraico-diferencial parcial acoplado. Véase la siguiente ilustración para visualizar cómo la radiación en forma de fotones interacciona con un elemento volumétrico según lo ya comentado. En base a ello se definirá más adelante la ecuación de transferencia de radiación (ETR). 15 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Figura 10: atenuación de radiación debido a absorción y dispersión [13] De esta manera, la radiación térmica es uno de los factores más importantes que tener en cuenta en numerosos ámbitos de la ciencia y la tecnología. Destaca especialmente aquellos que involucran altas temperaturas, [3] habida cuenta de la dependencia con la potencia a la cuarta de la temperatura que presenta la energía por radiación térmica emitida por un cuerpo a dicha temperatura, que convierte a la radiación en un eficiente mecanismo de transmisión de calor en hornos, cámaras de combustión, intercambiadores de calor o en la reentrada de una nave en la Tierra en régimen hipersónico. En este último caso, al friccionar el aire con la nave, la temperatura exterior a ésta se eleva considerablemente, y el material sufre una absorción considerable de energía térmica, para lo cual el correspondiente escudo térmico debe estar correctamente diseñado. La radiación térmica es tan importante que, en vista de que es la única manera en que el calor puede propagarse por el vacío, puede explicar fenómenos de la vida real, como, por ejemplo, [5] el enfriamiento que ocurre durante una noche con el cielo despejado, debido que a que el calor que irradian las personas o el acumulado por el suelo durante el día, se pierde al exterior. Es también la razón de que pueda existir congelación por encima de los 0º centígrados, por ejemplo, en las plantas de una floristería, cubiertas con plástico, cuando no existe ningún emisor de radiación térmica en los alrededores de las plantas, ni siquiera en la floristería; de que, en las misiones espaciales, las naves se abastezcan de energía solar, para aprovecharla a bordo, por ejemplo, haciendo funcionar sistemas sin tener que incrementar el peso a través de baterías o generadores; o incluso de que, en las misiones Apolo, la nave tuviera que reorientarse cada cierto tiempo para cambiar las partes de la nave expuestas a la radiación solar, y evitar su sobrecalentamiento. Por otro lado, como se verá a lo largo del trabajo, los efectos de la radiación antes mencionados dependen, no solo de su frecuencia, sino también del medio o medios con los que interacciona. En función del espesor óptico de estos, la radiación penetrará más o menos sin ser absorbida. Al mismo tiempo, en el caso de que la radiación esté compuesta por más de una longitud de onda, como es el caso de la luz solar, [5] puede estudiarse un recubrimiento de cierto espesor, que absorba tan solo ciertas frecuencias. Tal es el caso del oro. Como material presenta gran absorción, pero para ciertos espesores pequeños, se comporta como un recubrimiento semitransparente que deja sólo deja pasar la radiación visible. En consecuencia, es ampliamente utilizado para proteger a los astronautas de longitudes de onda dañinas, a modo de filtro en sus cascos. En nuestro trabajo, asumimos radiación monocromática. Introducción a la radiación térmica 16 16 Figura 11: reflejo de la luna en el casco de un astronauta que presenta un fino recubrimiento de oro [5] El equilibrio radiativo que ocurre en la atmósfera ilustra aquello que se comentó previamente que la radiación térmica podía hacer. De esta manera, las capas más bajas de la atmósfera se encargan de absorber el espectro de la luz solar más perjudicial para el ser humano, [3] atenuándose la energía de la radiación que llega a la superficie terrestre, la cual absorbe y refleja un cierto porcentaje de dicha radiación. A su vez, un cierto porcentaje de la radiación se dispersa en otras capas de la atmósfera, como por ejemplo ocurre en la formación de un arco iris, cuando sale el Sol mientras llueve; o en los atardeceres, cuando [4] se dispersan especialmente aquellas longitudes de onda menores, correspondientes al color azul, y permanecen aquellas tonalidades más rojizas. A este equilibrio también se suma la energía que emite la Tierra, de la cual, una parte no puede escapar al espacio porque es reflejada por la atmósfera, debido al efecto invernadero; otra parte es absorbida por la atmósfera, y finalmente, aquella radiación de longitud de onda grande sí se transmite al espacio exterior. Si bien, las propiedades radiativas de superficies y/o elementos de volumen pueden ser dependientes de la longitud de onda, en este trabajo se consideran independientes (se toma la hipótesis de medio gris). Figura 12: equilibrio radiativo en la Tierra [3] Para finalizar, se hace un breve comentario sobre el medio con el que se va a trabajar en este proyecto. Como se ha mencionado en el apartado anterior, se trata de las capas límite que se generan en una placa plana horizontal; y las que se generan por convección libre en una placa plana vertical. En ambos casos, se pone en juego un cierto flujo de un gas, típicamente aire. Dicha asunción no es baladí en un trabajo que versa sobre interacción de la radiación con un medio. No en vano, las hipótesis de medio gris e índice de refracción unidad, si bien son asumibles para un gas como el aire, restan generalidad a este proyecto, al menos desde el punto de vista de la radiación. Como muestra un botón. 17 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Figura 13: absortividad espectral de una mezcla isoterma de nitrógeno y dióxido de carbono [13] Como puede verse, dependiendo de la temperatura a la que ocurra el proceso, y dependiendo de la longitud de onda, la propiedad de la absortividad, así como el resto de las propiedades, que se mencionarán a continuación, ofrecen variaciones importantes. Aun así, cabe resaltar que este proyecto sí ofrece generalidad, por ejemplo, en el tratamiento numérico a la hora de resolver las capas límite, con la posibilidad de ser resueltas para flujos alrededor de otras formas que no sean placas planas. Una vez expuestas unas líneas maestras del potencial que tienen los estudios sobre radiación, y algunas de las principales hipótesis que se despliegan en este trabajo, se pasa a una fundamentación teórica para poner en pie el estudio que se hará en los próximos apartados. 2.2 Conceptos generales. Propiedades radiativas, cuerpo negro, y potencia emisiva Llegados a este punto, puede apreciarse que la radiación térmica, a diferencia de la conducción y la convección, requiere, para cada punto en el espacio, [13] un balance de energía que tenga en cuenta todos los fotones que puedan llegar a ese punto, sin importar lo lejos que estuviesen. Por ello, un balance de conservación de energía debe realizarse en un espacio cerrado por paredes lo suficientemente gruesas que impidan que las atraviese la radiación (es decir opacas). Nada impide, en cambio, que dichas “superficies” puedan absorber o emitir fotones. Estrictamente hablando, una superficie hace referencia a un objeto geométrico en dos dimensiones, en el cual no es posible la atenuación de la radiación. Sin embargo, en la práctica, el espesor de la placa es muy pequeño comparado con las dimensiones totales del espacio cerrado en estudio (del orden de Amstrongs en metales y de micrómetros en no metales). Por tanto, es adecuado hablar de superficie para referirnos a dichas condiciones de contorno. Centrando la atención en una superficie tal, y recordando la Figura 9, podemos tener las siguientes propiedades radiativas para una superficie: Reflectancia ρ= radiación reflejada radiación total que llega a la placa (2.3) Absortividad α= radiación absorbida radiación total que llega a la placa (2.4) Transmisibilidad τ= radiación transmitida radiación total que llega a la placa (2.5) Sabiendo que la radiación se comportará de alguna de estas tres maneras al llegar a la superficie: ρ+α+τ=1 (2.6) Si el medio es opaco, entonces 𝜏=0 y por tanto 𝜌+𝛼=1. Por otro lado, dicho medio puede emitir también radiación, pudiéndose definir: Introducción a la radiación térmica 18 18 Emisividad ε= energía emitida desde una superficie energía emitida por un cuerpo negro a la misma temperatura (2.7) De esta manera, una superficie que, además de ser opaca, no refleja nada de radiación, 𝜌=0, se trata de un cuerpo negro. Ningún otro cuerpo puede absorber tanta energía como éste. De hecho, también es el cuerpo que emite más energía. Es una de las conclusiones de la ley de Kirchhoff. La ley de Kirchhoff establece que, para un sistema en equilibrio termodinámico, como el de la Figura 14, se aplica, para cada una de las superficies del sistema, que la emisividad es igual a la absortividad. Si bien, como puede apreciarse, esta ley se dedujo bajo la condición de que una determinada superficie se encuentra en equilibrio termodinámico con su entorno, es, de hecho, [4] una relación general que se aplica para condiciones de no equilibrio. La razón es el hecho de que tanto la emisividad como la absortividad solo dependen de la naturaleza de la superficie y de su temperatura, mientras el índice de refracción sea la unidad. Como consecuencias de la ley de Kirchhoff, se tienen: • Para un cuerpo negro, que no refleja ni deja pasar la radiación, se tiene que, según la ecuación (2.6), 𝛼=1. Consecuentemente, según la ley de Kirchhoff, 𝜀=1. Por tanto, ningún cuerpo puede absorber ni emitir tanta energía como lo hace el cuerpo negro. • La radiación del cuerpo negro es isótropa o difusa. Independientemente de la dirección, la intensidad es la misma. En cualquier caso, a instancias prácticas, a veces se asume dicha isotropía de la radiación en general, aunque no se trate la del cuerpo negro. En este trabajo se asume igualmente. Figura 14: ley de Kirchhoff [13] En cuanto a las partículas que componen el medio, debe asumirse que se encuentran en equilibrio termodinámico local. Es decir, [4] el total de átomos y moléculas que contribuye a la emisión y absorción de fotones se encuentra en equilibrio. Esta asunción es válida cuando las transiciones de energía están controladas por colisiones moleculares en vez de por absorción y emisión de energía, como es el caso. A su vez, dicha asunción implica [16] que cada punto geométrico que describe el medio está caracterizado por una temperatura T definida. Por tanto, cada elemento de materia en cada punto del espacio se comporta como si formase parte de un entorno cerrado isotermo a dicha temperatura T. Por ello, la emisión es independiente de la radiación incidente en el elemento de volumen, y la ley de Kirchhoff, al igual que en superficies, también se cumple en elementos de volumen. Ello será de especial utilidad para establecer la ETR. Así pues, todos los medios en general sean superficies o elementos de volumen, emiten radiación electromagnética al azar en todas las direcciones, dependiendo de la temperatura local, de la frecuencia y de las propiedades del medio. Dicha radiación puede estudiarse a nivel de intensidad o a nivel de energía. A nivel conceptual, como veremos, la energía de los fotones se integra, bien en un hemisferio o en una esfera unidad, según si el elemento en cuestión es de una superficie o del fluido, para obtener el flujo de calor [W/𝑚2] con el que se trabaja en este proyecto. En un principio los fotones pueden tener cualquier frecuencia y depender de la dirección. Más adelante, en coherencia con la hipótesis de cuerpo gris e intensidad difusa, se trabaja con magnitudes totales y se obtiene una relación entre la intensidad y la energía. Más adelante en este apartado se demostrará que, en convección forzada, la transferencia de radiación en dirección del avance de la capa límite es despreciable bajo ciertas condiciones, por lo que nos centramos en la transferencia de radiación verticalmente. 19 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre Comenzamos por la energía. [13] El flujo de calor por radiación emitido desde una superficie también se llama potencia emisiva. Debe distinguirse entre potencia emisiva espectral 𝐸𝜈 (en una cierta frecuencia) y potencia emisiva total E (integrada en todo el espectro): • Eν (Energía emitida por unidad de tiempo, unidad de superficie y unidad de frecuencia). • 𝐸 (Energía emitida por unidad de tiempo y unidad de superficie). E=∫ Eνdν ∞ 0 (2.8) Eν1−ν2=∫ Eνdν ν2 ν1 (2.9) A la radiación térmica por unidad de tiempo, unidad de superficie y de una cierta frecuencia o longitud de onda se le denomina también radiación monocromática. Como puede verse, 𝐸𝜈𝑑𝜈 se trata de la energía emitida por unidad de tiempo y unidad de superficie correspondiente a las ondas situadas en ese diferencial de frecuencia. Dicha expresión nos permite en la práctica integrar en todo el espectro o entre dos frecuencias, en función de si la potencia que deseamos obtener realmente es en todo el espectro o entre dos frecuencias. Esta reflexión puede realizarse también equivalentemente en función de la longitud de onda. Como se acaba de referir, el cuerpo negro presenta importantes propiedades, que para este trabajo son de especial relevancia puesto que las placas planas que se consideran son cuerpos negros. En la práctica, también aporta una idea de que cómo 𝐸𝜈 varía con 𝜈, puesto que las tendencias que a continuación van a explicarse para la energía que emite un cuerpo negro se corrigen con la emisividad ya definida para obtener una emisión de energía real. De esta manera, aplicando la Segunda Ley de la Termodinámica, así como argumentos de mecánica cuántica, [4] Max Planck demostró que, a una cierta temperatura y longitud de onda, existe una máxima cantidad de radiación que puede ser emitida. Tal emisor es el cuerpo negro. Se trata de la ley de Planck: Ebν(T)= 2πhν3n2 c02⋅(ehν σT−1) (2.10) En esta ecuación, h y k son las constantes de Planck y de Boltzmann. La constante de Planck nos aporta a través de ℎ⋅𝜈 la energía de un fotón. Sin embargo, en este trabajo se ha seguido la misma premisa que siguió Viskanta en la referencia [16]: “dos niveles diferentes de definición son posibles en este estudio de la interacción de la radiación electromagnética con materia: por un lado, el análisis de la interacción microscópica detallada de las ondas electromagnéticas son varios tipos de partículas desde el punto de vista de la mecánica cuántica; y por otro, el estudio macroscópico de la interacción entre la materia y la radiación […] En este informe, se considera el último punto de vista, es decir, el estudio del campo de radiación a través de ecuaciones macroscópicas, escritas en términos del coeficiente de absorción volumétrico promedio […]”. Por tanto, nos ceñimos a la deducción de la energía total de emisión del cuerpo negro, así como al estudio de las hipótesis de espesor óptico delgado y grueso, para comprender cómo interacciona la radiación con las ecuaciones de capa límite en convección forzada y convección natural. Volviendo a la ecuación (2.10), el índice de refracción n se refiere al medio que rodea el cuerpo negro. En nuestro caso, n toma valor unidad, como ya se justificó. Por otro lado, dicha ecuación puede escribirse también en función de la longitud de onda: ν=c0 λ → dν=−c0dλ λ2 ; Ebλdλ=−Ebνdν → Ebλ= 2πhc02 λ5⋅(ehc0 σλT−1) (2.11) Si el medio fuese un cristal o un líquido (𝑛≠1), [4] 𝑐0 debería sustituirse por 𝑐0/𝑛. La Figura 15 muestra la potencia emisiva espectral como una función de la longitud de onda para cuerpos negros a diferente temperatura. Por un lado, a mayor temperatura, la potencia emisiva se incrementa para una misma longitud de onda. Por otro lado, esta evolución, fijada una temperatura, presenta un máximo para una cierta longitud de onda, que puede obtenerse resolviendo numéricamente una ecuación transcendental (ver página 9 de [13]). El resultado se conoce Introducción a la radiación térmica 20 20 como Ley del Desplazamiento de Wien: λm=2898 Tμm (2.12) Una vez fijada la temperatura, la potencia emisiva decae para longitudes de onda mayores o menores que 𝜆𝑚. Según la ecuación (2.12), dicha longitud de onda 𝜆𝑚 decae a medida que crece la temperatura. Figura 15: energía específica del cuerpo negro en función de su temperatura y de la longitud de onda [13] Esta ley explica, por ejemplo, [3]por qué cuando un cuerpo metálico es calentado su superficie cambia de color gradualmente desde el rojo hasta el amarillo, e incluso hasta el blanco. De esta forma, si introducimos (2.11) en (2.8), obtenemos la potencia emisiva total de un cuerpo negro cuya superficie se encuentre a una temperatura T, según la ley de Stefan-Boltzmann: Eb=n2σT4 (2.13) Donde σ es la constante de Stefan Boltzmann (ver página 10 de [13] para más información sobre su obtención). Los cuerpos reales no son cuerpos negros, por lo que para medir la potencia emisiva espectral de una superficie real respecto a la del cuerpo negro a la misma temperatura, se emplea su emisividad espectral: Eλ(T)=ελ(T)⋅Ebλ(T) (2.14) Para superficies de diferentes materiales, 𝜀𝜆(𝑇)≤1 debe determinarse experimentalmente. La emisividad total se define como la ratio de la potencia emisiva de un cuerpo real respecto a la de un cuerpo negro a la misma temperatura. ε(T)=E(T) Eb(T)=∫𝜀λ(T)⋅Ebλ(T) ∞ 0σT4 (2.15) Al igual que 𝜀𝜆(𝑇), es una propiedad del material. Sin embargo, como se especificó en el apartado anterior, se trabajará con la hipótesis de cuerpo gris. Es decir, se hace un promedio de las propiedades antes mencionadas con la frecuencia y se considera un cuerpo gris "equivalente" con propiedades igual al promedio e independientes de la longitud de onda. Por tanto, de aquí en adelante se trabajará con magnitudes y ecuaciones que ya han sido integradas para todas las longitudes de onda. Por otro lado, como ya se mencionó, 𝛼=𝜀=1. De esta manera, se da por cerrado este apartado correspondiente al estudio de las condiciones de contorno del sistema que va a estudiarse en este trabajo. 21 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 2.3 Intensidad de radiación y ecuación de transferencia de radiación A continuación, pasamos a un estudio más detallado sobre la manera en la que se va a tratar la energía de los fotones para comprender la influencia que tiene en la propagación de calor debido a la radiación. Para ello, se va a repasar el concepto de ángulo sólido. El ángulo sólido es por definición similar al ángulo plano, pero viene dado en un espacio de tres dimensiones. La necesidad de usarlo en el aparato matemático da una idea de la dependencia del fenómeno de la radiación con el espacio. Mientras que el ángulo plano es la ratio entre un arco y la distancia entre el vértice del ángulo y un punto de dicho arco, [5] el ángulo sólido es la ratio entre un área contenido en una esfera y el cuadrado del radio de dicha esfera. Por ejemplo, un hemisferio tiene una superficie de 2𝜋𝑅2, por lo que, siendo la unidad del ángulo sólido el estereorradián (sr), tiene asociado un ángulo sólido de 2𝜋 𝑠𝑟 sobre el centro de su base. El estereorradián se define, al igual que el radián, a partir de un elemento unidad, en este caso, una esfera de radio unidad: es el área de valor unidad contenido en una esfera de radio también unidad, que se corresponde con un ángulo sólido de valor 1 sr. Por tanto, el ángulo sólido de una esfera completa será 4𝜋 𝑠𝑟. Véase a continuación el sistema de referencia provisional, así como el diferencial de ángulo sólido: Figura 16: hemisferio mostrando las dimensiones de un diferencial de ángulo sólido [5] Ya sea para caracterizar la cantidad de energía de radiación que emite una superficie a lo largo de un cierto camino, o para caracterizar la cantidad de energía de radiación que atraviesa localmente un elemento del medio fluido, [4] el concepto de un solo rayo es inadecuado. A cambio, se emplea el concepto de intensidad de radiación, que consiste en un modelo geométrico de un haz de fotones. Para el primer caso (desde una superficie) se emplea como símbolo “i” y para el segundo caso (en un medio) se emplea “I”. Para ser consistentes con las definiciones anteriores, se tiene: • 𝑖𝜈 (energía emitida desde una superficie por unidad de tiempo, por unidad de superficie normal al haz de fotones, por unidad de ángulo sólido y por unidad de frecuencia). • 𝑖 (energía emitida desde una superficie por unidad de tiempo, por unidad de superficie normal al haz de fotones, y por unidad de ángulo sólido). i=∫ iνdν ∞ 0 (2.16) Introducción a la radiación térmica 28 28 Dada la ecuación del gradiente del flujo de calor por radiación, se va a obtener la expresión aplicando la hipótesis de τ<<1, también denominada hipótesis thin. En caso de que τ<<1 [7]: E2(τ)=1+O(τ); E3(τ)=1/2−τ+O(τ2) (2.52) Siendo estos los términos, hasta de segundo orden, de la serie de potencias de 𝐸2(𝜏) y 𝐸3(𝜏). Considerando superficies negras, (2.38) queda: qR=e1(1−2τ)−e2(1−2τ0+2τ)+2∫ e(T)dτ′−2∫ e(T)dτ′ τ0 τ τ 0 (2.53) Evaluando los dos integrandos en los límites establecidos, y derivando la expresión, se tiene: −dqR dτ =2e1+2e2−4e(T) (2.54) La interpretación de (2.54) es como sigue. El tercer sumando del segundo miembro representa la energía emitida por unidad de volumen desde un elemento, mientras que el primer y segundo sumandos del segundo miembro denotan la energía absorbida por un elemento proveniente de la placa inferior y superior respectivamente. Por tanto, bajo condiciones de espesor óptico delgado todos los elementos del medio intercambian radiación directamente con las placas planas, porque, al ser λp relativamente grande, los fotones atraviesan el medio en estudio, sin que este sea capaz de incrementar su energía interna, hasta llegar a las placas. En general, se trata de un medio cuyos elementos tienen una capacidad despreciable de absorber/atenuar energía emitida/dispersada por los propios elementos del medio. A continuación, dada la ecuación del gradiente del flujo de calor por radiación, se va a razonar la obtención de la expresión aplicando la hipótesis de τ>>1, también denominada hipótesis thick, haciendo un comentario que involucra al estudio de la convección libre. Si el medio presenta un espesor óptico grueso, [13] las integrales exponenciales involucradas (de orden 2 y 3) decaen rápidamente a lo largo de un corto espesor desde 𝜏′=𝜏. Además, puede asumirse que, en una gran extensión del medio en estudio, nos encontramos alejados de las dos placas, a una gran distancia óptica. Bajo estas condiciones, la influencia de las superficies puede despreciarse. Por otro lado, sabiendo que 𝐸2(𝜏) decae rápidamente, desde 𝜏′=𝜏; y que la función S variará muy poco una vez 𝐸2(𝜏) sea despreciable, puede desarrollarse S mediante una serie de Taylor (ver demostración en páginas 451 y 452 de [13]). Despreciando términos de orden superior, se llega a: qR(τ)=−4π 3dS dτ; qR(y)=−4π 3κRdS dy (2.55) Para un medio que no dispersa radiación, 𝑆=𝐼𝑏, y la ecuación (2.55) se reduce a: qR(τ)=−4π 3dIb dτ=−4σ 3d(T4) dτ ; qR(y)=−4σ 3κRd(𝑇4) dy =−16σ𝑇3 3κRdT dy (2.56) El gradiente en dirección “y” del flujo de calor quedaría entonces: dqR dy(y)=−4σ 3κRd2(T4) d𝑦2=−4σ 3κRd dy[d(T4) dy ]=−16σ 3κRd dy[𝑇3dT dy]=−d dy[16σ𝑇3 3κRdT dy] (2.57) Con 𝜅𝑅 el coeficiente medio de absorción de Rosseland. Haya o no haya dispersión de la radiación, el coeficiente que se obtiene al realizar esta aproximación se denomina siempre de Rosseland, ya que fue el astrofísico noruego Svein Rosseland el que realizó esta aproximación referida en [13]. Cabe destacar que (2.57) es del mismo tipo 29 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre que la ley de Fourier de difusión de calor. Es decir, el problema se ha simplificado a uno de conducción con gran dependencia de la Temperatura. La razón radica en que [7] la emisión de las moléculas del medio se someterá a una gran atenuación de manera que el proceso de transferencia de calor dependerá completamente de las condiciones locales. Es decir, del choque de partículas del medio entre sí. La solución a este problema, sin embargo, no es válida cerca de las condiciones de contorno, ya que se requiere pequeños gradientes de temperatura, algo que no ocurre cerca de la placa plana, en este caso. De hecho, según [8], debe modificarse el coeficiente 16/3 por 8/3 cuando se calculen los flujos de calor de la pared al fluido, como se ve más adelante. En efecto, la solución nos muestra que [16] existe una discontinuidad entre la placa y el medio en contacto con ella. Solo en el caso de espesor óptico tendiendo a infinito, la temperatura exhibe continuidad cerca de la placa. Este hecho es consistente con el modelo de transferencia de energía mediante radiación únicamente, como se verá en los siguientes apartados. En la práctica este método es útil solo para casos con espesor óptico realmente grueso, [5] para valores de τ mayores que 10. Dicho valor depende de la geometría y de las condiciones del problema y podría ser mayor en algunos casos. Como ya se ha comentado, en convección libre se estudiará la distribución de velocidades y de temperaturas en el caso de que la capa límite generada sea ópticamente gruesa. De esta manera, es aplicable el método numérico y las reflexiones generales sobre la interacción de la radiación con el medio, en especial, todo lo concerniente a un medio ópticamente grueso. Por definición de esta hipótesis, las radiosidades asociadas a las dos placas planas no llegan a afectar a los alrededores de éstas. Por ello, se plantea el problema de la capa límite generada por convección libre en los mismos términos en los que se plantea el problema análogo en convección forzada: imponiendo las condiciones de contorno en temperatura en la pared y en el flujo ideal, sin entrar a valorar la influencia que pueda tener la presencia de las placas planas en la transferencia de radiación. El principal cambio será el del sistema de referencia (ver Figura 3), y en el sistema de ecuaciones que debe resolverse. 2.5 Hipótesis generales aplicadas a capas límite y a la radiación En este apartado se recopilan el total de las hipótesis que se aplican en este trabajo, tanto las que involucran a las capas límite, como las que involucran a la radiación térmica. Se incluyen también otras hipótesis que anteriormente no se mencionaron para no romper el ritmo de exposición. Las hipótesis no son excluyentes; es decir, tanto las que se aplican al modelo de capa límite como las que se aplican a la radiación, se aplican a un mismo flujo en estudio. Las hipótesis aplicadas a las capas límite en convección forzada y libre son: • Régimen estacionario (ambas). • Régimen subsónico. En consecuencia, régimen incompresible (ambas). • En convección forzada, Re>>1, Pe>>1. Los esfuerzos viscosos y de transmisión de calor se concentran en una zona muy estrecha adyacente a la pared del objeto. • En convección forzada, Fr>>1. • Ec<<1, con lo que se desprecian los efectos de disipación en la capa límite térmica debida a los efectos viscosos (ambas). Ec=u∞ 2 Cp∆T<<1 (2.58) • En convección forzada, en la ecuación de cantidad de movimiento: 𝜕2𝑢 𝜕𝑠2~𝑈∞𝐿2<<𝜕2𝑢 𝜕𝑦2~𝑈∞𝛿2 ⁄ ⁄ (2.59) Mientras que, en la ecuación de la energía: Introducción a la radiación térmica 30 30 𝑘𝜕2𝑇 𝜕𝑠2~𝑘𝑇∞−𝑇𝑠 𝐿2<<𝑘𝜕2𝑇 𝜕𝑦2~𝑘𝑇∞−𝑇𝑠 𝛿𝑇2 (2.60) Es decir, 𝛿/𝐿<<1 lleva a que todas las variaciones respecto a la variable “s” se desprecian, excepto los gradientes de presión impuestos por el flujo exterior. • Tanto en convección forzada como en convección libre, la ecuación de la cantidad de movimiento en dirección “y” se reduce en primera aproximación (ver [1] para más información) a: 𝜕𝑝 𝜕𝑦=0 (2.61) • En convección libre, se aplica la hipótesis de Boussinesq, por la cual, la densidad se toma constante en las ecuaciones de Navier-Stokes, excepto en la ecuación de la cantidad de movimiento en dirección longitudinal de la placa vertical, en cuyo caso: −𝜕𝑝 𝜕𝑥+𝜌𝑓𝑚 󰇍 󰇍 󰇍 󰇍  =−(𝜌−𝜌∞)𝑔 𝒊 (2.62) Las hipótesis aplicadas a la radiación son: • Dado que la velocidad de la materia, respecto a la de la luz, es despreciable, tomamos un sistema de referencia en reposo respecto a las placas. • El fluido que se estudia constituye un cuerpo gris. Por tanto, las propiedades involucradas son independientes de la longitud de onda. • El índice de refracción se toma 1. Según [4], es una buena aproximación para gases. Por tanto, la elección de dicho gas depende de la elección que hagamos del número de Prandtl. • No hay dispersión de la radiación. • Radiación difusa o isotrópica. • El medio se encuentra en equilibrio termodinámico local. Es decir, aunque las variaciones de las magnitudes puedan ser muy importantes desde un punto de vista macroscópico, su efecto desde el punto de vista de las variables de estado, en especial la temperatura, es muy pequeño. Así pues, cada elemento de materia, a una cierta temperatura, situado en cada punto del espacio se comporta como si formase parte de un entorno cerrado isotermo a dicha temperatura T. • Las placas planas constituyen cuerpos isotermos y negros, que absorben y emiten de manera difusa, y opacos (transmisibilidad nula). • La transmisión de radiación térmica se realiza perpendicularmente a las placas, unidimensionalmente. Dicho de otro modo, se busca demostrar que la transmisión de calor por radiación en la dirección “x” es despreciable frente al calor por convección en la dirección “x”, sabiendo que este último, bajo la mayoría de las condiciones, es un término predominante en la ecuación de la energía: u∂T ∂x>> 1 ρCp∂qRx ∂x (2.63) La convección en la dirección longitudinal se estima tal que: u∂T ∂x~U∞(Tw−T∞) L (2.64) Para estimar la expresión del flujo de calor por radiación, se tiene en cuenta que el límite de espesor óptico grueso, como se verá más adelante, [4] sobredimensiona la transferencia de calor por radiación si se aplicase a condiciones que no cumplen dicho límite. Según (2.63), ello permite que el límite de 31 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre espesor delgado pueda abrazar las conclusiones de este análisis de órdenes de magnitud. Por tanto: qRx=−16σT3 3κR∂T ∂x; dqR dx =~16σT3 3κR(Tw−T∞) L2 (2.65) De (2.63) y sabiendo que 𝛼=𝑘/𝜌𝐶𝑝 U∞L 16σT3 3κRρCp>>1 → Pe 16σT3 3κRk>>1 (2.66) El denominador constituye el inverso de Ncr, que indica el grado de importancia de la transferencia de calor por radiación respecto a la transferencia de calor por conducción. Dado que este análisis se plantea en el límite de espesor óptico grueso, donde se emplea el modelo de difusión de la energía, se espera que dicho coeficiente sea mucho menor que uno. Queda demostrado entonces que puede despreciarse la transmisión de calor por radiación en la dirección “x” respecto a la convección en la dirección “x”, siempre que Pe>>1, lo cual era ya una hipótesis necesaria para garantizar una capa límite adherida. • Como se mencionó antes, no es necesario considerar de manera detallada el mecanismo de interacción entre átomos y moléculas y el medio físico donde ocurren los fenómenos radiantes para deducir la ecuación de transferencia de radiación. Solo son necesarias las interacciones a nivel macroscópico entre la materia y el medio. En consecuencia, es innecesario retener la formulación de la teoría cuántica. Es suficiente con considerar el correspondiente haz de radiación, y seguir los cambios de intensidad a través de un medio en el que la temperatura varía de un punto a otro en el espacio. • No es un canal. Las dos placas están suficientemente alejadas físicamente hablando para que el flujo ideal no se vea afectado y se asegure el equilibrio termodinámico; pero lo suficientemente cerca ópticamente hablando para que la energía emitida por ambas placas afecte a la distribución de temperatura, para estudiar y comprender dicho efecto. De esta manera, por ejemplo, la placa plana más alejada, modificando sus propiedades, jugaría el papel de condición de contorno a elevada temperatura en un flujo hipersónico. • Si bien en un mismo flujo, sea convección forzada o convección libre, se transmite el calor a través de tres mecanismos diferentes, [12] la temperatura 𝑇∞ de dicho flujo se toma constante en todo momento, y en todo punto de la placa. Hecha esta doble introducción con los aspectos más interesantes para este proyecto de la convección forzada y natural; y de la radiación térmica, nos centramos en el apartado 3, en el cual, se explican los métodos numéricos empleados, y se aplican a las ecuaciones que se pretenden resolver. 32 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 33 3 Métodos numéricos empleados y formulación de ecuaciones 3.1 Métodos numéricos empleados En este proyecto se emplean fundamentalmente el método de las diferencias finitas y el método de Newton Raphson para resolver sistemas de ecuaciones no lineales. También se realiza un tratamiento numérico de expresiones integro-algebraicas, del cual puede encontrarse más información en [3] y en [18]. Como ya se ha dejado claro en el apartado 1, la capa límite es uno de los pilares del proyecto. En consecuencia, la resolución de las ecuaciones en este trabajo bebe de este concepto. De este modo, [2] la solución numérica de las ecuaciones de capa límite están basadas en la asunción de que las ecuaciones diferenciales, involucren o no a más de una variable dependiente pueden aproximarse por expresiones en diferencias. Esta aproximación, entendida como discretización, se obtiene como un desarrollo en serie de las derivadas de la temperatura o de la velocidad en el sistema de referencia. Es lo que se conoce como método de las diferencias finitas. Para obtener las expresiones en diferencias, debe situarse una malla en el espacio geométrico donde cada partícula tiene una velocidad y una temperatura incógnitas a resolver. Esta malla se compone de líneas paralelas a los ejes coordenados. Dichas incógnitas se asociarían a puntos de malla, es decir, puntos del espacio geométrico de coordenadas entendidos como la intersección de las líneas coordenadas en cuestión. La separación entre dichas líneas puede escogerse como variable o constante. Si se esperan grandes variaciones en los campos de velocidad o de temperatura, sería necesario que dicha separación se redujese, para que el error correspondiente a la discretización fuese mínimo: Figura 20: malla de discretización [2] En mecánica de fluidos, se suele adimensionalizar en primer lugar para disminuir el número de parámetros involucrados, para luego implementar métodos numéricos. Nada impide, aun así, que se aplique el método de las diferencias finitas a variables dimensionales, como “x” e “y”. Sin embargo, nuestro propósito es aplicarlo sobre variables adimensionales. A partir del apartado 3.2, introducimos los cambios de variable necesarios. Dicho esto, en la dirección “x” habitualmente se emplean Nxi líneas, siendo la variable adimensional ξ, y en la dirección “y” se emplean Neta líneas, siendo la variable adimensional η. En ocasiones se renombran como Ntau líneas, porque se emplea una malla que toma como variable τ (espesor óptico). En ese caso se deben modificar las matrices de transformación que a continuación se explican (ver anexos). Este método se aplica, no solo a ecuaciones lineales, sino también a ecuaciones no lineales, en cuyo caso, el sistema debe resolverse de manera iterativa. En este proyecto, ello se consigue implementando el método de Newton Raphson. Métodos numéricos empleados y formulación de ecuaciones 34 34 A la hora de discretizar las derivadas parciales mediante cocientes, hay tres posibilidades diferentes: expresión en diferencias hacia delante, expresión en diferencias hacia atrás, y expresión en diferencias centrada. Ver el trabajo fin de grado [19] para más información. En este proyecto se implementan una expresión en diferencias hacia atrás, y dos expresiones en diferencias centradas. Sea H la variable a resolver, y “x” e “y” las variables independientes. Despreciándose errores de segundo orden, y sabiendo que 𝐻𝑖,𝑗 tiene la naturaleza de un vector con componentes j-ésimas y en la estación i-ésima, se necesitan los siguientes cocientes para este proyecto: (∂H ∂x)i,j=Hi,j−Hi−1,j ∆x ; (∂H ∂y)i,j=Hi,j+1−Hi,j−1 2∆y (3.1) (∂2H ∂y2)i,j=Hi,j+1−2⋅Hi,j+Hi,j−1 (∆y)2 (3.2) Teniendo en mente los sistemas de referencia en la Figura 3 y la Figura 18, este método resuelve todas las componentes j-ésimas (perpendicular a la placa plana o forma general) en la estación i-ésima (a lo largo de la placa plana o forma general). Por tanto, Hi−1,j se calcula en la estación anterior. Sabiendo eso, (∂H ∂x)i,jtendrá un tratamiento más detallado en el apartado 3. De otra parte, los otros dos cocientes conviene que sean tratados como una transformación matricial, para facilitar la posterior aplicación del método de Newton Raphson: (∂H ∂y)i,j= [ 01 2∆y … −1 2∆y ⋱ ⋱ ⋮ ⋱ ⋱ ] ⋅Hi,j=D1⋅Hi,j (3.3) (∂2H ∂y2)i,j= [ −2 (∆y)21 (∆y)2… 1 (∆y)2⋱ ⋱ ⋮ ⋱ ⋱ ] ⋅Hi,j=D2⋅Hi,j (3.4) Cuando hablamos de ∆x, en nuestro proyecto se refiere como hξ. Realmente, no se van a deducir ecuaciones con ∆𝑦, ya que está dentro de las matrices D1 y D2. En cualquier caso, nos referiríamos a ello con hη o hτ. En los anexos, estas dos matrices se definen con el comando sparse al ser matrices banda. De este modo, dado que la capa límite para convección forzada en ausencia de radiación se trata ampliamente en la asignatura de “Mecánica de Fluidos II”, no se implementa (1.17) en ausencia de radiación ni se muestran resultados gráficos en el apartado 4 del campo de velocidades para placa plana. En cambio, en el apartado 3.3.3, sí se implementa (1.24), porque se emplean unos cambios de variable diferentes a los habituales para el campo de velocidades (el cambio habitual se implementa en el script 1). Para que conste, también se explica a continuación cómo se implementa el término del gradiente de velocidades, para ser considerado en proyectos similares, así como la condición de desprendimiento. En convección libre sí se considera pertinente implementar las ecuaciones de la capa límite mecánica y térmica en presencia de radiación, así como aportar resultados gráficos del campo de velocidades, dado que esta vez, sí depende del campo de temperaturas, siendo una novedad respecto a la asignatura de “Mecánica de Fluidos II”. En general, este proyecto aporta resultados para placa plana. Uno de los motivos es que generalizar a cualquier forma el cálculo de capas límite en presencia de radiación incurre en un aumento importante de la complejidad de las expresiones involucradas. Además, como se verá en el apartado 3.3, desde el punto de vista físico, la condición de contorno en temperatura en el exterior de la capa límite es más complicada de calcular a medida que el espesor óptico de la misma disminuye y deja pasar parte de la radiación. Sin embargo, en caso de espesor óptico grueso, la condición de contorno es simplemente la del exterior, 𝑇∞. Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 35 Profundizando en la condición de contorno en velocidad, hay dos maneras de calcularla: a través del coeficiente de presiones si el flujo es real, ver (3.5), o a través de una expresión analítica si el flujo es ideal, menos frecuente. En este trabajo se calcula la capa límite en espesor óptico grueso para una placa plana, un cilindro ideal, ver (3.6), y una elipse ideal, para simplificar. Se entiende por ξ como la variable adimensional que sigue el perfil: pe(s)+ρue 2(s) 2=p∞+ρU∞ 2 2 ; Cp=pe(s)−p∞ 1 2ρU∞ 2→ Cp+ue 2(s) U∞ 2=1 (3.5) ue(s) U∞=2∙sinξ (3.6) Una vez obtenida la velocidad exterior a la capa límite, se le aplica el método de las diferencias finitas al gradiente de velocidades tal que así: duedxi(2:Nxi)=(ue(2:Nxi)-ue(1:Nxi-1))./hxi; y se incorpora al término independiente de la ecuación numérica de la cantidad de movimiento: Figura 21: resolución numérica de la capa límite mecánica para una forma genérica; y condición de desprendimiento Para que conste en acta, y a pesar de que en la Figura 31 se muestran capas límite que no llegaron a desprenderse, en la Figura 21 se comenta que la condición de desprendimiento se cumple si la segunda componente del vector de velocidades horizontal es negativa. Se deja para futuros proyectos involucrar esta situación con la temática que se está abordando. Pasamos a tratar el método de Newton Raphson. Tal método se aplica en este proyecto para resolver una ecuación del tipo: 𝐁⋅θi 󰇍 󰇍 󰇍  +𝐂⋅θi4 󰇍 󰇍 󰇍 󰇍  =b 󰇍  (3.7) Siendo B y C matrices, y 𝜃𝑖 el vector temperatura adimensional incógnita en la estación i-ésima. Las matrices B y C deben completarse según unas condiciones de contorno que en los siguientes apartados se razonan. Métodos numéricos empleados y formulación de ecuaciones 36 36 De este modo, si θi 󰇍 󰇍 󰇍  =𝜃𝑘 󰇍 󰇍 󰇍 󰇍  +∆𝜃𝑘 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  , y despreciando términos de orden superior: 𝐁⋅(𝜃𝑘 󰇍 󰇍 󰇍 󰇍  +∆𝜃𝑘 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  )+𝐂⋅(θ𝑘 4 󰇍 󰇍 󰇍 󰇍  +4⋅θ𝑘 3 󰇍 󰇍 󰇍 󰇍  ∗∆𝜃𝑘 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  )=b 󰇍  (3.8) Usándose “*” como el producto componente a componente, también conocido como producto de Hadamard. De esta manera, y con el propósito de despejar ∆𝜃𝑘 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  se tiene: 𝐒⋅∆θk 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  =Fk, 𝐒=𝐁+4⋅𝐂⋅[θ13⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ θN 3], Fk=b 󰇍  −(𝐁⋅θk 󰇍 󰇍 󰇍 󰇍  +𝐂⋅θk 4 󰇍 󰇍 󰇍 󰇍  ) (3.9) ∆θk 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  =𝐒−1⋅Fk (3.10) Este método se itera hasta que la norma de ∆𝜃𝑘 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  es menor que una cierta cantidad. Véanse los anexos. En adelante, el propósito consiste en la obtención de las matrices B y C partiendo de las ecuaciones de la capa límite y del flujo de radiación explicados en los anteriores apartados. Antes de pasar a la resolución numérica de los problemas, [3] es necesario realizar transformaciones matriciales sobre expresiones integro-algebraicas tales como: F(τ)=∫ f(τ′)⋅E1(|τ−τ′|) τ0 0⋅dτ′ (3.11) De manera que puedan escribirse como: [F1 ⋮ FNz]=[Aij][f1 ⋮ fNz] (3.12) Para cada intervalo [τj,τj+1], con j= (1… 𝑁𝑧) una vez estamos en una estación i-ésima, la contribución a (3.11) se puede aproximar por: fj+fj+1 2∫ E1(|τi−τ′|) τj+1 τj⋅dτ′≡fj+fj+1 2gij (3.13) Recordando (2.37), y sabiendo que en (3.11) tenemos un valor absoluto en el argumento de la función exponcial, [18] dos casos se destinguen: E1(|τi−τ′|)=−dE2(τ′−τi) dτ′ (τi≤τj) y E1(|τi−τ′|)=dE2(τi−τ′) dτ′ (τi≥τj+1) (3.14) Si (3.14) se sustituye en la integral, 𝑔𝑖𝑗 toma la siguiente expresión compacta: gij=sign(τi−τj+1+0+)[E2(|τi−τj+1|)−E2(|τi−τj|)] (3.15) La cantidad arbitrariamente pequeña positiva 0+ se introduce en (3.15) para que la función sign aporte el signo correcto cuando 𝜏𝑖=𝜏𝑗+1. Por tanto, (3.11) y (3.13) resultan en: Fi=∑fj+fj+1 2⋅gij=f1 2gi1+f2 2(gi1+gi2)+...+fj 2(gij−1+gij)+...+fNz 2giNz−1 Nz−1 j=1 (3.16) y se rescata (3.12): Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 37 Fi=∑Aij Nz j=1 ⋅fj (3.17) donde: A(i,1)=1 2gi1, A(i,j)=gij−1+gij 2, A(i,Nz)=giNz−1 2 (3.18) Estas son las principales herramientas a nivel de métodos numéricos necesarias para resolver los problemas relacionados con convección forzada y convección libre. Concretamente, (3.17) se implementa en los scripts 2 y 4, que involucran un medio que absorbe y emite radiación. 3.2 Aplicación a la convección forzada. Espesor óptico grueso según Viskanta & Gross [8] Desde un punto de vista físico, es coherente que, a continuación, nuestros pasos se encaminen a la deducción de las ecuaciones que describen el flujo de un fluido no viscoso que no transmite calor por conducción, sino puramente mediante radiación. Sin embargo, dicho flujo adquiere en este trabajo una importancia mayor que la que su simplicidad desde un punto de vista conceptual pudiera sugerir. Efectivamente, un flujo viscoso se comporta, lejos de las capas límite, como un flujo ideal que, en primera aproximación, no transmite calor por conducción, y que además en este trabajo, transmite calor por radiación. Ello nos deja dos consecuencias. La primera es que, conceptualmente, el flujo ideal debería exponerse cuando exista interés en cuantificar el valor de la condición de contorno exterior a la capa límite; la segunda es que ello ocurre cuando el espesor de la capa límite no sea grueso. Cuando se introduzca el cambio de variable que dependa, entre otros, del espesor óptico de la capa límite, se darán más detalles. De momento, resaltar que la hipótesis thick se acerca más a la situación, ya estudiada en la asignatura de “Mecánica de Fluidos 2”, en la que no se tiene en cuenta la influencia de la radiación, y en la que se capturan los efectos del calor por conducción, que se propaga por difusión. Por ello, con la intención de estructurar el aparato matemático asentándonos en conocimientos previos, se considera en primer lugar el caso de espesor óptico grueso. Se tiene: ρ𝑐p(𝑢x∂T ∂x+vy∂T ∂y)=∂ ∂y[(16σT3 3β +k)∂T ∂y] (3.19) Se opta por denominar β a la absortividad con el objetivo de distinguirla de la conductividad térmica “k”. Se aplican los siguientes cambios de variable: x=ξ⋅L, y=η⋅L, T=T∗⋅θ, 𝑢x=U∞⋅u, vy=U∞⋅v (3.20) Rescatando los números adimensionales y magnitudes físicas antes mencionados: Re=U∞⋅L ν, Pr=να, α= k ρcp (3.21) Se manipula la expresión: U∞T∗ L(u∂θ ∂ξ+v∂θ ∂η)=ν⋅T∗ ν⋅L2∂ ∂η[(Pr⋅16σT∗3 Pr⋅3βρcpθ3+k ρcp)∂θ ∂η] (3.22) Métodos numéricos empleados y formulación de ecuaciones 44 44 óptico delgado, y la segunda integral se evalua usando la condición de contorno exterior (3.52). Reteniendo únicamente términos de primer orden, y recordando (2.54), la ecuación de la energía de la capa límite queda: ux∂T ∂x+vy∂T ∂y=α∂2T ∂y2+2βσ ρCp(T14+T∞4−2T4) (3.55) Es cierto que en el apartado 3.3.4, las ecuaciones se complican ya que es necesario tratar una expresión integral a través de una transformación matricial. Sin embargo, los cambios de variable que es preciso introducir a continuación son en esencia los mismos que los aplicados previamente, y serán los mismos en lo que queda de convección forzada. Solo se introducen otros dos cambios de variable para el campo de velocidades. Tenemos: ux=U∞⋅u; vy=U∞ √Re⋅v; v=v⋅√2⋅Pr⋅Ncr; θ=T T∗; x=ρCpU∞ 2βσT∗3⋅ξ=Lc⋅ξ; η=y⋅√U∞ ν⋅Lc=τβ⋅√U∞⋅2βσT∗3 ν⋅ρCpU∞=τ⋅√4σT∗3 ν⋅β⋅k α⋅2 ⁄→η= τ √2⋅Pr⋅Ncr (3.56) Si nos fijamos en (1.18) y (3.33), 𝜂~𝑂(1) si y=δ. Por tanto, se le da a η un valor máximo de 10. Cabe destacar que no solo se realizan los dos cambios de variable en η y en τ. También se realizan otros dos cambios que afectan al campo de velocidades. Partiendo de (1.15) y (1.16): U∞ Lc∂u ∂ξ+U∞ √ReLc∂v ∂η√ReLc Lc =0→∂u ∂ξ+∂v ∂τ=0 (3.57) U∞ 2 Lcu∂u ∂ξ+U∞ 2 √ReLcv∂u ∂η√ReLc Lc =ν∙U∞∙ReLc Lc2∙∂2u ∂η2→u∂u ∂ξ+v∂u ∂τ=2⋅Pr⋅Ncr⋅∂2u ∂τ2 (3.58) Con condiciones de contorno, correspondientes a placa plana τ=0 ,u=v=0; τ→∞,u→1 (3.59) El campo de velocidades es conceptualmente el mismo que el que se programa en el script 1, aunque se exprese en función de Pr, Ncr y el espesor óptico. Se trata de un mero artilugio matemático, en forma de cambio de variable, para que la resolución del sistema de ecuaciones (1.15)-(1.17) sea consistente. Por brevedad, no se le aplica en lo que sigue los métodos matemáticos a las ecuaciones de capa límite mecánica. Los códigos resultantes se encuentran en el script 4 para su consulta. En el apartado 3.4 sí se dedica más espacio a las ecuaciones, en vista de que, al ser convección libre, el campo de velocidades depende de la temperatura. Centrándonos en la ecuación de la energía, y saltándonos por brevedad los pasos intermedios de los cambios de variable para implementar el método numérico, se tiene: u∂θ ∂ξ+v∂θ ∂τ=2⋅Ncr⋅∂2θ ∂τ2−2θ4+θ14+θ0 4 (3.60) Deben modificarse las condiciones de contorno de (3.29) porque se ha producido un cambio en el sistema de referencia. Ahora se trabaja con el espesor óptico. Las condiciones de contorno son entonces: θ=θ0=T0 T∗, ξ=0,τ; θ=θcc, ξ,τ→∞; θ=θw, ξ,τ=0 (3.61) Aplicando diferencias finitas se tiene: Matu⋅(θn−θn−1)+Matv⋅D1⋅θn=2⋅Ncr⋅D2⋅θn−2⋅θn 4+θ14+θ0 4 (3.62) Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 45 Ordenando, y sabiendo (3.7): B=Matu+Matv⋅D1−2⋅Ncr⋅D2; C=2∗speye(Ntau); b=Matu⋅θn−1+ θ14+θ0 4 (3.63) Igual que antes, para lograr que el método matemático reconozca (3.61), se necesita completar las matrices B y C de la siguiente forma: B(1,:)=0, B(1,1)=1, C(1,:)=0, b(1)=θw B(Nz,:)=0, B(Nz,Nz)=1, C(Nz,:)=0, b(Nz)=θcc (3.64) El espesor físico de la capa límite mecánica puede medirse nuevamente con η. Sin embargo, la capa límite térmica se mide con espesor óptico. Tras tomar un valor máximo de η de 10, 𝜏𝐶𝐿 es suficientemente grande para capturar la evolución de la temperatura. Por otro lado, como se dijo previamente, en el flujo ideal, η no toma un valor concreto. Puede ser tan grande como queramos. En consecuencia, 𝜏0 y Ncr son dos grados de libertad. Ello favorece que comprobemos a qué valores de Ncr corresponde cada hipótesis, según los valores de 𝜏0 que correspondan a espesor delgado, grueso e intermedio, representando las curvas para cada caso (ver apartado 4.2). Sin embargo, cuando se trata de flujo viscoso, éste se confina en una capa límite de un espesor 𝜂~10. Por tanto, se reduce en un grado de libertad el cambio para η. En primera aproximación, se tiene Ncr<<1 para espesor óptico delgado: Figura 24: valores de Ncr para diferentes sustancias según la temperatura [4] Para el apartado 4.4, se toma como valor de Ncr intermedio un valor similar a los que aparece en la Figura 24, mientras que el valor de Ncr correspondiente a espesor óptico delgado se toma como aquel que hace que la condición de contorno se corresponda con la segunda componente τ-ésima del vector de temperaturas del flujo ideal. De ese modo, se puede comparar con los resultados correspondientes a espesor óptico delgado, donde, como venimos comentando en este apartado, se aplica la condición de Howell como condición de contorno. 3.3.4 Flujo viscoso. Relajación de hipótesis relativas a 𝛕. Resumen de ecuaciones Llegados a este punto, es necesario resolver de nuevo las ecuaciones del flujo viscoso correspondiente a espesor óptico grueso, pero esta vez con el mismo cambio de variable en la dimensión “x”, dependiente de Lc, que se viene realizando en este apartado 3.3, para poder establecer una comparación entre las curvas de 𝑁𝑢 √𝑅𝑒 ⁄ correspondientes a cada uno de los tres casos de flujo viscoso en estudio. En ocasiones, puede deshacerse cambios de variable a la hora de representar curvas con la intención de Métodos numéricos empleados y formulación de ecuaciones 46 46 compararlas con resultados dependientes de otras variables, como precisamente se hace en este trabajo al representar los perfiles de temperatura de la capa límite térmica, en el apartado 4.1, respecto a 𝜂𝑣. Sin embargo, cuando una variable es independiente, como en el caso de ξ, sea cual sea el cambio empleado, no es posible deshacer el cambio de variable, porque por definición no depende de otras variables adimensionales que puedan manipularse en este caso para facilitar la resolución numérica de la ecuación. Por ello, es necesario resolver el caso de espesor óptico grueso en función de Lc. También se pretende incorporar en la discusión una relación cualitativa, dada por Bo, entre el calor por convección y el calor por radiación puesto en juego en la capa límite que se comenta en el apartado 4.4. Empezamos por resolver las ecuaciones para espesor grueso en función del espesor óptico. Se tiene: ρ𝑐p(𝑢x∂T ∂x+vy∂T ∂y)=∂ ∂y[(16σT3 3β +k)∂T ∂y] (3.65) Nuevamente se aplica (3.56). Tras aplicar los dos cambios de variable, para η y para τ, queda: 𝑢𝜕𝜃 𝜕𝜉+𝑣𝜕𝜃 𝜕𝜏=2⋅𝑁𝑐𝑟⋅𝜕2𝜃 𝜕𝜏2+2 3⋅∂2(θ4) ∂𝜏2 (3.66) Se aplican diferencias finitas, y se ordenan los términos: B=Matu+Matv⋅D1−2⋅Ncr⋅D2; 𝐶=−2/3⋅D2; 𝑏=𝑀𝑎𝑡𝑢⋅𝜃𝑛−1 (3.67) Las condiciones de contorno son (3.61), con el cambio de 𝜃𝑐𝑐=𝜃∞. Pasamos a resolver las ecuaciones generales, en cuanto a que relajan las hipótesis relacionadas con el espesor óptico. Se tiene: 𝜌𝐶𝑝(u𝑥∂T ∂x+v𝑦∂T ∂y)=k∂2T ∂y2−4𝛽𝜎𝑇4+2𝛽𝜎[𝑇14𝐸2(𝜏)+∫𝑇4(𝑥,τ′) 𝜏𝐶𝐿 0E1(|τ−τ′|)dτ′+ +𝑇∞ 4𝐸2(𝜏−𝜏𝐶𝐿)] (3.68) Como puede verse, se va a integrar en el entorno óptico donde se ubica la capa límite, entre 0 y 𝜏𝐶𝐿. Al igual que se comentó al principio de este apartado 3.3, se debería integrar entre 0 e ∞. Sin embargo, el flujo que ocupa el intervalo entre 𝜏𝐶𝐿 e ∞ es ideal, afectado por la presencia de la capa límite. Ello repercutiría en un incremento de la complejidad del proceso de resolución, porque la condición de contorno exterior de la capa límite dependería del flujo ideal, que, al mismo tiempo, estaría influenciado por la presencia de la propia capa límite. Se recuerda que, en la práctica, en convección forzada se estudian por separado flujo ideal y flujo viscoso, es decir, cada uno se estudia con su propio volumen de control a la hora de aplicar las ecuaciones. Si bien el espesor físico de una capa límite puede despreciarse bajo ciertas condiciones, a estas alturas del trabajo puede verse que 𝜏𝐶𝐿 puede tomar valores del orden de 𝜏0. Por tanto, ninguno de los dos intervalos de la integral debería despreciarse. Estaríamos por tanto ante un proceso iterativo que involucraría expresiones integro algebraicos diferenciales parciales no lineales, como el que se propone en [9]. En el apartado 6 se plantea como línea futura. En este trabajo, por simplicidad, la condición de contorno exterior de la capa límite con espesor óptico 𝜏𝐶𝐿, nuevamente 𝜃𝑐𝑐, se toma como la componente más cercana a 𝜏𝐶𝐿 del vector de temperaturas del flujo ideal en la estación correspondiente. En el apartado 4.4, se dan más detalles. Por otro lado, se opta nuevamente por incluir la radiosidad de la placa superior, porque, al igual que en flujo ideal, los resultados que aporta el método numérico son más estables. Para que conste, aplicando los dos cambios de variable, el sistema de ecuaciones (1.15)-(1.17) más general posible queda: 𝜕𝑢 𝜕𝜉+𝜕𝑣 𝜕𝜏=0 (3.69) Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 47 𝑢𝜕𝑢 𝜕𝜉+𝑣𝜕𝑢 𝜕𝜏=2⋅𝑃𝑟⋅𝑁𝑐𝑟⋅𝜕2𝑢 𝜕𝜏2 (3.70) 𝑢𝜕𝜃 𝜕𝜉+𝑣𝜕𝜃 𝜕𝜏=2⋅𝑁𝑐𝑟⋅𝜕2𝜃 𝜕𝜏2−2𝜃4+𝜃14⋅𝐸2(𝜏)+𝐴𝑖𝑗⋅𝜃4+𝜃∞ 4⋅𝐸2(𝜏−𝜏𝐶𝐿) (3.71) Centrándonos en (3.71), tras aplicar diferencias finitas y ordenar términos, se tiene: B=Matu+Matv⋅D1−2⋅Ncr⋅D2; 𝐶=2∗𝑠𝑝𝑒𝑦𝑒(𝑁𝑡𝑎𝑢)−𝐴𝑖𝑗; 𝑏=𝑀𝑎𝑡𝑢⋅ 𝜃𝑛−1+𝜃14⋅𝐸2(𝜏)+𝜃∞ 4⋅𝐸2(𝜏−𝜏𝐶𝐿) (3.72) Se comenta brevemente, al igual que ocurría con el flujo no viscoso, que inspeccionando (3.60) y (3.63); (3.66) y (3.67); (3.71) y (3.72), en los tres casos de estudio, la matriz B es idéntica. Se corresponde con la parte lineal de las ecuaciones, incluyéndose los campos de velocidad conocidos o el transporte de calor por conducción. En cambio, la parte no lineal, C, es diferente, según la hipótesis en cada momento, y según se comenta en el apartado 2. El término independiente b también es diferente, según si las radiosidades de las placas estén o no atenuadas por la función exponencial (general o espesor delgado), o según si no se contemplan radiosidades de las placas (espesor grueso). Antes de pasar a convección libre, se hacen dos comentarios necesarios para analizar los resultados de los apartados 4.3 y 4.4. Estamos interesados en el Nusselt local para estudiar la evolución del fluido a lo largo de la placa, así como en el gradiente de temperaturas. De este modo, según [4]: 𝑁𝑢 √𝑅𝑒=𝜉1/2(𝜕𝜃 𝜕𝜂)𝜂=0,𝜃′=(𝜕𝜃 𝜕𝜂)𝜂=0 (3.73) Siendo θ la temperatura adimensional en uso. Los resultados respecto a 𝑁𝑢 √𝑅𝑒 y 𝜃′se expresan en función de (4.10) siempre que se compare con resultados de [4]. Por otro lado, en el apartado 4.3, estamos interesados nuevamente en expresar unos resultados respecto a 𝜂𝑣, como hacemos en el apartado 4.1. Sin embargo, el cambio en η es diferente, ver (3.33). En ese caso, 𝜂𝑣 se calcularía tal que así: 𝜂𝑣=𝑦⋅[1 2⋅𝑈 𝜈∗⋅𝑥]1 2=𝜂⋅𝐿𝑐 √𝑅𝑒𝐿𝑐⋅[1 2⋅𝑈∞ 𝜈∗⋅𝜉⋅𝐿𝑐]1 2=𝜂 √𝑅𝑒𝐿𝑐⋅[𝐿𝑐2 2⋅𝑈∞ 𝜈∗⋅𝜉⋅𝐿𝑐]1 2=𝜂⋅[1 2⋅𝜉]1 2 (3.74) En [4], dicho cambio de variable es ligeramente distinto, ya que no se expresa en función de 2. El propósito de dejar ese 2, es para uniformar el cambio de variable que realizan Viskanta&Gross y Sparrow. En cualquier caso, se demuestra que el resultado en la Figura 41 es óptimo. 3.4 Formulación de ecuaciones en convección libre. Hipótesis de espesor óptico grueso Para terminar el apartado 3, se van a implementar los métodos numéricos en las ecuaciones de la convección libre, incluyendo la novedad de que el medio fluido en reposo tiene un espesor óptico grueso. Al igual que en el apartado 3.2, se resuelven las ecuaciones con una malla en η y ξ; y con unos cambios de variable en función de la longitud de la placa L. El cambio en θ también es el mismo, con el propósito de obtener unos resultados gráficos a partir del mismo marco de referencia propuesto en [8]. La principal referencia usada en este apartado, y en el apartado 5, es [6]. En dicho trabajo fin de grado se resolvía un flujo de Rayleigh-Benard como éste, pero sin radiación. Así pues, partimos de las ecuaciones (3.33), (3.38) y (3.35): ∂u𝑥 ∂x+∂v𝑦 ∂y=0 (3.75) Métodos numéricos empleados y formulación de ecuaciones 48 48 ρ(u𝑥∂u𝑥 ∂x+v𝑦∂u𝑥 ∂y)=gβ(T−T∞)+μ∂2u𝑥 ∂y2 (3.76) ρCp(u𝑥∂T ∂x+v𝑦∂T ∂y)=∂ ∂y[(16σT3 3β +k)∂T ∂y] (3.77) Con condiciones de contorno: ux(x,0)=vy(x,0)=ux(x,y→∞)=0;T(0,y)=T∞ ,T(x,y→∞)=T∞,T(x,0)=Tw (3.78) Se realizan los siguientes cambios de variable: ξ=xL, η=yL, ux=u0∙u, vy=u0∙v,con u0=√gβL(Tw−T∞) , θ=T T∗ (3.79) Nuevamente, Ncr=κk 4σT∗3. Sabiendo que 𝐺𝑟=gβΔT⋅L3 ν2=(𝑢0∙𝐿 𝜈)2, se aplican entonces dichos cambios: ∂𝑢 ∂ξ+∂v ∂η=0 (3.80) 𝑢𝜕𝑢 𝜕𝜉+𝑣𝜕𝑢 𝜕𝜂=𝜈 𝑈∞∙𝐿∙𝜕2𝑢 𝜕𝜂2+𝑔𝛽𝑇∗𝐿 𝑔𝛽𝐿(𝑇𝑤−𝑇∞)(𝜃−𝜃∞) (3.81) 𝑢𝜕𝜃 𝜕𝜉+𝑣𝜕𝜃 𝜕𝜂=𝜈 𝑢0⋅𝐿⋅∂ ∂η[(16σT∗3𝜃3 3β𝜌𝐶𝑝𝜈+𝛼 𝜈)∂θ ∂η] (3.82) Dando como resultado: ∂𝑢 ∂ξ+∂v ∂η=0 (3.83) 𝑢𝜕𝑢 𝜕𝜉+𝑣𝜕𝑢 𝜕𝜂=1 𝐺𝑟1/2∙𝜕2𝑢 𝜕𝜂2+𝜃−𝜃∞ 𝜃𝑤−𝜃∞ (3.84) 𝑢𝜕𝜃 𝜕𝜉+𝑣𝜕𝜃 𝜕𝜂=1 𝐺𝑟1/2∙𝑃𝑟(1 3∙𝑁𝑐𝑟𝜕2𝜃4 𝜕𝜂2+𝜕2𝜃 𝜕𝜂2) (3.85) Se implementan diferencias finitas en (3.84) y (3.85). De (3.83) se obtiene v, partiendo de los resultados de u, y empleando el método de integración de los trapecios. Es la manera estándar de tratar numéricamente la ecuación de la masa en este proyecto (ver script 1, script 4 y script 5): un−1un−un−1 hξ+vn−1⋅D1⋅un=1 Gr1/2∙D2⋅un+θn−1−θ∞ θw−θ∞ (3.86) unθn−θn−1 hξ+vn⋅D1⋅θn=1 3∙Ncr∙Gr1/2∙Pr∙D2∙θn 4+1 Gr1/2∙Pr∙D2∙θn (3.87) De (3.86) cabe destacar las siguientes aproximaciones, con el propósito de iterar las estaciones precedentes: 𝑢𝑗(𝑥𝑛)≈𝑢𝑗(𝑥𝑛−1); 𝑣𝑗(𝑥𝑛)≈𝑣𝑗(𝑥𝑛−1); 𝜃𝑗(𝑥𝑛)≈𝜃𝑗(𝑥𝑛−1). A continuación, se agrupa (3.86) con el propósito de resolver un sistema lineal tal que 𝐴∙𝑢𝑛=𝑟𝑛; y se agrupa (3.87) con el propósito de resolverlo mediante Newton Raphson. En (3.86) se definen Matu y Matv según los resultados del campo de velocidades de las estaciones precedentes. En (3.87) se redefinen Matu y Matv como los resultados del campo de velocidades Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 49 de la estación actual, para calcular el campo de temperaturas (ver Figura 22): (Matu+Matv⋅D1−D2 Gr1/2)⋅un=un−1 2 hξ+θn−1−θ∞ θw−θ∞ (3.88) (Matu+ Matv⋅D1−1 Gr1/2∙Pr∙D2)⋅θn−1 3∙Ncr∙Gr1/2∙Pr∙D2∙θn 4=Matu⋅θn−1 (3.89) Las condiciones de contorno para (3.88) en términos de A y r son: 𝐴(1,:)=0; 𝐴(1,1)=1; 𝑟(1)=0; 𝐴(𝑁𝑒𝑡𝑎,:)=0; 𝐴(𝑁𝑒𝑡𝑎, 𝑁𝑒𝑡𝑎)=1; 𝑟(𝑁𝑒𝑡𝑎)=0; (3.90) Las condiciones de contorno para (3.89) en términos de B, C y b son: B(1,:)=0, B(1,1)=1, C(1,:)=0, b(1)=θw B(Nz,:)=0, B(Nz,Nz)=1, C(Nz,:)=0, b(Nz)=θ∞ (3.91) Por otro lado, según [6], δv/L~ 1 Gr1/4≪1,δt/L~ 1 Gr1/4∙Pr≪1 si el medio es completamente opaco a la radiación, es decir, si no se tiene en cuenta la radiación en el modelo. Según (3.31) y (3.82): δT L~√1+4 (3⋅Ncr) ⁄ Gr1/2 (3.92) Como se verá en el apartado 5, se toman dos valores diferentes de Pr, Pr=1, Pr=100. En consecuencia, (3.92) no presenta dependencia con la raíz de Pr en el denominador. Por un lado, para resolver el campo de velocidades y el campo de temperaturas al mismo tiempo, se necesita un mismo valor máximo de η. Por otro lado, si se incluyese la raíz de Pr en el denominador de (3.92), el campo de velocidades podría no tener cabida en el rango de η, al no presentar su espesor físico dicha dependencia con Pr. En consecuencia, se opta por elegir como cota superior de 𝜂𝑚𝑎𝑥 para la malla de diferencias finitas la expresión (3.92). Finalmente, se debe tomar como condición inicial de unm1=0.001 por ejemplo. La condición de vnm1 se queda nula. El motivo es puramente numérico, para evitar errores de condicionamiento en Matlab. Métodos numéricos empleados y formulación de ecuaciones 50 50 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 51 4 Efectos de la radiación térmica en convección forzada 4.1 Flujo viscoso. Capa límite ópticamente gruesa En este apartado se van a exponer los resultados, estableciendo una comparación con la referencia [8], obra de Viskanta y Grosh en 1962, correspondientes al caso de espesor óptico grueso. Dicho artículo presenta sus resultados para valor unidad de Pr. El cambio de variable que se ha efectuado permite recoger dichos resultados en función de la variable adimensional η multiplicada por un factor corrector dependiente de ξ, para poder comparar con [8]. Concretamente, dicha variable adimensional que emplea [8] es como sigue: 𝜂𝑣=𝑦⋅[(𝑚+1) 2𝑈 𝜈∗⋅𝑥]1 2 (4.1) La variable “m” es la variable de la que depende lo que [8] denomina β a través de: β=2m (m+1) ⁄ (4.2) Estrictamente hablando, [8] calcula la distribución de temperaturas de la capa límite sobre una cuña cuyo ángulo es πβ: Figura 25: sistema de ejes coordenados empleado en [8] En este trabajo se está interesado en la distribución de temperaturas en una placa plana, por lo que “m”, así como β, se anulan, quedando: 𝜂𝑣=𝑦⋅[1 2⋅𝑈 𝜈∗⋅𝑥]1 2 (4.3) Efectos de la radiación térmica en convección forzada 52 52 A su vez, U se toma como una solución de semejanza, expresándose como una constante multiplicada por una potencia m-ésima de x. Hay que destacar una vez más que uno de los propósitos de este trabajo es generalizar un método numérico, que combina diferencias finitas y el método de Newton Raphson, para calcular la distribución de temperaturas para cualquier forma en presencia de radiación, a través del cálculo de la capa límite mecánica. Más adelante se presenta el procedimiento a seguir en el caso general de cualquier forma. De momento, dado que “m” se anula, y en vista de que 𝜈∗ es entendida como la viscosidad cinemática a la temperatura de referencia 𝑇∗, en este trabajo se le da a U el mismo carácter de velocidad de referencia, tomándose como 𝑈∞. Introduciendo Re y las variables adimensionales ξ y η nos queda: 𝜂𝑣=𝜂𝐿⋅[1 2⋅𝑈∞ 𝜈∗⋅𝜉𝐿]1 2=𝜂⋅[𝐿2 2⋅𝑈∞ 𝜈∗⋅𝜉𝐿]1 2=𝜂⋅[𝑅𝑒 2⋅𝜉]1 2 (4.4) Este vector 𝜂𝑣 depende de una componente del vector ξ, con vistas a representar el perfil de temperatura en dicha estación. Comenzamos la exposición de los resultados con el perfil de temperaturas para un caso de pared fría (Figura 26) y otro de pared caliente (Figura 27). La naturaleza física de los resultados [8] puede visualizarse cuando se presta atención al valor del parámetro Ncr, que representa, como ya se ha comentado, la importancia relativa del transporte de energía mediante conducción y mediante radiación. Para el caso de Ncr=10, predomina la conducción de calor; en el caso de Ncr=1, el transporte por conducción molecular y por radiación son del mismo orden de magnitud, y en el caso de Ncr=0.1, predomina el transporte por radiación, y la hipótesis del espesor óptico grueso presenta menos aplicabilidad (puede verse cómo se separa claramente de las otras dos curvas para Ncr=10 y Ncr=1). Por otro lado, si bien para 𝑁𝑐𝑟→∞ el perfil de temperaturas se corresponde con el caso de conducción pura, la referencia [8] contabiliza en un 2% la desviación de todos los casos para Ncr=10, respecto al caso de conducción pura. De esta forma, el estudio de los perfiles de temperatura deja como conclusiones, que el efecto de la radiación en el caso de espesor óptico grueso es hacer más denso el flujo de la capa límite, dando pie a la aproximación del flujo de energía de radiación por difusión (efecto similar al que produce disminuir Pr, ya que implica aumentar la difusividad térmica) y que, dejando fijo un cierto valor de Ncr, en el caso de pared caliente, se aprecia, en el correspondiente perfil de temperatura, una mayor diferencia con el caso puramente conductivo que en el caso de pared fría. A continuación, pasamos a analizar el gradiente de temperatura en la última estación de la placa plana, según la ordenada vertical respecto a la pared 𝜂𝑣, tanto para pared fría como para pared caliente (Figura 28). Hay que destacar que [8] para el caso de pared fría, el gradiente de temperatura, que da una idea de la transmisión de calor por conducción, es máximo en la pared, mientras que para pared caliente el gradiente de temperatura es máximo a una ordenada vertical ciertamente alejada de la pared. Así, mientras que, en el caso de pared fría, una disminución de Ncr provoca que el valor absoluto del gradiente de temperatura más elevado, para cada valor de Ncr, se dé en cualquier caso en la propia pared, en el caso de pared caliente, una disminución de Ncr provoca que el valor absoluto del gradiente de temperatura más elevado, para cada valor de Ncr, se dé a una mayor ordenada vertical. Estas tendencias se deben a la remarcable no linearidad que presenta la dependencia con la temperatura del flujo de radiación. Por otro lado, en la zona adyacente a la placa plana, caracterizada por velocidades casi nulas de la capa límite mecánica, la transmisión de calor es principalmente por conducción. En dicha zona, en el caso de pared fría, debe destacarse que el gradiente de temperatura es constante si Ncr es grande, aunque deja de ser constante, y aumenta, si Ncr toma valores más pequeños como Ncr=0.1, en cuyo caso, la aplicabilidad de la aproximación difusa sería limitada. De esta manera, según (4.5), para el caso de pared fría el flujo de calor por conducción es máximo en la pared, y el flujo de calor por radiación es máximo a una ordenada vertical ciertamente alejada de la pared, cuando la temperatura se incrementa para adaptarse a la temperatura del infinito. Si la pared es caliente, en cambio, y debido a que la temperatura nunca volverá a ser tan grande como en la pared, el flujo de calor por radiación es máximo en la superficie, y el flujo de calor por conducción es máximo en una ordenada vertical ciertamente alejada de la pared. Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 53 Figura 26: comparación con [8] de la distribución de temperaturas en función de η escalada, en la última estación de la placa plana. 𝜃𝑤=0.1; 𝜃0=1; 𝑁𝑐𝑟={0.1,1,10} Figura 27: comparación con [8] de la distribución de temperaturas en función de η escalada, en la última estación de la placa plana. 𝜃𝑤=1; 𝜃0=0.1; 𝑁𝑐𝑟={0.1,1,10} Los gradientes de temperatura en la pared en función de 𝜃𝑤, si la pared es fría, y en función de 𝜃0, si la pared es caliente, se muestran en la Figura 29 y en la Figura 30 respectivamente. Para el caso puramente conductivo y pared fría, los gradientes de temperatura en la pared se desvían menos de un 1% respecto a Ncr=10, por lo que solo se muestra el caso puramente conductivo si la pared es caliente. La principal conclusión es que, si el transporte de energía está dominado por la conducción, los gradientes de temperatura en la pared presentan una tasa lineal en ambos casos, mientras que, si la transmisión de calor por radiación se convierte en el modo dominante de transmisión de calor, cuando disminuye Ncr, los gradientes de temperatura en la pared en ambos casos dejan de tener una tendencia lineal. Concretamente, si la pared es fría, la radiación incrementa los gradientes de temperatura, salvo que 𝜃𝑤 aumente, mientras que, si la pared es caliente, la presencia de la radiación disminuye el salto de temperatura para cualquier valor de 𝜃0. Efectos de la radiación térmica en convección forzada 60 60 Figura 34: distribución de temperaturas para flujo no viscoso con Ncr=0 (a, b y c), y con Ncr=1 (d, e y f); para tres posiciones en la placa plana, hasta ξ=1, Nxi=2000, y para τ=20 Figura 35: distribución de temperaturas para flujo no viscoso con Ncr=0 (a, b y c), y con Ncr=1 (d); para tres posiciones en la placa plana, hasta ξ=1, Nxi=2000, y para τ=0.2 Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 61 Figura 36: distribución de temperaturas para flujo no viscoso con Ncr=0 (a, b y c), y con Ncr=1 (d); para tres posiciones en la placa plana, hasta ξ=1, Nxi=2000, y para τ=0.002 Desde un punto de vista físico, [4] e independientemente de que se aplique o no la hipótesis de espesor óptico grueso o hipótesis thick, el transporte de energía radiativo (dirección vertical)-convectivo (dirección longitudinal) resulta en la formación de una capa límite térmica adyacente a la placa plana, de manera que los gradientes de temperatura se encuentran confinados en esta zona. Como puede apreciarse en la Figura 34, el requisito para que el modelo general recoja dichos gradientes es un gran espesor óptico, que en este trabajo se toma como τ=20. Por otro lado, recordando (2.50), τ>>1 radica en λp<<1, para que cada molécula le pueda transmitir su energía a su vecindad, y se pueda generar dicha distribución de temperatura. En consecuencia, estos gradientes también pueden reproducirse a través de la hipótesis thick, tratada en este trabajo, junto con la hipótesis thin, para el caso de flujo ideal en el apartado 3.3.1. Dado que la hipótesis thick constituye una aproximación por difusión, su radio de aplicación se amplía hasta ξ>>1, especialmente cuando Ncr=1, ya que incorporar al modelo general conducción de calor refuerza la transmisión de energía por difusión. Esta idea se torna más decisiva si cabe en el marco de las capas límite viscosas, donde por definición de la viscosidad de un fluido, se transmite energía por conducción. De momento, en la Figura 34, con flujo ideal y Ncr=1, cuando ξ=1, la desviación de la hipótesis thick no supera el 1% respecto al caso general. Del mismo modo, la aproximación por difusión permite que, efectivamente, el salto de temperatura que refleja radiation slip sea menor, y que prácticamente desaparezca conforme avanzamos en el flujo. En cambio, en las primeras estaciones, la aplicabilidad de esta hipótesis se reduce ligeramente. En las otras dos gráficas se aprecia también que, conforme se reduce el espesor óptico, la hipótesis thick no permite capturar el comportamiento térmico del fluido, incluso aunque se incremente Ncr. De hecho, para τ=0.2 y τ=0.002, el flujo resuelto con hipótesis thick alcanza el régimen estacionario mucho antes de alcanzar ξ=1. Por tanto, la primera conclusión de este trabajo es que el espesor óptico supone la primera variable a determinar para pasar de hipótesis thick a thin. No en vano, en la Figura 35, puede verse que el flujo general parece adaptarse mejor a la hipótesis thin cuando τ=0.2, sobre todo en las primeras estaciones. Cuando τ=0.002, el flujo general se adapta a la hipótesis thin con una desviación de menos del 0.1%. La segunda conclusión, como puede verse en la Figura 34, es que Ncr, modula cuánta distancia recorre el flujo general para comportarse según la hipótesis thick. En los siguientes apartados se va en esta dirección, en el marco de las capas límite viscosas. Efectos de la radiación térmica en convección forzada 62 62 4.3 Flujo viscoso. Capa límite thin. Revisión de la capa límite thick Figura 37: comparación entre segunda componente del vector de temperaturas del flujo ideal designado, y condición de contorno exterior de Howell para el caso específico de hipótesis thin en capas límite viscosas A continuación, se comentan los resultados que involucran a la distribución de temperaturas en las capas límite viscosas, estando aplicada la hipótesis thin. También se revisa la capa límite thick, resuelta esta vez en función de (3.56). Al igual que, previamente, la componente del vector de temperatura cercana a la placa debía expresarse a través de (4.7), esta vez, es el vector de temperatura completo lo que debe expresarse en función de dicho cambio de variable, para comparar con los resultados de [4]: Θ(ξ,η)=𝑇(𝜉,𝜂)−𝑇𝑤 𝑇∞−𝑇𝑤=𝜃(𝜉,𝜂)−𝜃𝑤 𝜃∞−𝜃𝑤 (4.10) En la Figura 37 puede apreciarse la relativa exactitud con la que la condición de Howell sigue la condición de contorno correspondiente al flujo no viscoso no conductivo que se toma como flujo exterior a la capa límite viscosa. Concretamente se toma la segunda componente. Como se comentaba en el apartado 3.3.3, cuando el espesor óptico es delgado, la capa límite deja pasar la mayor parte de la radiación térmica, entrando en juego la energía radiada por las placas, y la propia emisión del medio, sin ser atenuadas por éste. Se recuerda además que el parámetro conducción-radiación, en este contexto, conviene que tome valores tal que Ncr<<1, para evitar cualquier indicio de difusión que implique alejarnos del caso thin. Por tanto, estando a merced de la exactitud del método numerico empleado, el espesor óptico de la capa límite que no puede ser más delgada se correspondería con la segunda componente del flujo ideal exterior. De este modo, es posible comparar esta componente con la condición de Howell, la cual, para este trabajo, se iteró unas 1000 veces para lograr una aproximación lo más exacta posible: Figura 38: código que implementa la condición de contorno de Howell Para que conste, en la pared fría se ha mostrado cómo son las iteraciones de “n” impar y “n” par. Puede verse, tanto para pared caliente como para pared fría, que la aproximación es válida hasta 𝜉~0.5. Sabiendo (3.56), y estimando tal que 𝜌=1.223 𝑘𝑔/𝑚3,𝐶𝑝=1005 𝐽/𝑘𝑔⋅𝐾,𝑇∗=289 𝐾,𝛽=0.5,𝐿𝑐∼1𝑘𝑚. En función de Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 63 la aplicación ingenieril en la que estemos interesados en aplicar estas estimaciones, puede valorarse si durante unos 500 metros es deseable que el flujo se comporte con espesor óptico delgado. Tabla 4: comparación de la distribución de temperaturas θ en algunas estaciones de la capa límite de espesor thin en el caso de Ncr<<1 y de Ncr=2 En el apartado 4.1 había interés en dar valores diferentes a Ncr para saber cómo el flujo de temperaturas evolucionaba con 𝜂𝑣 en una capa límite con espesor óptico grueso. De la misma forma, también es interesante hacer lo mismo en una capa límite con espesor óptico delgado. De ahí la razón de ser de la Tabla 4: vemos que, en este último caso, por más que aumentemos Ncr, la evolución es la misma en las componentes escogidas. Efectivamente, la hipótesis sobre el espesor óptico determina la dependencia con el parámetro Ncr. Así mismo, la hipótesis sobre el espesor óptico sería más o menos acertada según el valor de τ en la ecuación de la energía general del flujo, en este caso, de la capa límite, como se deja entrever al final del apartado 4.2. Figura 39: evolución de 𝑁𝑢 √𝑅𝑒 con ξ para pared fría. Comparación con [4] Para pared caliente no se obtienen curvas tan estables como las que se muestran en la Figura 39. No se descarta que, en ese caso, al igual que se comentó en el contexto de la capa límite de espesor thick, sea necesario tomar alguna consideración cerca de la placa, como, por ejemplo, relajar la hipotesis thin, y plantear la posibilidad de que el medio pudiese ser participativo, debido a la energía puesta en juego. Volviendo a la Figura 39, cabe destacar la veracidad del resultado hasta 𝜉=0.1, valor longitudinal que se encuentra dentro del rango en el que la condición de contorno de Howell era realista. Como se muestra en el apartado 4.4, hasta ese punto la solución que relaja las hipótesis thin y thick se adapta muy bien a la solución thin si el valor de Ncr es adecuado. Por otro lado, merece la pena verter unos comentarios acerca de las implicaciones físicas para la capa límite ópticamente delgada. Como puede verse en la Figura 39, para ξ≪1, [4] el valor límite de 𝑁𝑢/√𝑅𝑒∼0.332, y representa un intercambio de energía de radiación despreciable. Ello es lógico, sabiendo que, por (1.11), ξ es una medida relativa de la proporción de la radiación frente a la convección en capas límite viscosas. Por tanto, las primeras estaciones dejan a la convección como el modo dominante. Más adelante, la convección pierde terreno de forma cualitativa, dando lugar a un balance entre conducción y radiación. Según [4], la máxima transmisión de calor por conducción que podría lograrse con el intercambio de calor por radiación puede ser del doble respecto a la conducción de las primeras estaciones. Tal afirmación no se puede verificar con el resultado original, ya que a partir de 𝜉=0.1 se reduce la aplicabilidad de la hipótesis de espesor óptico delgado. En la Figura 46 se realizan estas reflexiones sobre 𝑁𝑢/√𝑅𝑒 frente a ξ para pared caliente y fría y cada hipótesis. η\ξ 100, Ncr<<1 100, Ncr=2 500, Ncr<<1 500, Ncr=2 1000, Ncr<<1 1000, Ncr=2 100 0,17777416605 0,17777416604 0,65482314726 0,65482314726 0,81982181382 0,81982181382 300 0,14950362516 0,14950362516 0,34892879892 0,34892879892 0,60058214654 0,60058214654 Efectos de la radiación térmica en convección forzada 64 64 Figura 40: evolución que describe el gradiente de temperaturas Θ𝑤 ′ en función de 𝜃𝑤/𝜃0 para cuatro estaciones diferentes de la placa plana. Comparación con [4]. Respecto a la variación de Θ𝑤 ′ en función de 𝜃𝑤/𝜃0 (una vez que se adimensionaliza cociente y denominador por la temperatura de referencia), cabe destacar que el resultado original calcula la capa límite para obtener ese resultado hasta ξ=1. El motivo es que, en la Figura 39, se demuestra que más allá de ξ=1 el resultado no reproduce la gráfica en [4], por lo que la comparación se prefiere realizar hasta el punto donde el resultado es fiable. Dicho esto, la Figura 40 se obtiene calculando en primer lugar Θ𝑤 ′ según 𝜃0=1, y dándole valores cada vez más grande desde 0.1 a 𝜃𝑤; y luego según 𝜃𝑤=1, y dándole valores cada vez más pequeños desde 0.9 a 𝜃0. Vemos que el máximo incremento de conducción ocurre entre 𝜃𝑤/𝜃0∼0.5 y 1. Según avanzamos en la placa, los gradientes, en general, tienden a disminuir, sin que tenga alguna incidencia el valor de Ncr, a diferencia del caso thick. En el resultado original, el gradiente no lleva a ser negativo, si bien la tendencia es que, conforme avanzamos hacia un ξ mayor, los gradientes tiendan a ser negativos. Es decir, la presencia de la radiación produce un incremento o un decremento de conducción según el valor de 𝜃𝑤/𝜃0, fruta de las no linealidades presentes. Figura 41: comparación de la distribución de temperaturas en capa límite viscosa con espesor thin, con Pr=1 y con Nxi=2000 resuelta con códigos originales y con consideraciones expuestas en [4] En la Figura 41 podemos apreciar la capa límite térmica para dos estaciones situadas al principio de la placa plana: para n=5 y para n=200. En el primer caso, ξ<<1, y en el segundo, ξ=0.1. En virtud de la condición de Howell por la que se opta en este trabajo para reproducir la condición de contorno exterior para el caso thin, que emplea [5], se tiene que el flujo evoluciona de manera análoga a cómo evoluciona en [4]. La temperatura exterior de la capa límite evoluciona desde la temperatura en el infinito hasta una temperatura que se acerca a la temperatura de la pared. Llegados a este punto, cabe destacar que el salto de temperaturas que ocurre en la pared en el flujo ideal nos permite, dentro del orden en el que se trabaja en este proyecto, poder relacionarlo con la Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 65 condición de contorno exterior de la capa límite que relaja hipótesis. De esta manera, se evita unas ligeras desviaciones que la Figura 41 muestra para la referencia [4] cuando η alcanza las condiciones exteriores. Para el siguiente apartado, es necesario recalcular el flujo de temperaturas para el caso thick en función de las mismas variables τ y ξ con el propósito de comparar dichos resultados con aquellos que se obtienen con hipótesis relajadas. Como puede verse, los resultados son idénticos al caso thick en función de η y ξ, con igual número de componentes verticales y horizontales en el método numérico. Tabla 5: comparación de la distribución de temperatura θ en capa límite viscosa thick para algunas estaciones resuelta en función de τ y ξ (ver apartado 3.3.4); y resuelta en función de Re y η (script 1) 4.4 Flujo viscoso. Condiciones para asumir un límite óptico. Conclusiones En este apartado se pretende una síntesis de las conclusiones más importantes en convección forzada en términos de τ, Ncr y espesor físico η. Comenzamos por describir el flujo ideal. Se decide 𝜏0=20, pared fría (𝜃𝑤= 0.1,𝜃0=1), Nxi=2000 y Ntau=1000 en script 2. Se muestra la segunda componente de θ de todas las estaciones en ξ de dicho flujo en la Figura 37. Se guardan los resultados de θ para cada 𝜏𝑗 y cada 𝜉𝑖 de la malla en ‘thetaccideal20fria.mat’, ver script 4. Dicho esto, se muestran ahora resultados análogos a los que se muestran en el apartado 4.2 con flujo ideal. Éstos tenían la razón de ser de justificar que la hipótesis depende de τ, y la evolución del flujo en función de la conducción de calor depende de Ncr. La Tabla 4 iba en la misma línea. Esta vez, el propósito es demostrar que en la capa límite con espesor óptico variable ocurre lo mismo, con el matiz de que este flujo conduce calor siempre. En consecuencia, las condiciones de espesor óptico delgado se ven comprometidas conforme el flujo avanza a lo largo de la placa, sea cual sea el valor de Ncr. Figura 42: distribución de temperaturas en función de τ de las tres posibles soluciones originales, relajando hipótesis (𝑓(𝜏)), hipótesis thick e hipótesis thin para Ncr=2, Ntau=1000 y Nxi=2000. En este caso se produce en torno al 0.01% de error entre flujo general y thick en ξ=1, ópticamente cerca de la placa n(τ ó η)\n(ξ) 100(con τ) 100(con η) 300(con τ) 300(con η) 1000(con τ) 1000(con η) 100 0,8674662645662 0,8674662645665 0,5662610141026 0,5662610141032 0,4462613819049 0,4462613819061 300 0,9999560087605 0,9999560087605 0,9513512970244 0,9513512970245 0,8590854875565 0,8590854875567 Efectos de la radiación térmica en convección forzada 66 66 Figura 43: distribución de temperaturas en función del espesor óptico de f(τ)), thick y thin, para Ncr=2⋅10−6, Ntau=1000 y Nxi=2000. En este caso se produce menos del 0.01% de error entre flujo general y thin en ξ=0.0045023. Figura 44: distribución de temperaturas en función del espesor óptico de f(τ), thick y thin para Ncr=0.1, Ntau=1000 y Nxi=2000; y superposición con la distribución de temperaturas de flujo ideal Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 67 De esta forma, [10] a medida que crece 𝜉≪1, la radiación tiende a reducir la temperatura en la capa límite y un gradiente de temperaturas se genera en la placa. La magnitud de dicho gradiente de temperaturas se controla por un balance entre el flujo de radiación desde la placa y el flujo de calor por conducción hacia la superficie. En esencia, para pequeños valores de 𝜉, y como se comenta anteriormente, se cumple la hipótesis de espesor óptico delgado, mientras que para 𝜉~𝑂(1) o mayor que uno, el flujo se asemeja a la hipótesis de espesor óptico grueso, ya que se da más espacio y distancia para que las partículas absorban la energía de los fotones, se exciten y se produzca difusión de energía. Esta conclusión es la misma que se extrae de flujo ideal, si bien, para un mismo valor de Ncr, el flujo necesita una mayor ξ, así como un τ grande, para que se alcance un nivel de difusión que case con la hipótesis por la ausencia de la capa límite. Evidentemente, el valor del espesor óptico predispone la opacidad de la capa límite frente a la radiación, y el valor de Ncr, si se alcanzan las condiciones ópticamente gruesas para un valor grande o pequeño de ξ. Pero la evolución que los perfiles de temperatura tienen a lo largo de la placa son independientes de las condiciones del experimento. Figura 45: distribución de temperaturas en función de η de f(τ) para los tres valores típicos de Ncr, Ntau=1000 y Nxi=2000 A continuación, se comenta la Figura 45, que versa sobre la evolución del espesor físico de las capas límite. A ese respecto, huelga decir que, para las estaciones más tempranas, la evolución, como puede verse, es similar. Para las siguientes estaciones, se aprecian dos ideas. La primera es que el espesor físico para hipótesis thin, es decir, Ncr<<1, es cualitativamente menor que el espesor físico para hipótesis thick, es decir, para Ncr>>1, y relajando hipótesis, porque, si bien en la capa límite thin los fotones pasan de largo con más frecuencia sin ser absorbidos por las partículas fluidas, se transmite energía en forma de radiación unos órdenes superior a la energía transmitida en forma de conducción, de forma que el flujo alcanza la temperatura exterior cuando aún hay partículas fluidas que no han alcanzado la velocidad 𝑈∞ ideal ( ver apartado 1.2). En cambio, el espesor físico para hipótesis thick es mayor, porque, esta vez, la energía por radiación se admite que se propaga por difusión, de manera que las partículas fluidas excitan por contacto a sus vecinas y la transmisión de energía se produce de una manera más suave. En este caso, cuando las partículas circulan a velocidad 𝑈∞ ideal, aún no han Efectos de la radiación térmica en convección forzada 68 68 alcanzado la temperatura exterior, que, por tener la capa límite un espesor óptico grande, tiende a ser 𝑇∞. La conclusión, por tanto, para Ncr<<1 es que la radiación se revela como un eficiente mecanismo de transmisión de energía. Esta idea puede verse desde otro punto de vista. Según (3.31), Si 𝑁𝑐𝑟→∞, es decir, si se tiene conducción pura, 𝛿𝑇/𝐿~𝑅𝑒−1/2𝑃𝑟−1/2<<1. Pero, si 𝑁𝑐𝑟→0, es decir, si tenemos tanto conducción como radiación, y cada vez se le otorga más importancia a la radiación, 𝛿𝑇/𝐿 crece. Aparentemente, si Ncr<<1, 𝛿𝑇/𝐿 es mayor que en el caso thick. Sin embargo, (3.31) es válida solo para condiciones thick, de tal forma, que usarla para Ncr<<1 implica admitir que la energía por radiación térmica se comporta, en general, de manera difusiva, afectando solo a la vecindad de las partículas fluidas, cuando, en realidad, como se puede ver, es mucho más energética si Ncr<<1. Por tanto, admitiendo Pr del orden unidad, según la Figura 26 y Figura 45 los órdenes de magnitud quedan, en función de los valores de Ncr: 𝛿/𝐿 (𝑚𝑒𝑐á𝑛𝑖𝑐𝑎)~𝛿𝑇/𝐿 (𝑁𝑐𝑟→∞)<𝛿𝑇/𝐿 (𝑁𝑐𝑟~𝑂(1))<𝛿𝑇/𝐿 (𝑁𝑐𝑟<1 𝑠𝑖 𝑡ℎ𝑖𝑐𝑘)> 𝛿𝑇/𝐿 (𝑁𝑐𝑟<<1) (4.11) Dicho de otro modo, cuando el parámetro Ncr se aleja de situaciones propicias para que se dé difusión térmica en el medio, y la radiación es varios órdenes más importantes que la conducción, entonces la radiación térmica se convierte en un eficiente mecanismo de transmisión de energía, al ser 𝛿𝑇/𝐿 menor que si se diese difusión térmica. Figura 46: evolución de 𝑁𝑢 √𝑅𝑒 con ξ en el caso de pared fría y pared caliente, y los tres valores típicos de Ncr Por otro lado, en la Figura 46, se pretende realizar un análisis similar al de la Figura 39. ¿Cómo varía cualitativamente la transmisión de calor por conducción a la placa (en caso de placa fría), y la transmisión de calor por conducción desde la placa (en caso de placa caliente)? En el primer caso, vemos que, si el medio es ópticamente grueso, la transmisión de calor por conducción a la placa crece lentamente. Conforme se reduce Ncr y pasamos a espesor delgado, la transmisión de calor por conducción crece más claramente. En el segundo caso, vemos que la transmisión de calor por conducción tiende incluso a disminuir conforme se reduce Ncr. Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 69 Figura 47: gradiente de temperaturas frente a 𝜂𝑣 calculado a partir de f(τ) en las estaciones indicadas en la leyenda, para los valores típicos de Ncr Finalmente, la Figura 47 pretende analizarse basándonos en la Figura 28 anterior. Por un lado, vemos que, para pared fría, se reproduce con cierta exactitud las curvas correspondientes. Sin embargo, para pared caliente, las curvas son notoriamente diferentes. Como se comenta en el apartado 4.1, aparentemente con pared caliente debe corregirse el flujo de radiación, porque se acentúa la discontinuidad propia del flujo con espesor óptico cercano a la placa. Aquí concluye el apartado 4, así como las reflexiones sobre convección forzada. En el apartado 5 no se pretende añadir más información sobre convección forzada, sino realizar un estudio introductorio sobre convección libre del medio típico de este proyecto con espesor óptico grueso. En el apartado 6 se escribirá una síntesis de las conclusiones más importantes del proyecto. Efectos de la radiación térmica en convección libre. Espesor óptico grueso 76 76 77 6 Motivación, conclusiones y futuras líneas La motivación última que ha impulsado al autor de este proyecto a su realización es la comprensión de los fenómenos que ocurren en la reentrada en la atmósfera terrestre de una cápsula en régimen hipersónico, destacando efectos de compresibilidad, radiación térmica, e ionización. Una vez comprendida la magnitud de la empresa, se partió de trabajos ya existentes para estudiar efectos de la radiación térmica en capas límite, generadas por convección forzada (en régimen subsónico) y por convección libre. De esta manera, se pretende valorar esta aportación al departamento de Mecánica de Fluidos de la Universidad de Sevilla como parte de un conocimiento en contínuo desarrollo que daría la oportunidad en futuros proyectos de comprender fenómenos con los que las agencias espaciales deben lidiar para, por ejemplo, permitir regresar a la Tierra misiones tripuladas en condiciones óptimas de seguridad, y con los que deberán lidiar para afrontar futuras misiones tripuladas a Marte, que también presenta una atmósfera. Dicho esto, este trabajo fin de grado deja como conclusiones las siguientes: • Tanto un flujo ideal como un flujo viscoso de materia llena de partículas, puede ser un medio ópticamente grueso u ópticamente delgado. • Tanto en un flujo ideal como en un flujo viscoso a través de una placa plana, se repite el siguiente patrón: en las primeras estaciones de la placa plana el flujo tiende a filtrar la radiación porque los fotones aún no han tenido espacio para ser absorbidos, incluso aunque 𝜆𝑝≪1; y en las siguientes estaciones, en general, los fotones cada vez serán absorbidos más frecuentemente. Esa es la razón por la que las figuras gráficas, que incluían perfiles de temperatura con los tres casos posibles, y proponía condiciones que favorecían la hipótesis thin, mostraban que, en las primeras estaciones, el flujo general se comportaba rápidamente como le correspondía; y las figuras gráficas que proponían condiciones que favorecían la hipótesis thick, mostraban que las últimas estaciones, era cuando el flujo general se comportaba como correspondía. Según si se transmitía el calor por conducción en cada caso, podía asemejarse más rápidamente. Por ejemplo, en la Figura 34 se muestra que, en flujo ideal, cuando Ncr=1, y se tenían condiciones thick, en las primeras estaciones, en el flujo general las partículas aún no transmitían calor por difusión tan claramente como en las últimas estaciones. • En una capa límite mecánica se transmite por necesidad calor por conducción debido a su naturaleza física. La consecuencia de dicha naturaleza es un flujo confinado que hace que Ncr y τ dependan uno del otro. Al mismo tiempo, analizar la opacidad de la capa límite está directamente relacionado con el valor que toma el parámetro Ncr. En cambio, en flujo ideal, Ncr y τ son dos grados de libertad, y puede estudiarse por separado aquellos valores de Ncr que casan con según qué valores de τ. La regla general, es que, si un medio es ópticamente grueso, las partículas tienden a transmitir la energía de manera difusiva. Por tanto, si en dicho medio, se favorece la transmisión de calor por conducción, de manera indirecta el medio se comportará como ópticamente grueso. Eso es casar un alto valor de τ y de Ncr. En cambio, si un medio es ópticamente delgado, se favorece la transmisión de calor por radiación: esta se filtra por el medio. Para que eso sea posible, no es posible que se “tropiece” con partículas que se encuentre por su camino. En consecuencia, casan un valor pequeño de τ y un valor pequeño de Ncr. • Finalmente, en caso de que Ncr<<1, es decir, la radiación tiende a filtrarse por un medio ópticamente delgado, la radiación se convierte en un óptimo transmisor de energía: el espesor físico es más pequeño, y puede hacer a su vez que la transmisión de calor por conducción se vea incrementada. Motivación, conclusiones y futuras líneas 78 78 Finalmente, estas serían futuras líneas de trabajo: • Programar códigos que contemplen flujo de radiación a través de medios ópticamente gruesos y delgados que bordeen figuras en 3D. • Estudiar régimen hipersónico y ionización. • Programar una condición de contorno en función del flujo de calor, en vez de la temperatura. • Convección libre con medio que absorbe y emite, del mismo modo que se ha estudiado en convección forzada en este proyecto. • Efecto de emisividades no unidad, y de dispersión, en capas límite. • Modelar de manera sinérgica las dos regiones en las que se divide un flujo viscoso, capa límite y flujo ideal, de tal forma, que una capa límite en condiciones de espesor óptico delgado se corresponda con un flujo ideal en las mismas condiciones, a través de un proceso iterativo como el de [9]. 79 Anexo: códigos Matlab %script 1: código que calcula el perfil de temperaturas para una capa %límite en placa plana, generada por un flujo viscoso que no dispersa radiación, y %confinado entre dos placas isotermas, las cuales irradian como cuerpos negros. %Dicha capa límite sí tiene capacidad de transferir calor por conducción,y %en ella es aplicable la hipótesis de espesor óptico grueso. En consecuencia, %la condición de contorno exterior de la capa límite será la temperatura de %la placa superior, enfrentada por la placa en estudio. clear all, close all,clc %parámetros Re=1e4; Pr=1; Ncr=10; Tw=0.1; Tinf=1; %discretización en xy ximin=0; ximax=1; Nxi=2000; etamin=0; etamax=10*sqrt((1+4/3/Ncr)/Re/Pr);Neta=2000; xi(1:Nxi)=linspace(ximin,ximax,Nxi); eta(1:Neta)=linspace(etamin,etamax,Neta); hxi=(ximax-ximin)/(Nxi-1); heta=(etamax-etamin)/(Neta-1); %matrices correspondientes a cocientes incrementales D1=sparse(Neta,Neta); D2=sparse(Neta,Neta); for j=2:(Neta-1), D1(j,j-1)=-1/2/heta; D1(j,j+1)=1/2/heta; D2(j,j-1)=1/heta^2; D2(j,j)=-2/heta^2; D2(j,j+1)=1/heta^2; end %primera estación capa límite mecánica unm1(1:Neta,1)=1; vnm1(1:Neta,1)=0; %primera estación capa límite térmica Thetanm1(1:Neta,1)=Tinf; %matriz C capa límite térmica C=-(1/3/Re/Pr/Ncr)*D2; C(1,:)=0; C(Neta,:)=0; for n=2:Nxi, % Capa límite mecánica en placa plana % Matriz A y vector r Matu=spdiags(unm1/hxi,0,Neta,Neta); Matv=spdiags(vnm1,0,Neta,Neta); r(1:Neta,1)=unm1.^2/hxi; A=-D2/Re+Matv*D1+Matu; % Condiciones de contorno en A y r A(1,:)=0; A(1,1)=1; r(1)=0; A(Neta,:)=0; A(Neta,Neta)=1; r(Neta)=1; % Solución un(1:Neta,1)=(A\r); % Cálculo de vn Anexo: códigos Matlab 80 80 dun(1:Neta)=(un-unm1)/hxi; vn(1,1)=0; for j=2:Neta, vn(j,1)=vn(j-1,1)-0.5*heta*(dun(j)+dun(j-1)); end % Actualizar: unm1=un; vnm1=vn; %Capa límite térmica con hipótesis de espesor óptico grueso %reasignamos Matu y Matv Matu=spdiags(un/hxi,0,Neta,Neta); Matv=spdiags(vn,0,Neta,Neta); %vector b en función de velocidades de la estación en la que nos %encontramos bt(1:Neta,1)=Matu*Thetanm1; b=[Tw; bt(2:Neta-1,1);Tinf]; %matriz B B=Matu+Matv*D1-1/Re/Pr*D2; %condiciones de contorno de B B(1,:)=0; B(1,1)=1; B(Neta,:)=0; B(Neta,Neta)=1; %método de newton raphson Tk=Thetanm1; for k=1:100, Tk3=Tk.*Tk.*Tk; Tk4=Tk.*Tk3; J=B+4*C*spdiags(Tk3,0,Neta,Neta); F=b-B*Tk-C*Tk4; dTk=J\F; Tk=Tk+dTk; if max(abs(dTk))<= 10^-12, Thetan=Tk; break end end %variable de semejanza respecto a la cual enfrentamos nuestro vector de %temperaturas para poder comparar con los resultados de Viskanta etaVisk=eta*sqrt(Re/2/xi(n)); %plot (Thetan, eta) %axis([0.4 1.1 etamin etamax]) plot(eta,Thetan) axis([etamin etamax 0 1.1]) %actualizar Thetanm1=Thetan; end %script 2: código que calcula el perfil de temperaturas para un flujo no %viscoso, que no dispersa radiación, y confinado entre dos placas isotermas, %las cuales irradian como cuerpos negros. Dicho flujo puede transferir calor %por conducción en este trabajo. Para que conste, se busca establecer una %comparación con el script 3, código basado en [2], que calcula igualmente %el perfil de temperaturas para dicho flujo, aunque en función de un aparato %matemático que permitía considerar dispersión, y emisividades de las placas %que no son unidad. A lo largo del código se muestran líneas que se %corresponden con algunas de las gráficas obtenidas para flujo ideal close all,clc,clear all %parámetros adimensionales Pr=1; Ncr=1; Tw=0.1; Tinf=1; %discretizar en x ximin=0; ximax=1; Nxi=2000;hxi=(ximax-ximin)/(Nxi-1);xi(1:Nxi)=linspace(ximin,ximax,Nxi); Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 81 %discretizar en espesor óptico Ntau=300; taumin=0; tau0=0.002; tauv(1:Ntau)=linspace(taumin,tau0,Ntau); htau=(tau0-taumin)/(Ntau-1); %integral exponencial correspondiente a exp(-X/mu); es decir, E2; %posteriormente se interpola dicho vector en tauv tauH=tau0; NX=10^3; Nmu=10^5; Xv=linspace(0,tauH,NX); mu=linspace(10^-10,1,Nmu); E2v(1)=1; E2Hv(1)=trapz(mu,exp(-tauH./mu)); for iX=2:NX, X=Xv(iX); E2v(iX)=trapz(mu,exp(-X./mu));%correspondiente a placa 1 E2Hv(iX)=trapz(mu,exp(-abs(tauH-X)./mu));%correspondiente a placa 2 end %matriz A, en función de tauv Mtau=tauv'*ones(1,Ntau)-ones(Ntau,1)*tauv; % Matriz que obtiene tau_i-tau_j, de tau-tau' Mabtau=abs(Mtau); % Matriz que obtiene |tau_i-tau_j| de |tau-tau'| ME2=interp1(Xv,E2v,Mabtau); %Matriz E2(|tau-tau'|) VE2=interp1(Xv,E2v,tauv);%vector E2(tauv), correspondiente a placa 1 VE2H=interp1(Xv,E2Hv,tauv);%vector E2(tau0-tauv), correspondiente a placa 2 Mstaup=sign(Mtau+10^-10); Mstaum=sign(Mtau-10^-10); g(1:Ntau,1:Ntau-1)=Mstaup(:,2:Ntau).*ME2(:,2:Ntau)-Mstaum(:,1:Ntau-1).*ME2(:,1:Ntau-1); Aa(1:Ntau,1)=0.5*g(:,1); Aa(:,2:Ntau-1)=0.5*(g(:,1:Ntau-2)+g(:,2:Ntau-1)); Aa(:,Ntau)=0.5*g(:,Ntau-1); %matrices correspondientes a cocientes incrementales Dtau1=sparse(Ntau,Ntau); Dtau2=sparse(Ntau,Ntau); for j=2:(Ntau-1), Dtau1(j,j-1)=-1/2/htau; Dtau1(j,j+1)=1/2/htau; Dtau2(j,j-1)=1/htau^2; Dtau2(j,j)=-2/htau^2; Dtau2(j,j+1)=1/htau^2; end %primera estación capa límite térmica Thetanm1(1:Ntau,1)=Tinf; %método de newton raphson, calculamos campo de temperaturas con radiación %en función de Ncr % contador % comp=[]; openfig('idncr01thinthickftau0002'); p=4; for n=2:Nxi xi(n) %matrices B y C con condiciones de contorno B=speye(Ntau)/hxi-2*Ncr*Dtau2;C=2*speye(Ntau)-Aa; C(1,:)=0; C(Ntau,:)=0; B(1,:)=0; B(1,1)=1; B(Ntau,:)=0; B(Ntau,Ntau)=1; %término independiente bact=Tw^4*VE2'+Thetanm1/hxi+Tinf^4*VE2H'; % b=[Tw; bact(2:Ntau-1,1) ;Tinf]; %Newton Raphson Tk=Thetanm1; for k=1:100, Tk3=Tk.*Tk.*Tk; Tk4=Tk.*Tk3; J=B+4*C*spdiags(Tk3,0,Ntau,Ntau); F=b-B*Tk-C*Tk4; dTk=J\F; Tk=Tk+dTk; Anexo: códigos Matlab 82 82 %condición para salir del método if max(abs(dTk))<= 10^-12, Thetan=Tk; break end end %obtención de matriz que contiene las componentes de temperaturas %calculadas en la malla dispuesta, con la intención de implementarlas %como condición de contorno en la capa límite de espesor óptico variable Thetanideal(1:Ntau,n)=Thetan; %perfil de temperaturas obtenido % figure(1) % plot (Thetan,tauv,'r') % axis([0 1.1 taumin tau0]) % pause(0.01) if n==100|n==1000|n==2000 figure(1) subplot(2,3,p) plot (Thetan,tauv,'r') axis([0 1.1 taumin tau0]) % legend(['ideal f(\tau)'],['ideal thick']); title(['n=',num2str(n)]); xlabel('\theta');ylabel('\tau'); filename=['idncr01thinthickftau0002'];saveas([1],filename); p=p+1; end if max(abs((Thetan-Thetanm1)/hxi)) <= 10^-7%||n==Nxi, % comp(2,1:4)=[Thetan(50),Thetan(70),Thetan(100),Thetan(150),Thetan(290)]; % xlswrite('compideal2.xlsx',comp,'Hoja1','B2'); % %condición que se emplea para llegar a la estacionariedad en % %términos de posición longitudinal a lo largo de la placa plana; % %empleada para comparar la curva estacionaria a la que llega el % %flujo ideal con hipótesis de espesor óptico grueso % [max(abs((Thetan-Thetanm1)/hxi))] figure(1) subplot(2,3,p) plot (Thetan,tauv,'r') axis([0 1.1 taumin tau0]) % legend(['ideal f(\tau)'],['ideal thick']); title(['n=',num2str(n)]); xlabel('\theta');ylabel('\tau'); filename=['idncr01thinthickftau0002'];saveas([1],filename); % filename=['pruebaidealthickk'];saveas([1],filename); p=p+1; % % % break; % else % [max(abs((Thetan-Thetanm1)/hxi))] % hold off end Thetanm1=Thetan; % figure(1) % com theta(n)=(Thetan(2)-Tw)/(Tinf-Tw); % semilogx(xi(2:n),theta(2:n),'g') % hold on end Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 83 % openfig('sparrow289') % hold on % semilogx(xi(2:Nxi),theta(2:Nxi),'k') % hold on % axis([0.01 4 0 1]); % xlabel('\xi') % ylabel('\Theta') % title('Salto de temperatura en la pared en flujo no viscoso. Pared caliente y pared fría'); % legend(['\theta_w=1, \theta_0=0.1'],['\theta_w=0.1, \theta_0=1'],['\theta_w=1, \theta_0=0.9'],['\theta_w=0.9, \theta_0=1']); % filename=['sparrow289']; % saveas([1],filename); %script 3: flujo ideal encerrado entre dos placas que permite transmisión% %de calor, por conducción y por radiación, por diferencia de temperaturas% %entre dichas placas; cálculo a través de matrices A y MG clear all; clc; close all; %parámetros adimensionales emiss1=1; emiss2=1; Ncr=0;Tw=1;Tinf=0.1; %discretizar en x ximin=0; ximax=1200; Nxi=2000;hxi=(ximax-ximin)/(Nxi-1); %discretizar en tau tau0=20; Ntau=300; taumin=0; tau0; tauv(1:Ntau)=linspace(taumin,tau0,Ntau); htau=(tau0-taumin)/(Ntau-1); %integral exponencial correspondiendo a exp(-X/mu) y exp(-abs(tauH-X)); %calculada con una malla muy fina; cuando sea necesario se interpolará tauH=tau0; NX=10^3; Nmu=10^5; Xv=linspace(0,1.5*tauH,NX); mu=linspace(10^-10,1,Nmu); E2v(1)=1; E2Hv(1)=trapz(mu,exp(-tauH./mu)); for iX=2:NX, X=Xv(iX); E2v(iX)=trapz(mu,exp(-X./mu)); E2Hv(iX)=trapz(mu,exp(-abs(tauH-X)./mu)); end E3_H=trapz(mu,mu.*exp(-tauH./mu)); %para radiosidades %Matriz A Mtau=tauv'*ones(1,Ntau)-ones(Ntau,1)*tauv; % Matriz de la discretización tau_i-tau_j de tau-tau' Mabtau=abs(Mtau); % Matriz de la discretización |tau_i-tau_j| de |tau-tau'| ME2=interp1(Xv,E2v,Mabtau); % Matriz de la discretizacíon E2(|tau-tau'|) Mstaup=sign(Mtau+10^-10); Mstaum=sign(Mtau-10^-10); g(1:Ntau,1:Ntau-1)=Mstaup(:,2:Ntau).*ME2(:,2:Ntau)-Mstaum(:,1:Ntau-1).*ME2(:,1:Ntau-1); % matrix g A(1:Ntau,1)=0.5*g(:,1); A(:,2:Ntau-1)=0.5*(g(:,1:Ntau-2)+g(:,2:Ntau-1)); A(:,Ntau)=0.5*g(:,Ntau- 1); %Vectores a y b: vE2(1:Ntau)=interp1(Xv,E2v,tauv); E2(1)=1; vE2H(1:Ntau)=interp1(Xv,E2Hv,tauv); E2H(Ntau)=1; a(1:Ntau)=0.5*[vE2(1:Ntau-1).*(tauv(2:Ntau)-tauv(1:Ntau-1)) 0]+0.5*[0 vE2(2:Ntau).*(tauv(2:Ntau)-tauv(1:Ntau-1))]; b(1:Ntau)=0.5*[vE2H(1:Ntau-1).*(tauv(2:Ntau)-tauv(1:Ntau-1)) 0]+0.5*[0 vE2H(2:Ntau).*(tauv(2:Ntau)-tauv(1:Ntau-1))]; %Matrices LGJ and RIb: LGJ(1:Ntau,1:Ntau)=eye(Ntau); LGJ(1:Ntau,Ntau+1)=-2*vE2(1:Ntau)'; LGJ(1:Ntau,Ntau+2)=- 2*vE2H(1:Ntau)'; Anexo: códigos Matlab 84 84 LGJ(Ntau+1,1:Ntau)=0; LGJ(Ntau+1,Ntau+1)=1; LGJ(Ntau+1,Ntau+2)=-2*(1-emiss1)*E3_H; LGJ(Ntau+2,1:Ntau)=0; LGJ(Ntau+2,Ntau+1)=0; LGJ(Ntau+2,Ntau+2)=1; %J2=0 RIb(1:Ntau,1:Ntau)=2*pi*A; RIb(Ntau+1,1:Ntau)=pi*emiss1*[1 zeros(1,Ntau-1)]+2*pi*(1-emiss1)*a; %RIb(Ntau+2,1:Ntau)=2*pi*(1-emiss2)*(1-omega)*b+pi*emiss2*[zeros(1,Ntau-1) 1]; RIb(Ntau+2,1:Ntau)=0; %J2=0% LGJi=inv(LGJ); LGJiRIb=LGJi*RIb; MG=LGJiRIb(1:Ntau,1:Ntau); aJ1=LGJiRIb(Ntau+1,1:Ntau); aJ2=LGJiRIb(Ntau+2,1:Ntau); %matrices correspondientes a cocientes incrementales D_tau=sparse(Ntau,Ntau); D_tau2=sparse(Ntau,Ntau); for j=2:(Ntau-1), D_tau(j,j-1)=-1/2/htau; D_tau(j,j+1)=1/2/htau; D_tau2(j,j-1)=1/htau^2; D_tau2(j,j)=-2/htau^2; D_tau2(j,j+1)=1/htau^2; end %condición de contorno en temperatura Tnm1(1:Ntau,1)=Tinf; comp=[]; for n=2:1000000 b=[Tw; Tnm1(2:Ntau-1,1); Tinf]; % Matrices del sistema y condiciones de contorno B=eye(Ntau)-hxi*Ncr*D_tau2; B(Ntau,:)=0; B(Ntau,Ntau)=1; B(Ntau,Ntau-1)=-1; C=hxi*(eye(Ntau,Ntau)-MG/4/pi); B(1,:)=0; B(1,1)=1; C(1,:)=0; B(Ntau,:)=0; B(Ntau,Ntau)=1; C(Ntau,:)=0; %Iteración Newton Raphson Tk=Tnm1; for k=1:100, Tk3=Tk.*Tk.*Tk; Tk4=Tk3.*Tk; Fk=B*Tk+C*Tk4-b; Hk=B+4*C*spdiags(Tk3,0,Ntau,Ntau); dTk=-Hk\Fk; if max(abs(dTk)) <= 10^-12, break else Tk=Tk+dTk; end end Tn=Tk; Tn4=Tk4; [max(abs((Tn-Tnm1)/hxi))] if max(abs((Tn-Tnm1)/hxi)) <= 10^-7, openfig('estacionario'); hold on plot(Tn,tauv,'b'); title('flujo uniforme en placa plana'); xlabel('\theta');ylabel('\tau'); legend(['código original'],['código [2]']); comp(2,1:5)=[Tn(50),Tn(70),Tn(100),Tn(150),Tn(290)]; xlswrite('compideal.xlsx',comp,'Hoja1','B2'); break; end Métodos numéricos sencillos para problemas de radiación en regímenes de convección forzada y libre 85 Tnm1=Tn; % if 5*round(n/5)==n, % figure(1) % plot(Tn,tauv) % hold off %figure(2) %plot(n,theta,'b') %hold on %pause(0.5) %end end %script4: capa límite viscosa sobre una placa plana, en presencia de radiación, %generada por un flujo ideal que no conduce calor y confinado entre dos placas %que se comportan como cuerpos negros. Dicha capa límite queda en función %del espesor óptico tauCL, de la condición de contorno correspondiente a %dicha tauCL en el flujo ideal, cargado en memoria, que se calcula con %script2, y del flujo de radiación que tiene en cuenta absorción y emisión %en el medio, relajándose hipótesis thin y thick clear all, close all,clc %parámetros adimensionales Pr=1; Ncr=0.1; Tw=0.1; Tinf=1;tau0=20; %discretizar en x ximin=0; ximax=1; Nxi=2000;hxi=(ximax-ximin)/(Nxi-1); xi(1:Nxi)=linspace(ximin,ximax,Nxi); %discretizar en espesor óptico Ntau=1000; taumin=0; tauCL=10*sqrt(2*Pr*Ncr); tauv(1:Ntau)=linspace(taumin,tauCL,Ntau); htau=(tauCL-taumin)/(Ntau-1);heta=htau/sqrt(2*Pr*Ncr); eta(1:Ntau)=tauv/sqrt(2*Pr*Ncr);%eta=tauadi de sparrow %discretizar en espesor óptico de flujo ideal tauvideal(1:Ntau)=linspace(taumin,tau0,Ntau); %variable diferencia para saber en qué posición de tauvideal se encuentra %la condición de contorno exterior de la capa límite dif=abs(tauvideal(1)-tauCL); p=1; for s=2:Ntau if abs(tauvideal(s)-tauCL)<dif dif=abs(tauvideal(s)-tauCL); p=s end end %p es la posición de tauvideal en la que se encuentra la condición de %contorno exterior de la capa límite pause %se carga el flujo ideal correspondiente a tau0=20, y pared fría load ('thetaccideal20fria','-mat'); CC=Thetanideal; %obtención de E2(tau), E2(|tau0-tau|), y Aij, siempre respecto a tauCL tauH=tauCL; NX=10^3; Nmu=10^5; Xv=linspace(0,tauH,NX); mu=linspace(10^-10,1,Nmu); E2v(1)=1; E2Hv(1)=trapz(mu,exp(-tauH./mu)); for iX=2:NX, Referencias 92 92 [17] M. Á. M. Medina, «Gaussianos,» 22 Junio 2011. [En línea]. Available: https://www.gaussianos.com/calcular-la-derivada-de-una-integral/. [Último acceso: 22 Mayo 2020]. [18] M. P.-S. Sánchez-Pastor, «A numerical method for transient, one-dimensional conduction-radiation- convection problems with applications to the thermal structure of simple models of planetary lower atmospheres,» Sevilla, 2020. [19] G. D. D. Reina, «Métodos numéricos sencillos para el cálculo de chorros y penachos térmicos laminares,» Sevilla, 2016. [20] S. Rosseland, Astrophysik auf Atom-Theoretischer Grundlage, Berlin, 1931.