Estudio numérico de la convección de Rayleigh-Benard en coronas circulares: aspectos fluido-mecánicos y energéticos

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Estudio numérico de la convección de Rayleigh-Benard en coronas circulares: aspectos fluido-mecánicos y energéticos Autor: Alejandro Guzmán Pozo Tutor: Miguel Perez-Saborid Sánchez-Pastor Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Estudio numérico de la convección de Rayleigh-Benard en coronas circulares: aspectos fluido-mecánicos y energéticos Autor: Alejandro Guzmán Pozo Tutor: Miguel Perez-Saborid Sánchez-Pastor Profesor Titular Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Grado: Estudio numérico de la convección de Rayleigh-Benard en coronas circulares: aspectos fluido-mecánicos y energéticos Autor: Alejandro Guzmán Pozo Tutor: Miguel Perez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos A mi abuela, mi madre, mi padre, mi hermana y Fabio. Gracias por apoyarme en todo momento. Gracias a ti también, Linda. A mi tutor, Miguel. Gracias por acompañarme y guiarme en todas las dificultades que han surgido en los meses de elaboración del proyecto. Alejandro Sevilla, 2024 I Resumen El objetivo de este trabajo es estudiar el fluido en el interior de una corona circular con gravedad radial cuya pared interior está a mayor temperatura que la pared exterior. De esta forma, el fluido está sometido a un gradiente de temperatura que provocará variaciones de densidad. Estas variaciones junto al efecto de la gravedad desembocan eventualmente en el conocido fenómeno de la convección de Rayleigh-Bénard. Para el análisis de la convección de Rayleigh-Bénard en una corona circular, se formula el problema matemáticamente partiendo de las ecuaciones general de Navier-Stokes y, tras una serie de simplificaciones e hipótesis, se obtienen las ecuaciones de Saltzman en coordenadas cilíndricas bidimensionales. Estas ecuaciones definen el problema de la convección de Rayleigh-Bénard en una corona circular. En cuanto al método numérico empleado para resolver las ecuaciones de Saltzman, se presenta el método de colocación mediante interpolación de Lagrange. Empleando este método de resolución numérica, se obtienen las soluciones de las ecuaciones de Poisson y de la Laplace en régimen no estacionario en una geometría anular. Posteriormente, se resuelven las ecuaciones de Saltzman y se analizan los resultados obtenidos. Para verificar que estos resultados son coherentes, se comparan con otros resultados clásicos encontrados en la literatura. Cabe destacar que no se emplea el desarrollo de Fourier para problemas periódicos, normalmente empleados en la literatura, ya que el problema de la convección de Rayleigh-Bénard es altamente no lineal (debido al término convectivo) y la proyección en las funciones de Fourier daría lugar a términos de acoplamiento complicados. Finalmente, se propone una serie de líneas de futuras extensiones, complementarias a los análisis y resultados expuestos en el presente trabajo. III XÍndice 5 Conclusiones y líneas futuras 79 5.1 Conclusiones 79 5.2 Líneas futuras 80 Índice de Figuras 83 Índice de Tablas 85 Índice de Códigos 87 Bibliografía 89 1 Introducción 1.1 Motivación y objetivos del proyecto La ingeniería siempre ha tratado de facilitar al ser humano el resolver problemas a partir de los recursos disponibles. Uno de los ejemplos más antiguos de los que se tiene registro es la invención de la fundición de metales como el cobre y el bronce que permitió la creación de herramientas y armas de metal, que hasta ese momento se hacían de materiales menos resistentes como la madera. Por suerte, la cantidad de recursos disponibles ha aumentado a lo largo de la historia. Sin embargo, siguen suponiendo una limitación que obliga al ser humano a pensar de manera más creativa y eficiente dando cabida al progreso. En el presente trabajo se llevará a cabo un exhaustivo análisis numérico de la convección de Rayleigh-Bénard dentro de un recinto en forma de corona circular con gravedad radial. Este fenómeno tiene un papel fundamental en numerosos procesos industriales y naturales. En efecto, su entendimiento resulta imprescindible a la hora de analizar los procesos que tienen lugar en sistemas de enfriamiento y calefacción (desde el punto de vista industrial), y en el manto terrestre, en atmósferas planetarias y en el interior de las estrellas (desde el punto de vista de procesos naturales). La resolución de problemas dentro de los campos de la Mecánica de Fluidos y Transferencia de Calor, campos donde se enmarca la convección de Rayleigh-Bénard, siempre ha supuesto un reto complejo ya que la formulación matemática de estos problemas resulta en sistemas de ecuaciones diferenciales en derivadas parciales que a menudo no tienen solución analítica. Tanto es así que, tradicionalmente, los estudiantes de las escuelas de ingeniería perdían la perspectiva de los conceptos fundamentales físicos de los problemas al someterse a largas horas de manipulaciones matemáticas tratando de resolver problemas que habían sido previamente simplificados con el objetivo de forzar la existencia de solución analítica. Actualmente, existen numerosas herramientas que permiten la resolución numérica de problemas ingenieriles. Así, en el caso del presente trabajo, se partirá de las ecuaciones de Navier-Stokes, que describen los movimientos de los fluidos, y bajo una serie de simplificaciones e hipótesis se obtendrá un sistema de ecuaciones diferenciales de cuarto orden en dos variables que se resolverá numéricamente mediante el programa MATLAB. De esta forma, se elaborará una rutina que permita implementar estas ecuaciones en MATLAB y obtener resultados relavantes de ellas que posibilite entender los aspectos fluido-mecánicos y energéticos de la convección de Rayleigh-Bénard en una corona circular. Previamente, se abordarán problemas 1 2Capítulo 1. Introducción menos complejos desde el punto de vista matemático, la ecuación de Poisson y la ecuación de Laplace en régimen no estacionario, que permitirán elaborar y desarrollar un método numérico de resolución posteriormente empleado para la convección de Rayleigh-Bénard. 1.2 Problema de convección natural de Rayleigh-Bénard Se denomina fenómeno de convección natural a los movimientos producidos en el seno de un fluido que está sometido a un campo gravitatorio y que presenta diferencias de densidad entre sus partes provocadas, en la mayoría de las ocasiones, por gradientes de temperatura o de concentración. Sin embargo, el fenómeno de convección natural también se da en fluidos con carga eléctrica cuando están sometidos a la acción de campos eléctricos o magnéticos. De todas las convecciones naturales posibles, este proyecto se centrará en el caso de un fluido sometido a un campo gravitatorio y en el que aparecen variaciones de densidad provocadas por gradientes de temperatura. Así, los movimientos que aparecen en el fluido bajo estas condiciones se denominan fenómeno de convección natural de Rayleigh-Bénard. La convección de Rayleigh-Bénard es un tipo de convección natural que se produce en una capa plana y horizontal de fluido que es calentado desde abajo. Como resultado, el fluido desarrolla un patrón regular de celdas de convección conocidas como celdas de Bénard. Figura 1.1 Celdas originadas debido a la convección de Rayleigh-Bénard. Este tipo de convección se manifiesta en numerosos ejemplos de la naturaleza. Así, en el manto terrestre se producen corrientes de convección natural debido a los gradientes de temperaturas existente entre el núcleo y la corteza de la Tierra. Por otro lado, en la atmósfera terrestre también se da el fenómeno de Rayleigh-Bénard, formándose distintas células de convección según la latitud: célula de Hadley (en el entorno del ecuador y responsable de los vientos alisios), célula de Ferrell (latitudes medias) y célula Polar. 1.2 Problema de convección natural de Rayleigh-Bénard 3 ((a)) Corrientes de convección en el manto terrestre. ((b)) Células de convección en la atomósfera terrestre. Figura 1.2 Ejemplos de convección de Rayleigh-Bénard en la naturaleza. Históricamente, en torno al año 1900, Henri Bénard fue el responsable de llevar a cabo los primeros estudios sistemáticos de la convección natural. Uno de ellos consistió en calentar desde abajo una capa delgada de esperma de ballena, que posee alta conductividad térmica, cuya superficie superior estaba expuesta al aire del ambiente. Mediante este experimento, fue capaz de obervar, tras alcanzarse una determinada temperatura, la aparición de movimientos en el fluido que, después de un régimen transitorio inicial, alcanzaba un régimen estacionario en el que eran apreciables las denominadas celdas de Bénard en la superficie superior del mismo. Geométricamente, tenían forma de hexágono regular, siendo similares a las celdas de los paneles de abeja, como se puede observar en la figura (1.1). El experimento, entre otros resultados, concluyó que el fluido ascendía por el interior de la celda y descendía a lo largo del perímetro de la misma. De esta forma, se explicaba la formación de estas celdas convectivas. Además, Bénard relacionó la visualmente perceptible depresión que aparecía en la superficie del fluido en contacto con el aire por la zona del centro de la célula con la tensión superficial existente en dicha superficie. Por otro lado, también analizó la dependencia entre el tamaño característico de las células con la profundidad de la capa de fluido, el flujo de calor y la temperatura. Figura 1.3 Patrón de movimiento de convección desarrollado en una capa de fluido al calentarla desde abajo. Figura 1.4 Vista lateral de las celdas de Bénard. 4Capítulo 1. Introducción Pasados unos años, en 1916, Lord Rayleigh fue el primer investigador que logró poner de manifiesto los mecanismos físicos involucrados en los fenómenos de convección natural mediante la publicación del artículo On Convection Currents in a Horizontal Layer of Fluid, When the Higher Temperature Is Under Side [ 9 ]. Así, cuando un fluido con densidad no uniforme se encuentra en presencia de la gravedad, la fuerza de flotabilidad (también conocida como empuje de Arquímedes) provoca que la porción de fluido que, en media, es más densa que su entorno tienda a descender y que, consecuentemente, la de menor densidad tienda a ascender. Además de la fuerza de flotabilidad, y como destacó Rayleigh, también se deben de considerar otros mecanismos que se encargan de contrarrestar su efecto como la fricción (originada por las fuerzas de viscosidad) y la conducción de calor (tiende a homogeneizar el campo de temperaturas y de densidades del fluido). De esta forma, al aplicarse una perturbación en un fluido inicialmente en reposo con una distribución de densidades creciente, se pueden producir dos situaciones diferentes en función de la importancia de los efectos de flotabilidad, fricción y conducción de calor: • Situación de equilibrio estable, en la que el fluido eventualmente alcanzará cierta estabilidad regresando a la posición de reposo inicial. • Situación de equilibrio inestable, en la que se inducirá un movimiento ascendente de carácter permanente en el fluido que originará un flujo convectivo responsable del fenómeno de convección natural. Rayleigh consideró una capa de fluido confinada entre dos placas paralelas e infinitas, dispuestas perpendicularmente a la dirección de la gravedad y a distinta temperatura (siendo la placa inferior la de mayor temperatura). En su análisis, linealizó las ecuaciones de Navier- Stokes, que gobiernan el movimiento del fluido entre las placas, en torno a la solución de equilibrio puesto que estaba interesado en analizar la estabilidad de pequeñas perturbaciones respecto a la situación de reposo. Como resultado, obtuvo un sistema lineal de ecuaciones en derivadas parciales en el que aparece un parámetro adimensional denominado número de Rayleigh: Ra =gβ∆TH3 να ,(1.1) donde: •βes el coficiente de dilatación térmica. •νes la viscosidad cinemática. •αes la difusividad térmica del fluido. •ges la aceleración de la gravedad. •H es el espesor de la capa fluida (o, en general, una longitud característica del dominio fluido). •∆Tes la diferencia de temperatura entre placas. El número de Rayleigh puede interpretarse físicamente como un parámetro adimensional que mide la importancia relativa entre los efectos de las fuerzas de flotabilidad y los efectos de las fuerzas de viscosidad y de conducción térmica. En efecto, existe un valor crítico de este parámetro adimensional, denominado Rac , a partir del cual el sistema se vuelve inestable y los movimientos de convección natural se hacen presentes a través de las denominadas celdas convecticas. De esta forma, el valor crítico de Rayleigh define el punto a partir del cual una 1.3 Estructura del proyecto 5 pequeña perturbación en el sistema desemboca en un régimen de convección natural en el fluido. A lo largo de la historia, una gran cantidad de autores han estudiado el valor del Rayleigh crítico para sistemas particulares. Entre ellos destaca Jeffreys (1926), quien obtuvo por primera vez el número de Rayleigh crítico para un sistema compuesto por placas infinitas ( Rac=1708 ). Por otro lado, Davis (1967) y Catton (1970), entre otros, proporcionaron resultados teóricos y experimentales de valores de Rac para sistemas con diferentes geometrías y conductividades en las paredes laterales. En secciones posteriores, se analizarán tanto la existencia de un valor de Rayleigh crítico como la existencia de un régimen turbulento de convección natural para números de Rayleigh muy grandes (en torno a 105), es decir, para Ra >> Rac. Además de los trabajos de Bénard y de Rayleigh, numerosos han sido los estudios que se han llevado a cabo relacionados con la convección natural en fluidos confinados. Se destaca el trabajo de Batchelor (1954), quien abarcó el problema de la convección natural en una cavidad rectangular de longitud L y altura H, analizándolo a partir de tres parámetros adimensionales: la relación de aspecto, A=L H , el número de Prandtl, Pr =ν α y el número de Rayleigh, cuya expresión es la indicada en la ecuación (1.1). El presente trabajo se centrará en el estudio bidimensional de la convección natural en una cavidad en forma de corona circular, con radio interior Ri y radio exterior Re , en la que existe una diferencia de temperatura entre sus paredes. En efecto, se produce un gradiente de temperaturas entre las paredes exterior e interior de la corona, siendo esta última la de mayor temperatura. Por tanto, el objetivo del trabajo será analizar el fenómeno de convección natural originado por variaciones de densidad asociadas a gradientes de temperatura en fluidos de composición homogénea. Sin embargo, existen también numerosos casos en los que la convección natural se debe a variaciones de densidad asociadas a gradientes de concentración en fluidos de composición heterogénea. Finalmente, de manera análoga al trabajo de Batchelor, los parámetros adimensionales que gobernarán el problema serán: • Relación de aspecto, λ=Re Ri<1 , que caracteriza las propiedades geométricas de la cavidad. •Número de Prandtl, Pr =ν α, que caracteriza las propiedades físicas del fluido. •Número de Rayleigh, Ra =gH3β∆T να . 1.3 Estructura del proyecto En este trabajo se abordará un estudio numérico y teórico del fenómeno de convección natural de Rayleigh-Bénard en una cavidad en forma de corona circular calentada desde la pared interior y con gravedad radial. Los análisis que se llevarán a cabo siempre estarán dentro del régimen laminar. Así, el movimiento del fluido que resulte de la convección será siempre ordenado y no turbulento. 6Capítulo 1. Introducción En el capítulo 2, partiendo de las ecuaciones generales de Navier-Stokes y asumiendo la simplificación del problema basada en la aproximación de Boussinesq, se obtendrán las ecuaciones de la convección de Rayleigh-Bénard en coordenadas cilíndricas bidimensionales, que forman un sistema de ecuaciones diferenciales de cuatro ecuaciones con cuatro incógnitas. Seguidamente, se transformará el sistema de ecuaciones obtenido en uno formado por dos ecuaciones con dos incógnitas: las denominadas ecuaciones de Saltzman. Para ello, se introducirá el concepto de función de corriente Ψ . Por último, se adimensionalizarán las ecuaciones de Saltzman mediante la introducción de una serie de variables adimensionales que definirán el problema. Esto facilitará la posterior resolución numérica del problema de convencción de Rayleigh-Bénard en MATLAB. Además, conllevarla la aparición de los parámetros adimensionales Ra y Pr en las ecuaciones del problema que permitirán modificar el régimen del fluido alojado en la cavidad bajo estudio. En el capítulo 3, se desarrollará el método numérico que se empleará para resolver el sistema de ecuaciones diferenciales en derivadas parciales que constituye las ecuaciones de Saltzman: el denominado método de colocación mediante interpolación de Lagrange. Para ello, se desarrollará el método para problemas unidimensionales y se presentará un procedimiento para cálcular matrices de diferenciación cuyo autor es el tutor del presente trabajo, Miguel Pérez- Saborid. Seguidamente, se extenderá el método a problemas bidimensionales y se ilustrará su aplicación para la resolución de la ecuación de Poisson y la ecuación de Laplace en régimen no estacionario mediante MATLAB. Finalmente, se detallará la implementación numérica de las ecuaciones de Saltzman adimensionales en MATLAB y se obtendrá una solución para el problema. En el capítulo 4, se presentarán los resultados de las ecuaciones de la convección natural de Rayleigh-Bénard en una corona circular para diferentes configuraciones definidas por parámetros como la relación de aspecto λ y el número de Rayleigh Ra . Para ello, se analizarán diferentes regímenes del fluido en función de parámetros adimensionales y se compararán los resultados obtenidos con información encontrada en la literatura. Así, se estudiará la evolución del fluido dentro de la cavidad para diferentes valores de Ra y una geometría fija. Esto permitirá definir el concepto de número de Rayleigh crítico Rac . Posterioremente, se evaluará la relación entre el número de Rayleigh crítico y la geometría de la cavidad. Por último, el cálculo del número de Nusselt ayudará a comprender el tipo de transferencia de calor que tiene lugar en el fluido dentro de la corona circular en cada instante de su evolución temporal en función de número de Rayleigh para una geometría fija. 2 Formulación del problema. Ecuaciones de Saltzman en coordenadas cilíndricas 2D 2.1 Ecuaciones de Navier-Stokes en coordendas cilíndricas bidimensionales Considérese un fluido en el interior de una cavidad en forma de corona circular de radio interior Ri y radio exterior Re , afectada, radialmente, por el campo gravitatorio terrestre ( g=−ger ) y por una diferencia de temperatura entre sus paredes interior ( T1 ) y exterior ( T2 ), siendo T1>T2 . Dicho sistema se estudiará bidimensionalmente según las coordenadas r y θ , como se muestra en la figura (2.1). El flujo vendrá caracterizado por un campo de velocidades v=vrer+vθeθ , densidades ρ(t,r,θ) presiones p(t,r,θ) y temperaturas T(t,r,θ) . Para el análisis del problema se emplearán las ecuaciones de Navier-Stokes en coordenadas cilíndricas bidimensinales (r,θ) . Este sistema, de forma genérica, está formado por las siguientes ecuaciones de conservación de la masa, de la cantidad de movimiento y de la energía, respectivamente: ∂ρ ∂t+∇·(ρv) = 0,(2.1) ρDv Dt =−∇p+∇·¯ ¯ τ′+ρg,(2.2) ρDe Dt =−p∇·v+¯ ¯ τ′:∇v+∇·(k∇T)+Qr.(2.3) En las ecuaciones (2.1)-(2.2), las incógnitas serán la densidad ρ , las dos componentes de la velocidad v ( vr y vθ ), la presión p , la energía interna e y la temperatura T . Así, para completar dicho sistema serán necesarias dos ecuaciones de estado, como por ejemplo: ρ=ρ(p,T)ye=e(T,p). 7 8Capítulo 2. Formulación del problema. Ecuaciones de Saltzman en coordenadas cilíndricas 2D Figura 2.1 Corona circular en coordenadas ryθ. De este modo, se tendrá un sistema compatible determinado que proporcionará el valor de las incógnitas bajo unas condiciones determinadas. 2.1.1 Simplificación del problema: aproximación de Boussinesq En una primera aproximación, las ecuaciones que gobiernan el movimiento del fluido dentro de una determinada cavidad son las ecuaciones de Navier-Stokes completas, correspondientes a un fluido compresible, viscoso y conductor de calor. Sin embargo, como la resolución de estas ecuaciones es muy compleja, debido principalmente a su carácter no lineal, es muy conveniente tratar de buscar simplificaciones que ayuden a facilitar el estudio del problema. En este contexto surge la denominada aproximación de Boussinesq. Esta aproximación, que constituye una simplificación bastante significativa, es aplicable siempre y cuando las variaciones relativas de densidad involucradas en el proceso de convección sean pequeñas. Es importante destacar que, para que esto ocurra, es necesario que las variaciones de la presión y temperatura sean también pequeñas en comparación con un estado de referencia definido por p0yT0(que son magnitudes características del problema). De este modo, para analizar cómo son las variaciones relativas de densidad, en primer lugar, se llevará a cabo el desarrollo en serie de Taylor de primer orden de la ecuación de estado del fluido ρ=ρ(p,T): ρ=ρ(p0,T0)+ ∂ρ ∂T0 (T−T0)+ ∂ρ ∂p0 (p−p0). Dividiendo entre ρ0 y despejando la expresión anterior se obtiene la variación relativa de la densidad con respecto al estado de referencia: ρ−ρ0 ρ0 =1 ρ ∂ρ ∂T0 (T−T0)+ 1 ρ ∂ρ ∂p0 (p−p0).(2.4) 2.1 Ecuaciones de Navier-Stokes en coordendas cilíndricas bidimensionales 9 Utilizando las definiciones de coeficientes de expansión térmica, β , y de compresibilidad, k , la expresión (2.4) puede escribirse como: ρ−ρ0 ρ0 =−β(T−T0)+k(p−p0).(2.5) En el caso de la convección de Rayleigh-Bénard, en la que las variaciones de temperatura desempeñan un papel fundamental en el movimiento del fluido, se puede demostrar que la influencia de las variaciones de presión en las variaciones de densidad es insignificante y, por tanto, despreciable. Así, a partir de la ecuación de cantidad de movimiento según la dirección r (perteneciente al sistema de ecuaciones de Navier-Stokes) y la ecuación de Bernoulli, se obtiene que los órdenes de magnitud de las velocidades y variaciones de presión características del movimiento son, respectivamente: v0∼pgβ(T−T0)(Re−Ri), p−p0∼ρ0v2 0∼ρ0gβ(T−T0)(Re−Ri). De este modo, la relación entre los dos términos del miembro derecho de la expresión (2.5) es la siguiente: k(p−p0) β(T−T0)∼kρ0g(Re−Ri).(2.6) Dada una cavidad con (Re−Ri)∼2m y una aproximación de los órdenes de magnitud de ρ , β y k para un líquido (agua) y un gas (aire), en los que (T−T0)∼30 K , empleando la expresión (2.6) se concluye que las variaciones de presión, según el fluido, tienen los siguientes órdenes de magnitud: Agua ⇒∆p∼60 N/m2<< 105N/m2(1atm). Aire ⇒∆p∼2N/m2<< 105N/m2(1atm). Tabla 2.1 Aproximación de valores (ρ,β,k) para agua y aire. Fluido ρ[kg/m3]β[K−1]k[N/m2]−1 Agua (líquido) 1000 1,5·10−40,5·10−9 Aire (gas) 1,2 1/300 10−5 Al analizar estos valores se deduce que las variaciones de presión son pequeñas y, por tanto, despreciables. Además, se puede demostrar que el término kρ0g(Re−Ri) es despreciable para valores de (Re−Ri)<< 11 m en gases y (Re−Ri)<< 100 m en líquidos. Consecuentemente, en la mayoría de las ocasiones, la ecuación de estado (2.5) puede aproximarse por: ρ−ρ0 ρ0≈−β(T−T0).(2.7) Una vez conseguida una aproximación de las variaciones relativas de densidad únicamente dependendiente de las variaciones de temperatura (2.7), para hacer que las variaciones relativas de densidad sean pequeñas bastará con hacer que el término de diferencia de temperatura sea pequeño también. En primer lugar, para el caso del aire (gas), que tiene β∼1/300 K−1 , 3 Método numérico de resolución. El método de colocación La resolución del problema de convección natural de Rayleigh-Bénard consiste en el cálculo de las incógnitas de un sistema de ecuaciones diferenciales en derivadas parciales. Por tanto, será necesario aplicar un método numérico que permita obtener la solución de dichas ecuaciones. En este contexto, se presenta el denominado método de colocación. El método de colocación consiste en expresar las funciones incógnitas de un problema como combinaciones lineales de funciones conocidas, como, por ejemplo, los interpolantes de Lagrange. El cálculo de los coeficientes de estas combinaciones lineales resulantes se hará sustituyendo estas expansiones, así como sus derivadas, en las ecuaciones diferenciales y en las condiciones de contorno que constituyen el problema, imponiendo, además, que estas se verifiquen en un determinado conjunto de puntos del dominio de integración, denominados puntos de colocación. Así, los coeficientes quedarán completamente determinados y la resolución del problema será posible. Es importante destacar que no se ha empleado el desarrollo de Fourier para problemas periódicos, normalmente empleados en la literatura, ya que el problema de la convección natural de Rayleigh-Bénard es altamente no lineal (debido al término convectivo) y la proyección en las funciones de Fourier daría lugar a términos de acomplamiento complicados. De esta forma, para la resolución numérica de las ecuaciones de Saltzman, se ha optado por emplear un método de colocación convencional en la dirección θ del problema e imponer condiciones de periodicidad haciendo que en θ=0 y θ=2π se tenga continuidad de la función de corriente Ψ∗ hasta las derivadas de orden 3 (debido a que el problema para la función de corriente es de orden 4 con respecto a θ∗ ) y hasta las derivadas de primer orden para la temperatura ϕ∗ (debido a que el problema para la temperatura es de orden 2 con respecto a θ∗). 3.1 Método de colocación aplicado a problemas unidimensionales Sea un problema diferencial en una sola variable independiente que se desea resolver en un determinado intervalo [a,b], siendo f1,f2,...,fN los valores que toma la función incógnita en N puntos del intervalo a=x1<x2<x3< ... < xN−1<xN=b . El objetivo es tratar de obtener una aproximación de la solución del problema, denominada f(x) , en el intervalo [a,b] mediante una expansión de la forma 17 18 Capítulo 3. Método numérico de resolución. El método de colocación f(x) = N ∑ i=1 Li(x)fi,(3.1) donde las funciones Li(x) serán elegidas como polinomios interpolantes de Lagrange para cada punto del dominio, los cuales vienen definidos por la siguiente expresión: Li(x) = N ∏ j=1 j=i x−xj xi−xj i=1,...,N.(3.2) De la definición (3.2) se deduce que Li(xj) = δji ∀i,j=1,..,N , siendo δi j la delta de Kronecker, definida de la siguiente forma: δi j =1,i=j. 0,i=j.(3.3) Además, la definición (3.2) permite obtener directamente la expresión de la derivada de Li(x): L′ i(x) = N ∑ k=1 (x−x1)...(x−xi−1)(x−xi+1)...(x−xk−1)(x−xk+1)...(x−xN) (xi−x1)...(xi−xi−1)(xi−xi+1)...(xi−xN).(3.4) De este modo, la expresión anterior permite obtener el valor de la derivada de f(x) en cualquier nodo del intervalo como f′(xi) = N ∑ j=1 dLj(x) dx x=xi fj= N ∑ j=1 L′j(xi)fji=1,...,N.(3.5) Análogamente, también se puede obtener el valor de la segunda derivada de f(x) en cualquier nodo del intervalo como f′′(xi) = N ∑ j=1 d2Lj(x) dx2x=xi fj= N ∑ j=1 L′′ j(xi)fji=1,...,N.(3.6) En efecto, se deduce que la derivada de orden s de f(x) en cualquier punto del intervalo, f(s)(xi), se puede calcular como f(s)(xi) = N ∑ j=1 L(s) j(xi)fji=1,...,N,(3.7) siendo L(s) j(xi) las derivadas de orden s de los interpolantes de Lagrange en un nodo del intervalo xi. La expresión anterior (3.7) puede escribirse en forma compacta en términos de vectores y matrices. Para ello, se definen 3.1 Método de colocación aplicado a problemas unidimensionales 19  f=     f1 f2 . . . fN      , f(s)=      f(s) 1 f(s) 2 . . . f(s) N       yDs(i,j) =       L(s) 1(x1)L(s) 2(x1)··· L(s) N(x1) L(s) 1(x2)L(s) 2(x2)··· L(s) N(x2) . . .. . ..... . . L(s) 1(xN)L(s) 2(xN)··· L(s) N(xN)       . Por tanto, las derivadas de cualquier orden s de f(x) en los puntos del intervalo [a,b] pueden obtenerse mediante  f(s)=Ds f.(3.8) 3.1.1 Nodos de Chebyshev En cuanto la selección de los N nodos x1,x2,...,xN pertenecientes al intervalo [a,b] , resulta muy importante elegirlos adecuadamente de manera que la función f(x) definida por la expresión (3.1) constituya una buena aproximación de la función incógnita del problema en todo el intervalo [a,b] . Así, el objetivo será tratar de conseguir que el error entre la función incógnita y su aproximación sea mínimo. Siendo F(x) la función incógnita y f(x) la aproximación de la misma, la teoría de la interpolación demuestra que existe un punto cdel intervalo [a,b]tal que para cualquier punto xde dicho intervalo se verifica la relación f(x)−F(x) = A(x) N!f(N)(c),(3.9) donde A(x) = (x−x1)(x−x2)...(x−xN) . De esta forma, para estimar el error cometido se necesitan cotas para la derivada f(N) y para el producto A(x) . Al ser este último un polinomio, su valor absoluto será máximo en algún punto del intervalo [a,b] y dicho valor máximo dependerá, entonces, de la elección de los nodos x1,x2,...,xN . Así, esta elección se hará de forma que el valor absoluto máximo sea lo menor posible. En efecto, definiendo la norma de A(x)como ∥A∥=m´ ax a≤x≤b|A(x)|,(3.10) el problema será encontrar la posición de los nodos de tal forma que la expresión anterior tenga valor mínimo. Chebyshev, quien consiguió resolver el problema, obtuvo que dicho polinomio de grado N debía ser el definido por la recurrencia TN(x) = 2xTN−1(x)−TN−2(x)∀N≥1,(3.11) siendo T1(x) = 1 y T2(x) = x . Así, los nodos x1,x2,...,xN deben coincidir con los ceros del polinomio TN(x), que vienen dados por: xi=a+b 2+a−b 2cosπ(i−1) N−1∀i=1,...,N.(3.12) De esta forma, los puntos x que verifiquen la expresión anterior (3.12) se denominarán nodos de Chebyshev. 20 Capítulo 3. Método numérico de resolución. El método de colocación Fenómeno de Runge De cara a analizar el efecto que tiene la posición de los nodos sobre la calidad de la interpolación de Lagrange de una función, se interpolará la función de Runge en el intervalo [−1,1] . La expresión de dicha función es f(x) = 1 1+25x2.(3.13) Para ello, se eligirán dos tipos de nodos de interpolación, llevándose a cabo dos interpolaciones diferentes. Por un lado, se calculará la interpolación empleando nodos equiespaciados. Por otro lado, se realizará una interpolación utilizando nodos de Chebyshev. Eligiendo una cantidad total de nodos ( N ) igual a 7 para ambas interpolaciones, se obtienen las siguientes gráficas: Figura 3.1 Interpolación con nodos equiespaciados y con nodos de Chebyshev para N=7. Se observa que el error cometido en el caso de la interpolación de Lagrange con nodos equiespaciados (línea verde) es considerablemente mayor que en el caso de interpolación con nodos de Chebyshev (línea roja). Esto se hace especialmente notable en las proximidades del extremos del intervalo, es decir, en torno a x=−1yx=1. Aumentando la cantidad N de nodos, se deduce que la interpolación con nodos de Chebyshev se convierte en una mejor aproximación, pues su gráfica coincide prácticamente en todos los puntos del intervalo con la de la función de Runge. La siguiente figura, donde se observa dicho suceso, muestra las gráficas de interpolación para el caso de 15 nodos de interpolación: 3.1 Método de colocación aplicado a problemas unidimensionales 21 Figura 3.2 Interpolación con nodos equiespaciados y con nodos de Chebyshev para N=15. En la figura (3.2) también destaca la enorme discrepancia entre la interpolación con nodos equiespaciados y la función de Runge, especialmente, en los nodos próximos a los extremos del intervalo, que es mucho mayor que en el caso de interpolación con N=7 . De esta forma, se deduce, tal y como descubrió Carl David Tolmé Runge, que cuanto mayor sea la cantidad de nodos equiespaciados utilizados para interpolar, lo que equivale a un polinomio de interpolación de mayor grado, mayor será el error cometido en la aproximación. Además, en la gráfica superior (3.2) también se puede observar que una posible solución a este fenómeno, denominado fenómeno de Runge, es la sustitución de los nodos equiespaciados por nodos de Chebyshev, consiguiendo que el error máximo disminuya al crecer el orden polinómico. Por último, es importante destacar que el uso de nodos de Chebyshev no será de utilidad para las interpolaciones de funciones en un intervalo con extremos en los que se den condiciones de periodicidad. Este será el caso de las interpolaciones en la coordenada θ de un sistema cilíndrico bidimensional, en las que bastará con utilizar nodos equiespaciados. Sin embargo, para las interpolaciones en la coordenada r de dicho sistema, sí será necesario emplear nodos de Chebyshev, pues dicha coordenada restá delimitada por dos extremos: RiyRe. 3.1.2 Cálculo de matrices de diferenciación empleando interpolación Lagrangiana En este apartado se presentará un procedimiento para el cálculo de matrices de diferenciación en mallas unidimensionales formadas por nodos arbitrarios utilizando la interpolación de Lagrange. Además, se incluirá su implementación en MATLAB. Para ello, se ha tomado como referencia el artículo A simple MATLAB program to compute differentiation matrices for arbitrary meshes via Lagrangian interpolation [ 8 ] escrito por el tutor del presente trabajo, Miguel Pérez-Saborid. En primer lugar, es importante definir el concepto de stencil. En matemáticas, concretamente en las áreas de análisis de las soluciones numéricas de las ecuaciones diferenciales parciales, un stencil consiste en una disposición geométrica de un conjunto de nodos que se relacionan con un determinado punto de interés mediante una técnica de aproximación numérica. Un 22 Capítulo 3. Método numérico de resolución. El método de colocación conocido ejemplo es el five-point stencil, que se detallará en el siguiente párrafo. Dada una cuadrícula, de una o dos dimensiones, de tal forma que tenga 4 cuadrados alrededor de cada vértice, un five-point stencil de un punto de dicha cuadrícula es un stencil compuesto por el propio punto de interés junto con sus cuatro puntos contiguos. Una de sus aplicaciones fundamentales es el cálculo numérico de derivadas en un punto de interés mediante diferencias finitas aplicadas en los puntos que forman parte del stencil asociado a dicho punto de interés. Figura 3.3 Representación de un five-point stencil en una y dos dimensiones, respectivamente. El procedimiento, que se detallará a continuación, se basa en el cálculo de matrices de diferenciación en una malla de N nodos arbitrariemente separados a partir de aquellas obtenidas mediante interpolación de Lagrange en diferentes stencils constituidos por un número arbitrario M≤Nde puntos contiguos del mallado. La elección de la interpolación de Lagrange como método de aproximación numérica, frente a, por ejemplo, las series de expansión de Taylor, radica en que el de Lagrange no da como resultado sistemas de ecuaciones de Vandermonde mal condicionados cuando el número de puntos del stencil es grande. Es decir, al utilizar la interpolación de Lagrange se asegura que un pequeño cambio en los coeficientes de las ecuaciones no da como resultado un cambio grande en la solución del problema. Matrices de diferenciación para un stencil Sea la interpolación de Lagrange de una función f(x) en un intervalo formado por un stencil de Mpuntos ¯x1,¯x2,...,¯xM: f(x) = M ∑ m=1 Lm(x)fm.(3.14) En la expresión anterior, ¯x1≤x≤¯xM , fm=f(¯xM) y Lm(x) es el interpolante de Lagrange asociado al punto del stencil ¯xM, tal que: Lm(x) = (x−¯x1)...(x−¯xm−1)(x−¯xm+1)...(x−¯xM) (¯xm−¯x1)...(¯xm−¯xm−1)(¯xm−¯xm+1)...(¯xm−¯xM).(3.15) A menudo, resulta muy conveniente expresar estos interpolantes de Lagrange (3.15) en la denominada forma baricéntrica: Lm(x) = a(x) (x−¯xm)am∀m=1,..,M.(3.16) 3.1 Método de colocación aplicado a problemas unidimensionales 23 En la forma baricéntrica de los interpolantes de Lagrange (3.16), el polinomio a(x) y la constante amse definen, respectivamente, de la siguiente forma: a(x) = M ∏ k=1 (x−¯xk),(3.17) am≡a′(¯xm) = ∏ k=m (¯xm−¯xk).(3.18) La forma baricéntrica es de gran utilidad cuando se quiere calcular las derivadas de Lm(x) en cualquier punto ¯xidel stencil. En efecto, reescribiendo la expresión (3.16) como a(x) = am(x−¯xm)Lm(x),(3.19) y tomando su derivada a′(x) = am(x−¯xm)L′ m(x)+ amLm(x),(3.20) se obtiene que para cualquier punto ¯xidistinto de ¯xmque L′ m(¯xi) = ai am(¯xi−¯xm)i=m,(3.21) puesto que, como se detalló previamente en este capítulo, Lm(¯xi) = 0 ya que m=i (según la relación entre el interpolante de Lagrange y la delta de Kronecker). Además, también se ha tenido en cuenta la definición (3.20) tal que a′(¯xi) = ai. Es importante destacar que la fórmula de la interpolación de Lagrange reproduce de manera exacta cualquier polinomio de grado r menor o igual que M−1 . En efecto, en estos casos, se verifican las siguientes relaciones, que serán muy útiles en lo que sigue a continuación: M ∑ m=1 Lm(x) = 1,(3.22) M ∑ m=1 L′ m(x) = 0.(3.23) Ambas relaciones se pueden demostrar fácilmente tomando, por ejemplo, la función f(x) = xr y su respectivo polinomio de Lagrange, teniéndose en cuenta que se debe de cumplir que 0≤r≤M−1. Para calcular L′ i(¯xi) , se empleará la expresión (3.23) particularizada en el punto ¯xi junto con (3.21), lo que da como resultado L′ i(¯xi) = − M ∑ m=1 m=i L′ m(¯xi) = − M ∑ m=1 m=i ai am(¯xi−¯xm).(3.24) Las derivadas de mayores órdenes de los interpolantes de Lagrange pueden ser calculadas emplenado la siguiente relación de recurrencia, que se deduce de tomar la derivada de orden s de la ecuación (3.19): a(s)(x) = am(x−¯xm)L(s) m(x)+ samL(s−1) m(x).(3.25) 24 Capítulo 3. Método numérico de resolución. El método de colocación Óbservese que cuando (3.25) se particulariza para x=¯xmda como resultado a(s)(¯xm) = samL(s−1) m(¯xm)∀m=1,...,M.(3.26) No obstante, para el caso de cualquier punto ¯xi=¯xm , la expresión (3.25) da como resultado: L(s) m(¯xi) = s (¯xi−¯xm)ai am L(s−1) i(¯xi)−L(s−1) m(¯xi)i=m.(3.27) De manera análoga a como se obtuvo el valor de L′ i(¯xi) , reflejado en la expresión (3.24), se obtiene que el valor de L(s) i(¯xi)es L(s) i(¯xi) = − M ∑ m=1 m=i L(s) m(¯xi).(3.28) Obsérvese que el caso s=1 , relativo a (3.21) y (3.24), está incluido en las expresiones anteriores (3.27)-(3.28) ya que L(0) m(¯xi)≡δim, donde δim es la delta de Kronecker. Estas cuatro últimas expresiones, mencionadas en el párrafo anterior, permiten el cálculo numérico de la derivada de orden sde cualquier interpolante de Lagrange en cualquier punto del stencil. En efecto, estas derivadas constituyen la matriz diferencial de orden s del stencil, denotada como ¯ Ds, que se define como ¯ Ds(i,m) = L(s) m(¯xi)∀i,m=1,...,M,(3.29) Expresando (3.27) y (3.28) en función de ¯ Ds, se obtienen, respectivamente: ¯ Ds(i,m) = s (¯xi−¯xm)ai am ¯ Ds−1(i,i)−¯ Ds−1(i,m)i=m,(3.30) ¯ Ds(i,i) = − M ∑ m=1 m=i ¯ Ds(i,m).(3.31) Finalmente, de acuerdo con (3.14) y (3.29), el valor de la derivada de orden s de f(x) en un punto ¯xidel stencil puede aproximarse como f(s)(¯xi) = M ∑ m=1 ¯ Ds(i,m)f(¯xm).(3.32) A continuación, se presenta una función de MATLAB denominada Dbs.m que implementa las ecuaciones (3.30)-(3.31) con el objetivo de obtener los elementos involucrados en las derivadas de f(x) en un punto ¯xi del stencil; es decir, los elementos de la fila i de las matrices ¯ Ds (siendo s=1,2,3 ). La salidas de la función serán tres vectores columnas con la siguiente forma: ¯ ds≡¯ Ds(i,1 : M)s=1,2,3.(3.33) Por último, es importante destacar que, debido a los errores de redondeo derivados del cálculo numérico, el cálculo de lo cocientes ai/am , que aparecen en la relación de recurrencia (3.30), debe realizarse cuidadosamente. Esto es especialmente relevante en el caso de stencils muy grandes y cuyos puntos están separados por distancias diferentes. 3.1 Método de colocación aplicado a problemas unidimensionales 25 Código 3.1 Cálculo de matrices de diferenciación para un stencil de M puntos en un punto ¯xi del mismo. 1function [db_1 db_2 db_3]=Dbs(i,xb,M) 2 3% Cálculo de los factores a_i/a_m y almacenamiento de los valores para cada m diferente a i 4 5aiam=zeros(1,M); 6F_i=xb(i)−xb; F_i(i)=1; n_i=length( find (F_i<0)); 7for m=1:M, 8if abs(m−i)>0, 9F_m=xb(m)−xb ; F_m(m)=1; n_m=length(find(F_m <0)); 10 aiam(m)=(−1)^(n_i−n_m)∗prod(sort(abs(F_i)) ./ sort (abs(F_m))); 11 end 12 end 13 14 % Matrices de diferenciacion : 15 16 Db_1=zeros(M,M); Db_2=zeros(M,M); Db_3=zeros(M,M); 17 18 % De primer orden, s =1: 19 20 for m=1:M 21 if abs(m−i)>0, 22 Db_1(i,m)=aiam(m)/(xb(i)−xb(m)); 23 end 24 end 25 Db_1(i,i)=−sum(Db_1(i,:)); 26 27 % De segundo orden, s=2 28 29 for m=1:M 30 if abs(m−i)>0, 31 Db_2(i,m) =2/(xb(i)−xb(m))∗(aiam(m)∗Db_1(i,i)−Db_1(i,m)); 32 end 33 end 34 Db_2(i,i)=−sum(Db_2(i,:)); 35 36 % De tercer orden, s=3: 37 38 for m=1:M 39 if abs(m−i)>0, 40 Db_3(i,m)=3/(xb(i)−xb(m))∗(aiam(m)∗Db_2(i,i)−Db_2(i,m)); 41 end 42 end 43 Db_3(i,i)=−sum(Db_3(i,:)); 44 45 % Salida 46 47 db_1=Db_1(i,:); db_2=Db_2(i,:); db_3=Db_3(i,:); 48 end 32 Capítulo 3. Método numérico de resolución. El método de colocación para problemas bidimensionales (en este caso, para x y z ), obteniéndose un procedimiento similar al caso unidimensional. 3.3 Ejemplos de aplicación en problemas bidimensionales 3.3.1 La ecuación de Poisson Sea la ecuación de Poisson en coordenadas cilíndricas bidimensionales ∇2φ=∂2φ ∂r2+1 r ∂φ ∂r+1 r2 ∂2φ ∂θ2=f(r,θ),(3.68) para la cavidad en forma de corona circular de radio interior Ri y radio exterior Re mostrada en la figura (2.1). El término fuente f(r,θ) es conocido y la ecuación debe resolverse sujeta a las siguientes cuatro condiciones de contorno: φ(r=Ri,θ) = Risin2θ,φ(r=Re,θ) = Resin2θ,(3.69) φ(r,θ=0) = φ(r,θ=2π),∂φ ∂θ (r,θ=0) = ∂φ ∂θ (r,θ=2π).(3.70) Las condiciones de contorno (3.69) son las relativas a la coordenada r , mientras que las (3.70) son las asociadas a θ . Es importante destacar que se impondrá que tanto la función φ como su derivada con respecto a θ tomen el mismo valor cuando se evalúen en θ=0 y θ=2π , debido a la periodicidad de la coordenada θ en dichos valores. Para ello, se aplican las condiciones de contorno (3.70), denominadas condiciones de periodicidad de la coordenada θ . En cuanto al mallado de la cavidad, se introducirá una distribución de Nr nodos de Chebyshev para la coordenada r (puesto que dicha coordenada está delimitada por dos extremos: Ri y Re ), mientras que para la coordenada θ se empleará una distribución de Nθ nodos equiespaciados (ya que en los valores extremos de dicha coordenada, θ=0 y θ=2π , se aplican condiciones de periodicidad). De esta forma, se obtiene un mallado con Nt=Nr×Nθ nodos (θi,rj) , similar a (3.4), que se muestra en la siguiente figura: Figura 3.6 Representación del mallado bidimensional en coordenadas θyr. 3.3 Ejemplos de aplicación en problemas bidimensionales 33 Haciendo uso de las expresiones (3.52) y (3.67), pero en este caso aplicándolas a las variables θyr, respectivamente, la ecuación de Poisson (3.68) puede discretizarse como D(2) r+1 rDr+1 r2D(2) θ φ= f,(3.71) donde los vectores  φ y  f tienen como componentes los valores de las funciones φ(r,θ) y f(r,θ)en los puntos I= (i−1)Nr +jdel mallado,  φx=       φ(1) . . . φ(I) . . . φ(Nt)        y f=       f(1) . . . f(I) . . . f(Nt)        .(3.72) Además, los términos 1 r y 1 r2 son matrices de dimensiones Ntx Nt , cuyas filas son todas iguales y están formadas por los valores de la coordenada r en cada nodo I del mallado. Así, definiendorcomo un vector fila de Ntelementos, r=r(1)··· r(I)··· r(Nt),(3.73) las expresiones de las matrices 1 ry1 r2son 1 r=       1\r(1)··· 1\r(I)··· 1\r(Nt) . . . 1\r(1)··· 1\r(I)··· 1\r(Nt) . . . 1\r(1)··· 1\r(I)··· 1\r(Nt)        ,(3.74) 1 r2=       1\r2(1)··· 1\r2(I)··· 1\r2(Nt) . . . 1\r2(1)··· 1\r2(I)··· 1\r2(Nt) . . . 1\r2(1)··· 1\r2(I)··· 1\r2(Nt)        .(3.75) De este modo, denominando como A al término que está multiplicando a φ en (3.71), que pasará a denominarse la matriz del sistema, la solución se obtendría resolviendo  φ=A\ f.(3.76) No obstante, antes de resolver el sistema algebraico (3.76) será necesario aplicarle las condiciones de contorno (3.69)-(3.70). Para aplicar las relativas a r , es decir, las condiciones ( 3.69 ), se modificarán la matriz del sistema A y el vector  f de forma que si el valor del índice I coincide con el de un nodo situado en un contorno horizontal, entonces todos los elementos de la fila I de la matriz A se anulan excepto el correspondiente a la diagonal, que tomará el valor unidad A(I,I) = 1 , y la componente f(I) de  f se impondrá igual a Risin2θ (si I está situado en el contorno inferior) o Resin2θ (si I está situado en el contorno superior). Por otro lado, para aplicar las condiciones de contorno 34 Capítulo 3. Método numérico de resolución. El método de colocación de periodicidad relativas a θ , se volverán a modificar la matriz del sistema A y el vector  f de manera que si el índice Icoincide con el de un nodo situado en: • El contorno vertical izquierdo, asociado a θ=0 , entonces todos los elementos de la fila I de la matriz A se anulan excepto el correspondiente a la diagonal, que tomará el valor unidad, A(I,I) = 1 , y el situado a Nt−Nr columnas a la derecha del elemento diagonal, que tomará el valor menos uno, A(I,I+Nt −Nr) = −1 . Por otro lado, la componente f(I)de  fse anulará. De esta forma, se consigue la igualdad de valores de la función φ cuando θes igual a 0y2π. • El contorno vertical derecho, asociado a θ=2π , entonces en la fila I de la matriz A se almacenará la diferencia entre los elementos correspondientes a la fila I y a la fila I−(Nθ−1)Nr de la matriz Dθ . Por otro lado, como en el caso anterior, la componente f(I) de  f se anulará. De esta forma, se consigue la igualdad de valores de la derivada parcial con respecto a θde la función φcuando θes igual a 0y2π. A continuación, se recopilan las expresiones de los índices de los nodos del mallado, I , donde se aplicarán las condiciones de contorno de este problema: Figura 3.7 Índices de los nodos donde se aplican las condiciones de contorno del problema. •Contorno vertical izquierdo →I=j,∀j=2,...,Nr−1. •Contorno vertical derecho →I= (Nθ−1)Nr+j,∀j=2,...,Nr−1. • Contorno horizontal superior (pared relativa a Re ) →I= (i−1)Nr+Nr , ∀i=1,...,Nθ . • Contorno horizontal inferior (pared relativa a Ri ) →I= (i−1)Nr+1 , ∀i=1,...,Nθ . La función fuente f(r,θ)elegida para la resolución del problema es f(r,θ) = 1 rsin2θ+2 rcos2θ,(3.77) que proporciona una solución analítica de la ecuación de Poisson (3.68), junto con las condiciones de contorno (3.69)-(3.70), igual a φ(r,θ) = rsin2θ.(3.78) Para la resolución numérica del problema se ha elaborado un programa de MATLAB, mostrado a continuación. Dicho programa implementa las funciones Dbs.m yDsmesh.m, mencionadas previamente, que permiten calcular las matrices de diferenciación necesarias en (3.71) y en la condición de contorno asociada a la derivada parcial con respecto a θde φ. 3.3 Ejemplos de aplicación en problemas bidimensionales 35 Código 3.3 Resolución numérica y analítica de la ecuación de Poisson en una corona circular. 1% Número de puntos del mallado 2Ntheta=30; 3Nr=30; 4Nt=Ntheta∗Nr; 5% Nodos equiespaciados para coordenada theta (debido a periodicidad en theta ) 6thetals =linspace (0,2∗pi ,Ntheta) ; 7% Nodos de Chebyshev para coordenada r (debido a las dos fronteras existentes ) 8Ri=2; Re=5; 9rch (1:Nr)=(Re+Ri)/2−(Re−Ri)/2∗cos(((1:Nr)−1)∗pi/(Nr−1)); 10 % Matrices de diferenciacion 11 M=9; % valor impar de M (y menor que Nr y Ntheta) 12 [ dtheta d2theta d3theta]=Dsmesh(Ntheta, thetals ,M); 13 [dr d2r d3r]=Dsmesh(Nr,rch,M); 14 Ltheta=eye(Ntheta); 15 Lr=eye(Nr); 16 for i=1:Ntheta 17 for j=1:Nr 18 I=(i−1)∗Nr+j; 19 r(I)=rch(j); 20 theta (I)= thetals (i); 21 for m=1:Ntheta 22 for n=1:Nr 23 K=(m−1)∗Nr+n; 24 % Matrices de diferenciacion de primer orden 25 Dtheta(I ,K)=dtheta(i ,m)∗Lr(j,n); 26 Dr(I,K)=Ltheta(i ,m)∗dr(j,n); 27 % Matrices de diferenciacion de segudno orden 28 Dtheta2(I ,K)=d2theta(i ,m)∗Lr(j,n); 29 Dr2(I,K)=Ltheta(i ,m)∗d2r(j,n); 30 31 end 32 end 33 end 34 end 35 rmat=zeros(Nt); 36 for i=1:Nt 37 rmat(i ,1: Nt)=r(1:Nt); 38 end 39 rmatelm1=1./rmat; % matriz 1/r(k) 40 rmatelm2=1./(rmat.^2) ; % matriz 1/( r^2)(k) 41 A=(rmatelm1.∗Dr)+Dr2+(rmatelm2.∗Dtheta2); 42 f =((1./ r (1:Nt).∗( sin( theta (1:Nt)) .^2) ) +(2./ r (1:Nt).∗cos(2.∗theta (1:Nt)))) '; 43 % CONDICIONES DE CONTORNO 44 % Condiciones de contorno relativas a r −contornos horizontales 45 for i=1:(Ntheta) 46 % Contorno inferior 47 I=(i−1)∗Nr+1; A(I ,:) =0; A(I,I)=1; f(I ,1)=Ri∗(sin( theta (I))) .^2; 48 % Contorno superior 49 I=(i−1)∗Nr+Nr; A(I,:)=0; A(I,I)=1; f(I ,1)=Re∗(sin( theta (I))) .^2; 50 end 51 % Condiciones de contorno de periodicidad relativas a theta −contornos 52 % verticales 36 Capítulo 3. Método numérico de resolución. El método de colocación 53 for j=2:(Nr−1) 54 % Contorno izquierdo 55 I=j; A(I ,:) =0; A(I,I)=1; A(I,I+Nt−Nr)=−1 ; f(I ,1)=0; 56 % Contorno derecho 57 I=(Ntheta−1)∗Nr+j; A(I ,:) =Dtheta(I ,:) −Dtheta(I−(Ntheta−1)∗Nr,:) ; f(I ,1)=0; 58 end 59 % Solución numérica de phi 60 phi=A\f; 61 for i=1:Ntheta 62 for j=1:Nr 63 thetamatrep (i , j)= thetals (i); 64 rmatrep(i , j)=rch(j); 65 phirep(i , j)=phi(( i−1)∗Nr+j,1); 66 end 67 end 68 % Representación de la soluci ón numérica en el mallado rectangular 69 figure (1) 70 mesh(thetamatrep,rmatrep, phirep) 71 zlim([0 5]) 72 colorbar 73 title ( 'Solución numérica en el mallado rectangular ') 74 % Representación de la soluci ón numérica en la corona circular 75 xp(1:Ntheta ,1: Nr)=rmatrep(1:Ntheta ,1: Nr).∗cos( thetamatrep (1:Ntheta ,1: Nr)); 76 yp(1:Ntheta ,1: Nr)=rmatrep(1:Ntheta ,1: Nr).∗sin ( thetamatrep (1: Ntheta ,1: Nr)); 77 figure (2) 78 mesh(xp,yp,phirep) 79 zlim([0 5]) 80 colorbar 81 title ( 'Solución numérica en la corona circular ') 82 % Solución analí tica 83 phianalitica =(r (1:Nt).∗( sin( theta (1:Nt)) .^2) ) '; 84 85 for i=1:Ntheta 86 for j=1:Nr 87 phianaliticarep (i , j)= phianalitica (( i−1)∗Nr+j,1); 88 end 89 end 90 % Representación de la soluci ón analí tica en el mallado rectangular 91 figure (3) 92 mesh(thetamatrep,rmatrep, phianaliticarep ) 93 colorbar 94 title ( 'Solución analí tica en el mallado rectangular ') 95 % Representación de la soluci ón analí tica en la corona circular 96 figure (4) 97 mesh(xp,yp, phianaliticarep ) 98 colorbar 99 title ( 'Solución analí tica en la corona circular ') 3.3 Ejemplos de aplicación en problemas bidimensionales 37 La geometría de la cavidad seleccionada para el programa de MATLAB está definida por un radio interior Ri=2 y un radio exterior Re=5 . En cuanto al mallado empleado, se ha realizado uno de 900 nodos con Nθ=30 y Nr=30 . Para la creación del mismo, se han tomado nodos equiespaciados en la variable θ y nodos de Chebyshev en la variable r ya que, como se ha explicado anteriormente, mientras que en la coordenada r existen dos fronteras que la delimitan (cuando r=Ri y r=Re ), en la coordenada θ no existe ninguna (cuando θ=0 y θ=2π se dan condiciones de periodicidad). Por otro lado, se ha elegido un valor de M impar e igual a 9 , con lo que M<Nθ,Nr. Tras ejecutar el programa se obtienen 4 gráficas diferentes. La primera y segunda gráfica representan la solución numérica del problema en el plano θ - r (mallado rectangular) y en la corona circular, respectivamente: ((a)) Solución numérica en el plano θr. ((b)) Solución numérica en la corona circular. Figura 3.8 Solución numérica de la ecuación de Poisson en una corona circular. Por otro lado, la tercera y cuarta gráfica representan la solución analítica, obtenida a partir de (3.78), en el plano θr y en la corona circular, respectivamente: ((a)) Solución analítica en el plano θr. ((b)) Solución analítica en la corona circular. Figura 3.9 Solución analítica de la ecuación de Poisson en una corona circular. 38 Capítulo 3. Método numérico de resolución. El método de colocación Comparando las figuras (3.8) y (3.9), se observa a simple vista que la solución numérica calculada mediante el método de colocación es prácticamente idéntica a la solución analítica del problema. Es importante destacar que cuanto mayor sea el número de puntos del mallado, mejor será la aproximación de la solución numérica a la analítica. 3.3.2 La ecuación de Laplace en régimen no estacionario Sea la ecuación de Laplace para una función φ(t,r,θ) en coordenadas cilíndricas bidimensionales y en régimen no estacionario, es decir, teniéndose en cuenta las variaciones temporales de la función φ mediante su derivada temporal. Considérese, además, la cavidad en forma de corona circular de radio interior Ri y radio exterior Re mostrada en la figura (2.1), que es la misma que se ha empleado para la resolución de la ecuación de Poisson (3.3.1). Así, el problema consistirá en resolver la ecuación ∇2φ=∂2φ ∂r2+1 r ∂φ ∂r+1 r2 ∂2φ ∂θ2=∂φ ∂t,(3.79) junto con las siguientes cuatro condiciones de contorno: φ(t,r=Ri,θ) = sin(θ),φ(t,r=Re,θ) = cos(θ),(3.80) φ(t,r,θ=0) = φ(t,r,θ=2π),∂φ ∂θ (t,r,θ=0) = ∂φ ∂θ (t,r,θ=2π).(3.81) Además, al aparecer en (3.79) la derivada parcial de φ con respecto a t , la resolución del problema requerirá una condición inicial: φ(t=0,r,θ) = 0.(3.82) Discretizando la ecuación (3.79) en el tiempo, es decir, aplicándola a un determinado intervalo de tiempo tn, se obtiene ∂2φn ∂r2+1 r ∂φn ∂r+1 r2 ∂2φn ∂θ2=φn−φn−1 ∆t,(3.83) donde ∆t es el tiempo transcurrido entre los instantes tn−1 y tn . En dicha discretización temporal se ha empleado la diferencia centrada, definida por ∂f ∂t=fn−fn−1 ∆tyf(t) = fn+fn−1 2,(3.84) donde f es una función que depende del tiempo, fn es el valor que toma dicha función en el instante tn, y fn−1es el valor que toma la función en el instante anterior tn−1. A continuación, discretizando la expresión anterior en los Nt puntos de colocación del mallado (3.6), de manera análoga a como se hecho en el caso de la ecuación de Poisson (3.3.1), se obtiene D(2) r+1 rDr+1  r2D(2) θ φn= φn− φn−1 ∆t.(3.85) La resolución del problema se va a llevar a cabo mediante una secuencia cíclica que permitirá iterar el problema en pequeños intervalos de tiempo. En primer lugar, se despeja la temperatura en un determinado instante n: 3.3 Ejemplos de aplicación en problemas bidimensionales 39 I−∆tD(2) r+1 rDr+1  r2D(2) θ φn= φn−1.(3.86) En segundo lugar, se aplicarán las condiciones de contorno (3.80)-(3.81). En este contexto, se introduce el vector −−→ FBC , denominado factor de condiciones de contorno, que permite imponer las condiciones de contorno en un vector que denotaremos como r sin que afecte a los valores de φdel instante anterior,  φn−1. Así, la ecuación (3.86) resultaría de la siguiente forma: I−∆tD(2) r+1 rDr+1  r2D(2) θ | {z } A  φn=−−→ FBC  φn−1+(1−−−→ FBC)r | {z }  b ,(3.87) donde −−→ FBC se define como un vector de Nt unos con elementos que toman el valor cero si sus índices coinciden con los de los nodos en los que se aplican las condiciones de contorno. De esta forma, denotando el término que está multiplicando a  φn en el miembro izquierdo de (3.87) como matriz del sistema A , la imposición de las condicones de contorno en este problema se llevará a cabo de manera similar al caso de la ecuación de Poisson (3.3.1). En efecto, se irá modificando la matriz Ay el vector r. Finalmente, se realiza un avance en el tiempo a la vez que se calcula la solución  φn para cada instante tn:  φn=A\ b.(3.88) Por ejemplo, la primera iteración consistirá en despejar el vector  φn correspondiente al instante posterior a t=0 utilizando el dato de φ proporcionado por la condición inicial (almacenado en el vector  φn−1 ). A continuación, para la segunda iteración, se almacenará el recién calculado  φn como  φn−1 y se repetirá el proceso. Así se procederá sucesivamente hasta obtener la evolución completa del problema en el intervalo de tiempo deseado. El siguiente programa de MATLAB implementa la resolución numérica de la ecuación de Laplace (3.79), junto con las condiciones de contorno (3.80)-(3.81). Dicho programa emplea las funciones Dbs.m yDsmesh.m para calcular las matrices de diferenciación involucradas en la resolución del problema. Código 3.4 Resolución numérica de la ecuación de Laplace en régimen no estacionario en una corona circular. 1% Número de puntos del mallado 2Ntheta=30; 3Nr=30; 4Nt=Ntheta∗Nr; 5% Nodos equiespaciados para coordenada theta (debido a periodicidad en theta ) 6thetals =linspace (0,2∗pi ,Ntheta) ; 7% Nodos de Chebyshev para r (debido a las dos fronteras existentes ) 8Ri=2; Re=5; 9rch (1:Nr)=(Re+Ri)/2−(Re−Ri)/2∗cos(((1:Nr)−1)∗pi/(Nr−1)); 10 % Matrices de diferenciacion 11 M=9; % valor impar de M (y menor que Nr y Ntheta) 12 [ dtheta d2theta d3theta]=Dsmesh(Ntheta, thetals ,M); 13 [dr d2r d3r]=Dsmesh(Nr,rch,M); 40 Capítulo 3. Método numérico de resolución. El método de colocación 14 Ltheta=eye(Ntheta); 15 Lr=eye(Nr); 16 for i=1:Ntheta 17 for j=1:Nr 18 I=(i−1)∗Nr+j; 19 r(I)=rch(j); 20 theta (I)= thetals (i); 21 for m=1:Ntheta 22 for n=1:Nr 23 K=(m−1)∗Nr+n; 24 % Matrices de diferenciacion de primer orden 25 Dtheta(I ,K)=dtheta(i ,m)∗Lr(j,n); 26 Dr(I,K)=Ltheta(i ,m)∗dr(j,n); 27 % Matrices de diferenciacion de segundo orden 28 Dtheta2(I ,K)=d2theta(i ,m)∗Lr(j,n); 29 Dr2(I,K)=Ltheta(i ,m)∗d2r(j,n); 30 end 31 end 32 end 33 end 34 rmat=zeros(Nt); 35 for i=1:Nt 36 rmat(i ,1: Nt)=r(1:Nt); 37 end 38 rmatelm1=1./rmat; % matriz 1/r(k) 39 rmatelm2=1./(rmat.^2) ; % matriz 1/( r^2)(k) 40 dt=0.01; % incremento de tiempo en cada iteraci ón (caso no estacionario ) 41 A=eye(Nt)−dt∗(rmatelm1.∗Dr+Dr2+rmatelm2.∗Dtheta2); 42 rizq=zeros(Nt,1) ; % vector para incorporar las condiciones de contorno 43 FBC(1:Nt,1)=1; % vector de factor de condiciones de contorno 44 % CONDICIONES DE CONTORNO 45 % Condiciones de contorno relativas a r −contornos horizontales 46 for i=1:(Ntheta) 47 % Contorno inferior 48 I=(i−1)∗Nr+1; A(I ,:) =0; A(I,I)=1; rizq (I ,1)=sin( theta (I)) ; FBC(I,1)=0; 49 % Contorno superior 50 I=(i−1)∗Nr+Nr; A(I,:) =0; A(I,I)=1; rizq (I ,1)=cos( theta (I)); FBC(I,1)=0; 51 end 52 % Condiciones de contorno de periodicidad relativas a theta −contornos 53 % verticales 54 for j=2:(Nr−1) 55 % Contorno izquierdo 56 I=j; A(I ,:) =0; A(I,I)=1; A(I,I+Nt−Nr)=−1 ; rizq(I ,1)=0; FBC(I,1)=0; 57 % Contorno derecho 58 I=(Ntheta−1)∗Nr+j; A(I ,:) =Dtheta(I ,:) −Dtheta(I−(Ntheta−1)∗Nr,:) ; FBC(I,1)=0; 59 end 60 % Solución numérica del vector phi 61 % Caso no estacionario : se tienen en cuenta las variaciones temporales de T 62 phinm1(1:Nt,1)=0; % condición inicial 63 Ntime=2000; 64 for n=1:Ntime 65 t=dt∗n 66 b=FBC.∗phinm1+(1−FBC).∗rizq; 67 phi=A\b; 3.3 Ejemplos de aplicación en problemas bidimensionales 41 68 for i=1:Ntheta 69 for j=1:Nr 70 thetamatrep (i , j)= thetals (i); 71 rmatrep(i , j)=rch(j); 72 phirep(i , j)=phi(( i−1)∗Nr+j,1); 73 end 74 end 75 % Representación en el mallado rectangular 76 figure (1) 77 mesh(thetamatrep,rmatrep, phirep) 78 colorbar 79 title ( 'Representación en el mallado rectangular ') 80 % Representación en la corona circular 81 xp(1:Ntheta ,1: Nr)=rmatrep(1:Ntheta ,1: Nr).∗cos( thetamatrep (1:Ntheta ,1: Nr)); 82 yp(1:Ntheta ,1: Nr)=rmatrep(1:Ntheta ,1: Nr).∗sin ( thetamatrep (1: Ntheta ,1: Nr)); 83 figure (2) 84 mesh(xp,yp,phirep) 85 colorbar 86 title ( 'Representación en la corona circular ') 87 pause 88 hold off 89 phinm1=phi; % actualizaci ón del vector phinm1 90 end 48 Capítulo 3. Método numérico de resolución. El método de colocación De manera análoga, para la variable ϕse tienen las expresiones ∂ϕ∗ ∂r∗I = Nt ∑ K=1 Dr∗(I,K)ϕ∗(K,t∗)≡(Dr ∗phi)I,(3.101) ∂2ϕ∗ ∂r∗2I = Nt ∑ K=1 D(2) r∗(I,K)ϕ∗(K,t∗)≡(Dr2∗phi)I,(3.102) ∂ϕ∗ ∂θ∗I = Nt ∑ K=1 Dθ∗(I,K)ϕ∗(K,t∗)≡(Dtheta ∗phi)I,(3.103) ∂2ϕ∗ ∂θ∗2I = Nt ∑ K=1 D(2) θ∗(I,K)ϕ∗(K,t∗)≡(Dtheta2∗phi)I.(3.104) A continuación, se definen las matrices 1 r∗ y 1 r∗2 , de manera similar a como se hizo en (3.74) y (3.75): 1 r=       1\r∗(1)··· 1\r∗(I)··· 1\r∗(Nt) . . .. . .. . .. . .. . . 1\r∗(1)··· 1\r∗(I)··· 1\r∗(Nt) . . .. . .. . .. . .. . . 1\r∗(1)··· 1\r∗(I)··· 1\r∗(Nt)       ≡rmatelm1, 1 r2=         1\r∗2(1)··· 1\r∗2(I)··· 1\r∗2(Nt) . . .. . .. . .. . .. . . 1\r∗2(1)··· 1\r∗2(I)··· 1\r∗2(Nt) . . .. . .. . .. . .. . . 1\r∗2(1)··· 1\r∗2(I)··· 1\r∗2(Nt)         ≡rmatelm2.(3.105) De esta forma, las expresiones discretizadas en el espacio de los operadores laplaciano, ∇2 , y biarmónico, ∇4, son, respectivamente, DL =Dr2+rmatelm1.∗Dr +rmatelm2.∗Dtheta2,(3.106) DL2=DL ∗DL.(3.107) Por tanto, se obtienen las expresiones ∇2Ψ∗I = Nt ∑ K=1D(2) r∗(I,K)+ 1 rDr∗(I,K)+ 1 r2D(2) θ∗(I,K)Ψ∗(K,t∗)≡(DL ∗psi)I,(3.108) ∇4Ψ∗I = Nt ∑ K=1D(2) r∗(I,K)+ 1 rDr∗(I,K)+ 1 r2D(2) θ∗(I,K)2 Ψ∗(K,t∗)≡(DL2∗psi)I,(3.109) 3.4 Implementación numérica de las ecuaciones de Saltzman adimensionales 49 ∇2ϕ∗I = Nt ∑ K=1D(2) r∗(I,K)+ 1 rDr∗(I,K)+ 1 r2D(2) θ∗(I,K)ϕ∗(K,t∗)≡(DL ∗phi)I,(3.110) Finalmente, se procede a discretizar espacialmente las ecuaciones (3.93)-(3.94), que ya habían sido discretizadas temporalmente, empleando las expresiones que se acaban de obtener. Así, se obtiene el siguiente sistema de ecuaciones en forma compacta, es decir, generalizado a todos los puntos del mallado: (DL ∗psin)−∆t∗∗rPr Ra ∗(DL2∗psin)+∆t∗∗rmatelm1∗(Dtheta ∗psin) = DL ∗psin−1+∆t∗∗NLpsin−1, phin+∆t∗∗1 ln(λ)∗rmatelm2∗(Dtheta ∗psin)−∆t∗∗1 √PrRa ∗(DL ∗phin) = phin−1+∆t∗∗,NLphin−1, donde NLpsin−1 y NLphin−1 son los términos no lineales de las ecuaciones (3.93) y (3.94), respectivamente, cuyas expresiones son NLpsin−1=rmatelm1∗[(Dr ∗DL ∗psin−1).∗(Dtheta ∗psin−1)− (Dtheta ∗DL ∗phin−1).∗(Dr ∗psin−1)], NLphin−1=rmatelm1∗[(Dtheta ∗psin−1).∗(Dr ∗phin−1)− (Dr ∗psin−1).∗(Dtheta ∗phin−1)]. Una vez obtenido el sistema de ecuaciones de Saltzman adimensional discretizado tanto espacial, en el mallado de Nt nodos (θi,rj) , como temporalmente, en un intervalo de tiempo definido por t∗ n , se plantea un sistema de matrices agrupando los coeficientes de los términos lineales, asociados al instante t∗ n , en una matriz que irá multiplicando a un vector columna con las variables a resolver ( Ψ∗ n , ϕ∗ n ), mientras que la matriz correspondiente al término independiente incorporará, además de los componentes lineales asociados al instante t∗ n−1 , los componentes no lineales NLpsin−1yNLphin−1. Así, el sistema matricial quedaría de la siguiente forma: Apsi Aphi Bpsi Bphi psin phin=DL ∗psin−1+∆t∗∗NLpsin−1 phin−1+∆t∗∗NLphin−1, donde Apsi =DL −∆t∗∗rPr Ra ∗DL2yAphi =DL∆t∗∗rmatelm1∗Dtheta, Bpsi =∆t∗∗1 ln(λ)∗rmatelm2∗Dtheta yBphi =I−∆t∗∗1 √PrRa ∗DL, 50 Capítulo 3. Método numérico de resolución. El método de colocación siendo I una matriz identidad de orden Nt. Lo siguiente será definir las condiciones de contorno del problema. Para ello, es necesario introducir el factor de condiciones de contorno, que ya se empleó para la resolución de la ecuación de Laplace no estacionaria (3.3.2). Sin embargo, en este caso, al tratarse de un sistema de dos ecuaciones con dos incógnitas, se introducirá un FBC para Ψ , −−−→ FBCΨ , y otro para ϕ , −−−→ FBCϕ. Apsi Aphi Bpsi Bphi psin phin="−−−→ FBCΨ −−−→ FBCϕ#.∗DL ∗psin−1+∆t∗∗NLpsin−1 phin−1+∆t∗∗NLphin−1+"(1−−−−→ FBCΨ) (1−−−−→ FBCϕ)#.∗rΨ rϕ, donde rΨ y rϕ son vectores columnas que permitirán almacenar las condiciones de contorno. De este modo, al tratarse de un problema bidimensional, será necesario definir condiciones de contorno tanto para Ψ∗como para ϕ∗. En primer lugar, se detallarán las de Ψ∗ . En el sistema de ecuaciones del problema a resolver aparece el biarmónico de Ψ∗ , ∇4Ψ∗ , que involucra la derivada cuarta de Ψ∗ con respecto a r∗ y la derivada cuarta de dicha variable con respecto a θ∗ , como se puede ver en la definición del operador biarmónico (3.96). Por tanto, para definir el problema son necesarias cuatro condiciones de contorno de Ψ∗ en r∗ y otras cuatro de dicha variable en θ∗ . Así, se impondrá que: Ψ∗(t∗,r∗=1,θ∗) = 0,Ψ∗t∗,r∗=Re Ri =λ,θ∗=0,(3.111) ∂Ψ∗ ∂r∗(t∗,r∗=1,θ∗) = 0,∂Ψ∗ ∂r∗t∗,r∗=Re Ri =λ,θ∗=0,(3.112) Ψ∗(t∗,r∗,θ∗=0) =Ψ∗(t∗,r∗,θ∗=2π),(3.113) ∂Ψ∗ ∂θ∗(t∗,r∗,θ∗=0) = ∂Ψ∗ ∂θ∗(t∗,r∗,θ∗=2π),(3.114) ∂2Ψ∗ ∂θ∗2(t∗,r∗,θ∗=0) = ∂2Ψ∗ ∂θ∗2(t∗,r∗,θ∗=2π),(3.115) ∂3Ψ∗ ∂θ∗3(t∗,r∗,θ∗=0) = ∂3Ψ∗ ∂θ∗3(t∗,r∗,θ∗=2π).(3.116) Las condiciones de contorno (3.111)-(3.112) son las asociadas a r∗ , mientras que (3.113)- (3.116) se corresponden con la coordenada θ∗ . La expresión (3.111) hace referencia a que la función de corriente adimensionalizada se impone nula en la pared interior y exterior de la corona circular. Además, en dichas paredes también se impone nula la componente en la dirección θ∗ de la velocidad , vθ∗ , mediante la expresión (3.112) (hay que recordar que la relación entre la función de corriente Ψ∗ y la componente de la velocidad según la dirección θ∗es vθ∗=∂Ψ∗/r∗). 3.4 Implementación numérica de las ecuaciones de Saltzman adimensionales 51 Por otro lado, las condiciones de contorno en θ∗ son las denominadas condiciones de periocidad de Ψ∗ , análogas a las aplicadas en los problemas de la ecuación de Poisson y Laplace no estacionaria, pero, en este caso, incluyendo, además, las igualdades de las derivadas segundas (3.115) y terceras (3.116), debido a la presencia del biarmónico de la función de corriente adimensionalizada en las ecuaciones del problema. En cuanto a las condiciones de contorno de ϕ∗ , se impondrán dos condiciones de contorno de dicha variable en r∗ y otras dos, de la misma variable, en θ∗ , para poder resolver el problema. Esto se debe a que, en el sistema de ecuaciones del problema a resolver, aparece el laplaciano de ϕ∗ , ∇2ϕ∗ , que involucra la derivada segunda de ϕ∗ con respecto a r∗ y la derivada segunda de dicha variable con respecto a θ∗ , como se puede ver en la definición del operador laplaciano (3.95). Así, se impondrá que: ϕ∗(t∗,r∗=1,θ∗) = 0,ϕ∗t∗,r∗=Re Ri =λ,θ∗=0,(3.117) ϕ∗(t∗,r∗,θ∗=0) = ϕ∗(t∗,r∗,θ∗=2π),(3.118) ∂ϕ∗ ∂θ∗(t∗,r∗,θ∗=0) = ∂ϕ∗ ∂θ∗(t∗,r∗,θ∗=2π),(3.119) Las condiciones de contorno (3.117) son las asociadas a r∗ , mientras que (3.118)-(3.119) se corresponden con la coordenda θ∗ . La expresión (3.117) hace referencia a que la componente de la temperatura adimensionalizada asociada al movimiento se impone nula tanto en la pared interior como en la exterior de la corona circular. Por otro lado, las condiciones de contorno en θ∗ son las denominadas condiciones de periodicidad de ϕ∗ , que, en este caso, no involucran las derivadas segundas y terceras con respecto a θ∗ ya que en las ecuaciones del problema aparece el laplaciano de ϕ∗y no el biarmónico de dicha función. Lo siguiente será detallar la implementación de estas condiciones de contorno para la resolución numérica del problema en MATLAB. Para ello, se irán modificando las submatrices Apsi , Aphi , Bpsi y Bphi , y los vectores rΨ y rϕ . De este modo, para aplicar las condiciones de contorno asociadas a r∗de: •Ψ∗=0 y ϕ∗=0 , (3.69) y (3.117), respectivamente, se emplearán los contornos horizontales superior e inferior del mallado. Así, para el caso de Ψ∗=0 , si el valor del índice I coincide con el de un nodo situado en un contorno horizontal, entonces todos los elementos de la fila I de la submatriz Apsi se anulan excepto el correspondiente a la diagonal, que tomará el valor unidad, Apsi(I,I) = 1 . Además, también se anularán todos los elementos de la fila I de la submatriz Aphi y la componente rΨ(I) de rΨ . Para el caso de ϕ∗=0 , si el valor del índice I coincide con el de un nodo situado en un contorno horizontal, entonces todos los elementos de la fila I de la submatriz Bphi se anulan excepto el correspondiente a la diagonal, que tomará el valor unidad, Bphi(I,I) = 1 . Además, también se anularán todos los elementos de la fila I de la submatriz Bpsi y la componente rϕ(I)de rϕ. •∂Ψ∗ ∂r∗=0 , (3.112), se emplearán los subcontornos horizontales superior e inferior del mallado. Así, si el valor del índice I coincide con el de un nodo situado en un subcontorno horizontal, entonces en la fila I de la submatriz Apsi se almacenarán los elementos correspondientes a la fila I−1 (si se trata del subcontorno inferior) o I+1 (si se trata del 52 Capítulo 3. Método numérico de resolución. El método de colocación subcontorno superior) de la matriz Dr . Además, también se anularán todos los elementos de la fila Ide la submatriz Aphi y la componente rΨ(I)de rΨ. Por otro lado, para aplicar las condiciones de contorno asociadas a θ∗ de periodicidad, en los valores extremos 0y2π, para: •Ψ∗ , (3.113), y ϕ∗ , (3.118), se empleará el contorno vertical izquierdo del mallado. Así, para la variable Ψ∗ , si el valor del índice I coincide con el de un nodo situado en dicho contorno, entonces todos los elementos de la fila I de la submatriz Apsi se anulan excepto el correspondiente a la diagonal, que tomará el valor unidad, Apsi(I,I) = 1 , y el situado a Nt−Nr columnas a la derecha del elemento diagonal, que tomará el valor menos uno, Apsi(I,I+Nt−Nr) = −1 . Además, también se anularán todos los elementos de la fila I de la submatriz Aphi y la componente rΨ(I) de rΨ . Para el caso de la variable ϕ∗ , si el valor del índice I coincide con el un nodo situado en el contorno vertical izquierdo, entonces todos los elementos de la fila I de la submatriz Bphi se anulan excepto el correspondiente a la diagonal, que tomará el valor unidad, Bphi(I,I) = 1 , y el situado a Nt−Nr columnas a la derecha del elemento diagonal, que tomará el valor menos uno, Bphi(I,I+Nt−Nr) = −1 . Además, también se anularán todos los elementos de la fila I de la submatriz Bpsi y la componente rϕ(I)de rϕ. •∂Ψ∗ ∂θ∗ , (3.114), y ∂ϕ∗ ∂θ∗ , (3.119), se empleará el contorno vertical derecho. Así, para el caso de la variable Ψ∗ , si el valor del índice I coincide con el de un nodo situado en dicho contorno, entonces en la fila I de la submatriz Apsi se almacenará la diferencia entre los elementos correspondientes a las filas IyI−(Nθ−1)Nrde la matriz Dtheta. Además, se anularán los elementos de la fila I de la submatriz Aphi y la componente rΨ(I) de rΨ . Para el caso de la variable ϕ∗ , si el valor del índice I coincide con el de un nodo situado en el contorno vertical derecho, entonces en la fila I de la submatriz Bphi se almacenará la diferencia entre los elementos correspondientes a las filas I y I−(Nθ−1)Nr de la matriz Dtheta . Además, se anularán los elementos de la fila I de la submatriz Bpsi y la componente rϕ(I)de rϕ. •∂2Ψ∗ ∂θ∗2 , (3.115), y ∂3Ψ∗ ∂θ∗3 , (3.116), se emplearán el subcontorno vertical izquierdo y derecho, respectivamente. Así, para el caso de la derivada de segundo orden, si el valor del índice I coincide con el de un nodo situado en el subcontorno vertical izquierdo, entonces en la fila I de la submatriz Apsi se almacenará la diferencia entre los elementos correspondientes a las filas I y I+(Nθ−2)Nr de la matriz Dtheta2 . Además, se anularán los elementos de la fila I de la submatriz Aphi y la componente rΨ(I) de rΨ . Para el caso de la derivada de tercer orden, si el valor del índice I coincide con el de un nodo situado en el subcontorno vertical derecho, entonces en la fila I de la submatriz Apsi se almacenará la diferencia entre los elementos correspondientes a las filas I y I−(Nθ−2)Nr de la matriz Dtheta3 . Además, se anularán los elementos de la fila I de la submatriz Aphi y la componente rΨ(I)de rΨ. A continuación, se recopilan las expresiones de los índices de los nodos del mallado, I , donde se aplicarán las condiciones de contorno del problema: 3.4 Implementación numérica de las ecuaciones de Saltzman adimensionales 53 Figura 3.18 Índices de los nodos donde se aplican las condiciones de contorno. •Contorno vertical izquierdo →I=j,∀j=2,...,Nr−1. •Subcontorno vertical izquierdo →I=j+Nr,∀j=3,...,Nr−2. •Contorno vertical derecho →I= (Nθ−1)Nr+j,∀j=2,...,Nr−1. •Subcontorno vertical derecho →I= (Nθ−2)Nr+j,∀j=3,...,Nr−2. • Contorno horizontal superior (pared relativa a r∗ e ) →I= (i−1)Nr+Nr , ∀i=1,...,Nθ . •Subcontorno horizontal superior →I= (i−1)Nr+(Nr−1),∀i=2,...,Nθ−1. • Contorno horizontal inferior (pared relativa a r∗ i ) →I= (i−1)Nr+1 , ∀i=1,...,Nθ . •Subcontorno horizontal inferior →I= (i−1)Nr+2,∀i=2,...,Nθ−1. Finalmente, al ser un problema no estacionario bidimensional, es necesario aplicar una condición inicial para cada variable involucrada en el problema (Ψ∗yϕ∗): Ψ∗(t∗=0,r∗,θ∗) = 0,05 sin2(θ∗)(r∗−r∗ i)2(r∗−r∗ e)2,(3.120) ϕ∗(t∗=0,r∗,θ∗) = 0,05 sin(θ∗)(r∗−r∗ i)(r∗−r∗ e).(3.121) Se puede comprobar fácilmente que ambas condiciones iniciales cumplen las condiciones de contorno del problema, pues han sido seleccionadas de tal forma que produzcan una pequeña perturbación con respecto a la situación de equilibrio hidrostático del problema, permitiendo estudiar la evolución de un estado inicialmente en equilibrio. A continuación, se presenta un programa de MATLAB que resuelve numéricamente las ecuaciones de Saltzman adimensionales, implementadas según lo desarrollado en este apartado: Código 3.5 Resolución numérica de las ecuaciones de Saltzman adimensionales en una corona circular. 1% Propiedades del fluido 2Pr=0.733; % Número de Prandtl 3Ra=1000; % Número de Rayleigh 4% Relación de aspecto −lambda 5lambda=2; 54 Capítulo 3. Método numérico de resolución. El método de colocación 6% Número de puntos del mallado 7Ntheta=40; 8Nr=40; 9Nt=Ntheta∗Nr; 10 % Nodos equiespaciados para theta (debido a periodicidad en theta ) 11 thetals =linspace (0,2∗pi ,Ntheta) ; 12 % Nodos de Chebyshev para r (debido a las dos fronteras existentes ) 13 re=lambda; 14 ri=1; 15 rch (1:Nr)=(re+ri )/2−(re−ri)/2∗cos (((1: Nr)−1)∗pi/(Nr−1)); 16 % Matrices de diferenciacion 17 M=9; % valor impar de M (y menor que Nr y Ntheta) 18 [ dtheta d2theta d3theta]=Dsmesh(Ntheta, thetals ,M); 19 [dr d2r d3r]=Dsmesh(Nr,rch,M); 20 Ltheta=eye(Ntheta); 21 Lr=eye(Nr); 22 Dtheta=sparse(Nt,Nt); Dtheta2=sparse(Nt,Nt);Dtheta3=sparse(Nt,Nt); 23 Dr=sparse(Nt,Nt); Dr2=sparse(Nt,Nt); 24 DL=sparse(Nt,Nt); DL2=sparse(Nt,Nt); 25 Apsi=sparse(Nt,Nt); Aphi=sparse(Nt,Nt); Bpsi=sparse(Nt,Nt); Bphi=sparse(Nt,Nt); 26 Asistema=sparse(2∗Nt,2∗Nt); 27 Dtheta=kron(dtheta ,speye(Nr)); Dr=kron(speye(Ntheta) ,dr); 28 Dtheta2=kron(d2theta ,speye(Nr)); Dr2=kron(speye(Ntheta),d2r) ; 29 Dtheta3=kron(d3theta ,speye(Nr)); 30 for i=1:Ntheta % recorro eje horizontal 31 for j=1:Nr % recorro eje vertical 32 % Calculamos I (punto de colocación) −existen Nt puntos de 33 % colocación 34 I=(i−1)∗Nr+j; 35 r(I)=rch(j); % vector de Nt filas 36 theta (I)= thetals (i); % vector de Nt filas 37 end 38 end 39 rmat=zeros(Nt); 40 for i=1:Nt 41 rmat(i ,1: Nt)=r(1:Nt); 42 end 43 relm1=(1./ r) ';% vector 1/r(k) de dimensiones (Nt,1) 44 rmatelm1=1./rmat; % matriz 1/r(k) 45 rmatelm2=1./(rmat.^2) ; % matriz 1/( r^2)(k) 46 dt=0.01; % incremento de tiempo en cada iteracci ón (caso no estacionario ) 47 DL=rmatelm1.∗Dr+Dr2+rmatelm2.∗Dtheta2; % operador laplaciano (en coordenadas cilí ndricas 2D) 48 DL2=DL∗DL; % operador biarmónico (en coordenadas cilíndricas 2D) 49 % Matrices Apsi, Aphi, Bpsi, y Bphi 50 Apsi=DL−(dt∗sqrt(Pr/Ra)∗DL2); 51 Aphi=dt∗rmatelm1.∗Dtheta; 52 Bpsi=dt/log(lambda)∗rmatelm2.∗Dtheta; 53 Bphi=eye(Nt)−(dt/ sqrt (Pr∗Ra)∗DL); 54 FBCpsi(1:Nt,1)=1; % vector de factor de condiciones de contorno para variable psi 55 FBCphi(1:Nt,1)=1; % vector de factor de condiciones de contorno para variable phi 56 % CONDICIONES DE CONTORNO 57 % Condiciones de contorno de psi y phi relativas a r en contornos horizontales 58 for i=1:(Ntheta) 3.4 Implementación numérica de las ecuaciones de Saltzman adimensionales 55 59 % Contorno inferior ( j=1,r=1) 60 I=(i−1)∗Nr+1; Apsi(I ,:) =0; Apsi(I, I)=1; Aphi(I ,:) =0; FBCpsi(I,1)=0; 61 Bpsi(I ,:) =0; Bphi(I ,:) =0; Bphi(I, I)=1; FBCphi(I,1)=0; 62 % Contorno superior ( j=Nr,r=Re/Ri=lambda) 63 I=(i−1)∗Nr+Nr; Apsi(I ,:) =0; Apsi(I ,I)=1; Aphi(I ,:) =0; FBCpsi(I,1)=0; 64 Bpsi(I ,:) =0; Bphi(I ,:) =0; Bphi(I, I)=1; FBCphi(I,1)=0; 65 end 66 % Condiciones de contorno de derivada psi con respecto a r en subcontornos horizontales 67 for i=2:(Ntheta−1) 68 % Subcontorno inferior 69 K=(i−1)∗Nr+2; Apsi(K,:)=Dr(K−1,:); Aphi(K,:)=0; FBCpsi(K,1)=0; 70 % Subcontorno superior 71 K=(i−1)∗Nr+(Nr−1); Apsi(K,:)=Dr(K+1,:); Aphi(K,:)=0; FBCpsi(K,1)=0; 72 end 73 % Condiciones de contorno de periodicidad de psi y phi relativas a theta −contornos 74 % verticales y subcontornos verticales 75 for j=2:(Nr−1) 76 % Contorno vertical izquierdo 77 I=j; Apsi(I ,:) =0; Apsi(I , I)=1; Apsi(I , I+Nt−Nr)=−1; Aphi(I,:)=0; 78 FBCpsi(I,1)=0; Bpsi(I ,:) =0; Bphi(I ,:) =0; Bphi(I, I)=1; Bphi(I, I+Nt−Nr)=−1; 79 FBCphi(I,1)=0; 80 % Contorno vertical derecho 81 I=(Ntheta−1)∗Nr+j; Apsi(I ,:) =Dtheta(I ,:) −Dtheta(I−(Ntheta−1)∗Nr,:) ; Aphi(I ,:) =0; FBCpsi(I,1)=0; Bpsi(I ,:) =0; Bphi(I ,:) =Dtheta(I ,:) −Dtheta(I−(Ntheta−1)∗Nr,:) ; FBCphi(I,1)=0; 82 end 83 for j=3:(Nr−2) 84 % Subcontorno vertical izquierdo 85 I=j+Nr; Apsi(I ,:) =Dtheta2(I−Nr,:)−Dtheta2(j+(Ntheta−1)∗Nr,:) ; Aphi(I ,:) =0 ; 86 FBCpsi(I,1)=0; 87 % Subcontorno vertical derecho 88 I=(Ntheta−2)∗Nr+j; Apsi(I ,:) =Dtheta3(I+Nr,:)−Dtheta3(j ,:) ; Aphi(I ,:) =0 ; 89 FBCpsi(I,1)=0; 90 end 91 % Matriz del sistema 92 Asistema=[Apsi, Aphi; Bpsi, Bphi]; 93 [L,U,P]=lu(Asistema); 94 % Solución numérica de psi y phi 95 % Caso no estacionario : se tienen en cuenta las variaciones temporales 96 for i=1:Ntheta 97 for j =1:Nr 98 I=(i−1)∗Nr+j; 99 psinm1(I,1)=0.05∗sin( thetals (i))^2∗(rch(j)−ri)^2∗(rch(j)−re)^2; % condición inicial para psi (cumple con las condiciones de contorno) 100 phinm1(I,1)=0.05∗sin( thetals (i))∗(rch(j)−ri)∗(rch(j)−re); % condición inicial para phi (cumple con las condiciones de contorno) 101 end 102 end 103 Ntime=100000; 104 A1=rmatelm1.∗Dtheta; 105 B1=Dr∗DL; 106 A2=A1∗DL; 107 for n=1:Ntime 108 t=dt∗n 56 Capítulo 3. Método numérico de resolución. El método de colocación 109 bpsi=dt∗((B1∗psinm1).∗(A1∗psinm1))−dt∗((A2∗psinm1).∗(Dr∗psinm1))+(DL∗psinm1); 110 bphi=dt∗((A1∗psinm1).∗(Dr∗phinm1))−dt∗((Dr∗psinm1).∗(A1∗phinm1))+phinm1; 111 b=[bpsi;bphi ].∗[ FBCpsi;FBCphi]; 112 b1=L\(P∗b); 113 sol=U\b1; 114 psi (1:Nt,1)=sol (1:Nt,1) ; 115 phi(1:Nt,1)=sol(Nt+1:2∗Nt,1); 116 % Organizar soluci ón para pintar 117 if round(n/10000)∗10000 == n 118 velocidad_r=−rmatelm1.∗Dtheta∗psi; 119 for i=1:Ntheta 120 for j=1:Nr 121 thetamatrep (i , j)= thetals (i); 122 rmatrep(i , j)=rch(j); 123 psimatrep(i , j)=psi (( i−1)∗Nr+j,1); % valores de psi en el instante t 124 phimatrep(i , j)=phi(( i−1)∗Nr+j,1); % valores de phi en el instante t 125 Tmatrep(i, j)=1/log(lambda)∗(−log(rmatrep(i , j)))+phimatrep(i , j); % valores de T en el instante t 126 velocidad_rmatrep(i , j)=velocidad_r (( i−1)∗Nr+j,1); 127 end 128 end 129 % Representaciones 130 xp(1:Ntheta ,1: Nr)=rmatrep(1:Ntheta ,1: Nr).∗cos( thetamatrep (1:Ntheta ,1: Nr)); 131 yp(1:Ntheta ,1: Nr)=rmatrep(1:Ntheta ,1: Nr).∗sin ( thetamatrep (1: Ntheta ,1: Nr)); 132 subplot (2,2,1) 133 plot ( ri∗cos(0:0.01:2∗pi) , ri∗sin (0:0.01:2∗pi), 'black ','LineWidth',2) 134 hold on 135 plot (re∗cos(0:0.01:2∗pi) , re∗sin (0:0.01:2∗pi) , 'black ','LineWidth',2) 136 hold on 137 contour(xp, yp, psimatrep) 138 xlim([−lambda lambda]) 139 ylim([−lambda lambda]) 140 xlabel ( 'x^∗') 141 ylabel ( 'y^∗') 142 hold off 143 colorbar 144 title ( 'Función de corriente \Psi^∗') 145 subplot (2,2,2) 146 plot ( ri∗cos(0:0.01:2∗pi) , ri∗sin (0:0.01:2∗pi), 'black ','LineWidth',2) 147 hold on 148 plot (re∗cos(0:0.01:2∗pi) , re∗sin (0:0.01:2∗pi) , 'black ','LineWidth',2) 149 hold on 150 contour(xp, yp, Tmatrep) 151 xlim([−lambda lambda]) 152 ylim([−lambda lambda]) 153 xlabel ( 'x^∗') 154 ylabel ( 'y^∗') 155 hold off 156 colorbar 157 title ( 'Campo de temperatura T^∗') 158 subplot (2,2,3) 159 plot (Tmatrep(:,round(Nr/2)), thetals (:) , 'r') 160 xlim([−1 0]) 161 ylim([0 2∗pi]) 3.4 Implementación numérica de las ecuaciones de Saltzman adimensionales 57 162 xlabel ( 'T^∗') 163 ylabel ( '\ theta^∗') 164 title ( 'Temperatura T^∗en r^∗/2 ') 165 hold off 166 subplot (2,2,4) 167 plot (velocidad_rmatrep (:, round(Nr/2)), thetals (:) , 'r') 168 ylim([0 2∗pi]) 169 xlim([−1 1]) 170 ylabel ( '\ theta^∗') 171 xlabel ( 'v_r^∗') 172 title ( 'Velocidad v_r^∗en r^∗/2') 173 hold off 174 pause(0.000001) 175 end 176 % Actualizaci ón de variables 177 psinm1=psi; 178 phinm1=phi; 179 end 64 Capítulo 4. Convección de Rayleigh-Bénard en una corona circular el instante de transición de 4 rollos a 8 rollos es en torno a t∗=50. También es perceptible que, a diferencia de la situación con Ra =1000 , en este caso, la función de corriente va aumentando su orden de magnitud a medida que trascurre el tiempo: su orden de magnitud para los instantes entre t∗=25 y t∗=70 es de 10−5 , y para los instantes posteriores se estabiliza en un orden de magnitud de 10−2. A partir de la evolución cronológica del campo de temperatura, se deduce que en un instante entre t∗=75 y t∗=500 se produce el salto de la convección en la corona, pues mientras que en t∗=75 las líneas isotermas son concéntricas, en t∗=500 ya no lo son. Así, a continuación, se analizará el campo de temperaturas en el radio medio de la corona en el intervalo de tiempo [75,500]para determinar con mayor exactitud el instante en el que salta la convección. (a) Solución para t∗=75. (b) Solución para t∗=100. (c) Solución para t∗=125. Figura 4.7 Evolución temporal de la temperatura T∗ en el radio medio de la corona para Ra =2500 cuando λ=2(I). De la figura (4.7) se deduce que el instante en el que salta la convección es, aproximadamente, t∗=100 , puesto que en los instantes previos, la desviación de la gráfica de la temperatura en el radio medio con respecto a una recta vertical (que es la forma de la gráfica característica de la situación en la que la convección no se ha producido) es prácticamente imperceptible. Analizando los instantes de tiempo posteriores, se observa que la gráfica del campo de temperatura en el radio medio de la corona va curvándose cada vez más, adoptando una forma periódica, hasta que en un determinado instante adquiere un cierto patrón, momento en el que se alcanza el régimen permanente. (a) Solución para t∗=500. (b) Solución para t∗=1000. (c) Solución para t∗=3000. Figura 4.8 Evolución temporal de la temperatura T∗ en el radio medio de la corona para Ra =2500 cuando λ=2(II). 4.1 Número de Rayleigh crítico 65 En efecto, como el campo de temperaturas en el instante t∗=500 , que está representado en la figura (4.6d), es prácticamente idéntico al asociado al instante t∗=1000 (y, en definitiva, al de los instantes posteriores), se deduce que para t∗=500 el fluido ya se encuentra en régimen permanente. Este fenómeno es también perceptible a partir de las gráficas de la temperatura en el radio medio de la corona, como se muestra a continuación. Según los resultados obtenidos hasta el momento, tras una pequeña perturbación inicial, cuando Ra =1000 , el movimiento del fluido se va amortiguando a medida que transcurre el tiempo hasta que finalmente regresa al estado de reposo inicial. Sin embargo, cuando Ra =2500 , transcurrido un determinado intervalo de tiempo en el que los valores de la función de corriente van aumentando, se produce la convección y, posteriormente, el fluido alcanza un régimen estacionario siguiendo un determinado patrón. De este modo, existirá un valor de Ra mínimo para el cual salta la convección en el fluido, alcanzando posteriormente el régimen estacionario y, por tanto, no regresando al estado inicial de reposo. Como se mencionó previamente, dicho valor es el denominado Rayleigh crítico, denotado como Rac. Así, para determinar el Rac , se analizarán las gráficas de función de corriente y campo de temperatura para distintos números de Ra en un intervalo de tiempo suficientemente grande que permita analizar la evolución del fluido y, por tanto, determinar si se produce convección natural en el mismo. A continuación, se recopilan los distintos valores de Ra que se han analizado para la determinación del Rac . Además, se incluyen los instantes de tiempo, aproximados, en los que salta la convección, siempre y cuando el fluido no tienda a amortiguarse y, por tanto, a volver a su situación de reposo inicial. Se recuerda que, tras los análisis llevados a cabo en este apartado, el Racestará dentro del intervalo [1000,2500]. Tabla 4.1 Búsqueda del Rayleigh crítico (Rac)para λ=2(I). Número de Rayleigh, Ra ¿Se produce la convección? Instante de tiempo, t∗ 2000 Sí 400 1900 Sí 1060 1800 No, el movimiento es amortiguado - En la tabla anterior (4.1) se puede observar que a medida que se va disminuyendo el Ra , el tiempo en el que salta la convección es cada vez mayor, hasta que finalmente, para Ra =1800 , ya no se produce la convección, si no que el movimiento se va amortiguando hasta que finalmente el fluido vuelve a su estado inicial de reposo. De este modo, se concluye que el Rac estará dentro del intervalo [1800,1900] . Así, se ha conseguido acotar el intervalo de búsqueda del Rayleigh crítico. Tabla 4.2 Búsqueda del Rayleigh crítico (Rac)para λ=2(II). Número de Rayleigh (Ra) ¿Se produce la convección? Instante de tiempo t∗ 1830 No, el movimiento es amortiguando - 1831 No, el movimiento es amortiguando - 1832 No, el movimiento es amortiguando - 1833 No, el movimiento es amortiguando - 1834 Sí 61000 66 Capítulo 4. Convección de Rayleigh-Bénard en una corona circular Del conjunto de números de Rayleigh que aparecen en la tabla (4.2), es destacable que a medida que se toma un Ra de mayor valor, el amortiguamiento del movimiento en la cavidad es cada vez más lento, hasta que finalmente, para un determinado Ra , el movimiento deja de amortigurse y se produce la convección. El primer valor de Ra , mayor que 1800 y menor que 1900 , para el que se da esta situación es 1834 , por lo que para una relación de aspecto λ=2 , el número de Rayleigh crítico es Rac=1834. A cotinuación, se presentan las gráficas asociadas al número de Rayleigh crítico ( Ra =1834 ) en el régimen estacionario. Figura 4.9 Solución para λ=2,Ra =1834 e instante de tiempo t∗=100000. En la figura anterior (4.9) se observa el hecho de que para el número de Rayleigh crítico, el régimen de convección natural alcanzado en el fluido es muy débil. 4.2 Tamaño de celdas según el número de Rayleigh En este apartado, se evaluará la variación del tamaño de las celdas convectivas con el número de Rayleigh en condiciones supercríticas, es decir, para Ra >Rac=1834. Para ello, se analizarán el número de rollos que aparecen en el fluido cuando se alcanza el régimen permanente para diferentes valores de Ra (mayores que el valor crítico, calculado en el apartado anterior). Los resultados se muestran en la tabla (4.3). 4.3 Relación entre el número de Rayleigh crítico y la relación de aspecto 67 Tabla 4.3 Número de rollos en régimen estacionario al variar el número de Rayleigh para λ=2 . Ra 1834 2000 2100 2200 3000 5000 Número de rollos 10 10 8 8 8 8 La anterior tabla muestra que para valores de Rayleigh ligeramente superiores a el crítico, el número de rollos en la corona circular disminuye de 10 a 8 . Cabe destacar que una disminución del número de rollos está estrechamente relacionada con un aumento del tamaño de las celdas de convección. Así, para regímenes moderadamente supercríticos, las celdas convectivas tienden a aumentar su tamaño debido al incremento en la intensidad del flujo convectivo que provoca la fusión de celdas más pequeñas en estructuras de mayor tamaño, organizándose en patrones más grandes, estables y ordenados. No obstante, en regímenes altamente supercríticos ( Ra >>> Rac=1834 ), el flujo puede volverse turbulento provocando la fragmentación del mismo en estructuras convectivas más pequeñas y caóticas. En efecto, existirá un tamaño máximo de las celdas a partir del cual el flujo se vuelva turbulento y las celdas disminuyan su tamaño, aumentando de nuevo el número de rollos. Las figuras (4.10a) y (4.10b) muestran las funciones de corrientes para Ra =2000 y Ra = 5000, respectivamente. (a) Ra =2000 y10 rollos (b) Ra =5000 y8rollos Figura 4.10 Función de corriente para λ=2 cuando Ra =2000 y Ra =5000 , respectivamente. 4.3 Relación entre el número de Rayleigh crítico y la relación de aspecto A continuación, se analizará la relación existente entre el número de Rayleigh crítico y la geometría del reciento, que está definida por la relación de aspecto λ . De este modo, se determinarán los números de Rayleigh críticos asociados a diferentes valores de λ mediante un procedimiento de prueba y error - análogamente al apartado anterior donde se obtuvo el Rac para λ=2 . Para llevar a cabo este estudio, resulta muy conveniente recuperar los resultados obtenidos por Lord Rayleigh, quien concluyó que, dado un problema de una capa de fluido horizontal- 68 Capítulo 4. Convección de Rayleigh-Bénard en una corona circular mente infinita, existe una diferencia de temperatura crítica entre las dos bases del recinto para la cual tiene lugar la convección natural. Obtuvo que dicho valor crítico se correspondía a un Rac=1708 , independientemente del fluido bajo estudio (es decir, válido tanto para líquido como para gas). Es importante destacar que dicho valor obtenido por Lord Rayleigh está referido a un número de Rayleigh basado en el espesor de la capa fluida, es decir, en la diferencia de alturas entre las bases - que en el caso de una corona circular es la diferencia entre el radio exterior y el radio interior: H=Re−Ri . Sin embargo, el Ra utilizado en el presente trabajo está basado en el radio interior de la corona. Por tanto, se establecerá una relación entre ambos números de Rayleigh basados en diferentes longitudes características para posteriores comparaciones de resultados: RaH=gβ∆T(Re−Ri)3 να RaRi=gβ∆T R3 i να    RaH=RaRiRe −Ri Ri 3 =RaRiRe Ri−13 ,(4.1) siendo RaRi el número de Rayleigh, basado en el radio interior de la corona, que se ha implementado en el programa de Matlab para la resolución del problema - que hasta el momento se ha detonado simplemente como Ra. Se recuerda que la relación de aspecto se ha definido como λ=Re Ri =r∗ e r∗ i , de forma que las expresiones del radio exterior e interior de la corona circular en coordenadas adimensionales se pueden expresar respectivamente como r∗ i=1yr∗ e=λ.(4.2) La relación de aspecto deberá ser siempre mayor que 1; en el caso de que λ=1, el recinto sería una corona sin área, es decir, un círculo. Así, cuando λ esté muy próximo a 1, el fluido dentro de la corona se asemejará a una capa horizontalmente infinita y se comprobará si su RaHes próximo al resultado obtenido por Lord Rayleigh. A continuación, se presentan diferentes geometrías del recinto que se obtendrán variando la relación de aspecto. En relación a las expresiones (4.2), el radio interior de la corona se mantiene fijo y es el exterior el que cambia al variar λ. 4.3 Relación entre el número de Rayleigh crítico y la relación de aspecto 69 (a) Geometría para λ=1.2(b) Geometría para λ=2(c) Geometría para λ=3 Figura 4.11 Geometrías del recinto variando la relación de aspecto λ. Introduciendo la expresión de la relación de aspecto en (4.1) se obtiene: RaH=RaRi(λ−1)3.(4.3) Los resultados obtenidos de los números de Rayleigh críticos - basados en el radio interior de la corona - para distintas geometrías del recinto se presentan en las siguientes dos tablas: Tabla 4.4 Números de Rayleigh críticos RaRicen función de la geometría del problema (I). Relación de aspecto (λ)1.25 1.6 1.7 1.8 2 2.5 3 Rayleigh crítico RaRic108000 7889 5225 3516 1834 573 247 Tabla 4.5 Números de Rayleigh críticos RaRicen función de la geometría del problema (II). Relación de aspecto (λ)3.5 4 5 6 7 8 9 Rayleigh crítico RaRic134 84 36 20 12 8 6 Para el análisis de la relación entre el número de Rayleigh crítico y la relación de aspecto se han empleado 14 geometrías diferentes. Ambas tablas (4.4) y (4.5) demuestran que a medida que la relación de aspecto va disminuyendo, lo que se traduce en una distancia menor entre el radio interior r∗ i y exterior r∗ e de la corona (es decir, un menor espesor de la capa fluida), el número de Rayleigh critico RaRicva aumentando significativamente. 70 Capítulo 4. Convección de Rayleigh-Bénard en una corona circular Figura 4.12 Variación de RaRicen función de λ. Empleando la relación (4.3) se pueden obtener los valores de Rayleigh críticos referidos a H a partir de los valores de Rayleigh críticos referidos al radio interior de la corona, RaRic , y las relaciones de aspecto λ correspondientes a cada caso. Estos valores se muestran en las siguientes dos tablas. Además, dichas tablas también recogen el número de rollos que aparecen en las gráficas de la función de corriente para cada geometría. Tabla 4.6 Números de Rayleigh críticos RaHc y número de rollos en función de la geometría del problema (I). Relación de aspecto (λ)1.25 1.6 1.7 1.8 2 2.5 3 Rayleigh crítico RaHc1686 1704 1792 1800 1834 1934 1976 Nºrollos 24 14 12 10 10 8 6 Tabla 4.7 Números de Rayleigh críticos RaHc y número de rollos en función de la geometría del problema (II). Relación de aspecto (λ)3.5 4 5 6 7 8 9 Rayleigh crítico RaHc2094 2268 2304 2500 2592 2744 3072 Nºrollos 6644444 A partir de las tablas (4.6) y (4.7) se aprecia que cuando la relación de aspecto se acerca a valores próximos a la unidad, los valores de RaHc tienden a valores cercanos a 1708 , que es el valor que obtuvo Lord Rayleigh experimentalmente para la convección entre dos placas planas horizontales infinitas. En efecto, cuando una corona circular tiene una diferencia entre el diámetro exterior e interior muy pequeña, es decir, para los casos en los que la relación de aspecto λ es próxima a 1, el espacio entre las dos superficies cilíndricas es pequeño en comparación con los diámetros y la cavidad se puede aproximar a una disposición semejante a dos placas planas horizontales infinitas. De esta forma, para casos de λ próximo a 1, la curvatura 4.3 Relación entre el número de Rayleigh crítico y la relación de aspecto 71 de las superficies es mínima y el gradiente de temperatura se desarrolla en la dirección radial. Bajo estas condiciones, el comportamiento del flujo y la transferencia de calor en la corona se asemejan a los asociados a dos pacas planas infinitas dispuestas horizontalmente con respecto a la dirección de la gravedad. Esta misma tendencia es la que obtuvo Pablo Ruiz Contreras [ 2 ] para la convección de Rayleigh-Bénard en una cavidad rectangular con paredes verticales adiabáticas. Cabe destacar que en el presente trabajo la gravedad ha sido tomada en dirección radial. Esto resulta muy favorable a la hora de realizar la comparación entre la corona circular estrecha y las dos placas planas horizontales infinitas. Numerosos autores han analizado este caso límite obteniendo correlaciones entre el número de Rayleigh crítico y la relación de aspecto de la cavidad. Un ejemplo de ello es la siguiente relación propuesta por Walton en 1980 [6]: Rac∼1707.8+516.8 A,(4.4) siendo Ala relación de aspecto definida como: A=Di Re−Ri =2Ri Re−Ri =2 λ−1.(4.5) Es evidente que para λ próximo a 1, el denominador de la expresión (4.5) se vuelve 0 y, por tanto, A tiende a infinito. De esta forma, la relación de Walton (4.4) proporciona que Rac∼1707.8cuando λ∼1. Por otro lado, cuando la relación de aspecto aumenta, los valores de RaHc también aumentan ya que los recintos donde se aloja el fluido requieren de mayor diferencia de temperatura para que tenga lugar la convección al tener una capa fluida con mayor espesor H=r∗ e−r∗ i . Cabe destacar que la única relación de aspecto que cumple que sus valores RaRic y RaHc son iguales es λ=2. A continuación, se presentan las gráficas de función de corriente Ψ∗ y campo de temperatura T∗para varias relaciones de aspecto con sus respectivos números de Rayleigh críticos RaHc: Figura 4.13 Solución para λ=5con RaHc=2304. 72 Capítulo 4. Convección de Rayleigh-Bénard en una corona circular Figura 4.14 Solución para λ=3con RaHc=1976. Figura 4.15 Solución para λ=2.5con RaHc=1934. Figura 4.16 Solución para λ=1.8con RaHc=1800. 4.4 Número de Nusselt 73 Figura 4.17 Solución para λ=1.7con RaHc=1792. Aparte de las observaciones previamente mencionadas, cabe destacar que el número de rollos en el régimen estacionario en situaciones críticas (es decir, cuando Ra =Rac ) es mayor cuando las relaciones de aspecto son menores, es decir, para coronas circulares más estrechas. 4.4 Número de Nusselt El número de Nusselt ( Nu ) se define como un parámetro adimensional que representa la relación entre la transferencia de calor por convección (es decir, cuando existe movimiento convectivo en el fluido) y la transferencia de calor por conducción (es decir, cuando el fluido se encuentra en equilibrio hidrostático) a través de una superficie en un sistema fluido. Así: Nu =Transferencia de calor por convección Transferencia de calor por conducción =qpconvecci ´on qpconducci ´on .(4.6) A continuación, se obtendrá la expresión del número de Nusselt en la superficie interior de la corona con r∗=r∗ i. Para ello, es necesario definir el calor aportado por una pared de coordenada r al fluido en dirección radial, cuya expresión viene dada por la ley de Fourier: qp=−k∂T ∂rr .(4.7) Realizando el cambio de variables para implementar la adimensionalización desarrollada en los apartados anteriores, se obtiene la siguiente expresión: qp=−k∂T ∂rr =−k∆T Ri ∂T∗ ∂r∗r∗ =−k∆T Ri ∂ ∂r∗−log(r∗) log(λ)+ϕ∗r∗ (4.8) La expresión anterior engloba la transmisión de calor por conducción y por convección al tener en cuenta ambas componentes de la temperatura adimensional, es decir, la componente de la temperatura hidrostática y la componente de la temperatura asociada al movimiento del fluido. Por tanto, comparando este calor con el que se transferiría al fluido si estuviese en reposo, se obtiene la siguiente expresión, que coincide con el número de Nusselt (Nu): 80 Capítulo 5. Conclusiones y líneas futuras iniciales. La resolución de las ecuaciones permite visualizar la evolución de la función de corriente y el campo de temperatura en la corona circular desde el régimen inicial transitorio hasta alcanzar un régimen estacionario. De esta forma, se puede determinar el número de Rayleigh crítico y el número de rollos asociados a cada relación de aspecto. Mediante el número de Rayleigh crítico, que estará asociado a una relación de aspecto concreta, se puede determinar el valor del número de Rayleigh a partir del cual el proceso térmico que tiene lugar en el fluido no está dominado por la conducción si no por el fenómeno de la convección. Al comparar estos valores de Rayleigh críticos con otros existentes en la literatura para cavidad rectangulares, se observa que la evolución de la tendencia para ambas configuraciones geométricas es similar: a medida que el espesor de la capa fluido se hace más pequeño, el número de Rayleigh crítico va tendiendo al resultado experimental para placas infinitas, cuyo valor es aproximadamente 1708. Cabe destacar que como en el presente trabajo se ha tomado la gravedad en dirección radial (perpendicular a las parades de la corona circular), no es descabellado asumir de antemano que los resultados numéricos obtenidos para geometría anular estrecha y los resultados experimentales obtenidos para las placas horizontales infinitas deben ser similares. Adicionalmente, se observa que, para una relación de aspecto fija, al aumentar el número de Rayleigh, el número de rollos en régimen estacionario disminuye, lo que lleva a un aumento del tamaño de las celdas de Bénard. Esto está estrechamente relacionado con el hecho de que a mayor número de Rayleigh, mayor será la intensidad del flujo convectivo y esto llevará a la fusión de celdas más pequeñas en estructuras de mayor tamaño, originándose estructuras más ordenadas y estables. El número de Nusselt, que representa la relación entre la transferencia de calor por convección y por conducción, permite obtener el tiempo que el fluido tarda en alcanzar el régimen estacionario. Mediante la comparación de las evoluciones temporales del Nusselt medio para una misma relación de aspecto y diferentes valores de Rayleigh, se obtiene que cuanto mayor es el número de Rayleigh menos tiempo tarda el fluido en alcazar el régimen permanente y mayor es el valor del Nusselt medio asociado a dicho régimen. Para validar este último hecho se emplean los resultados experimentales obtenidos por Silverstone en 1958 y por Koschmieder y Pallas en 1974. 5.2 Líneas futuras Una línea de mejora del presente trabajo podría consistir en la resolución del problema mediante un programa de CFD, como ANSYS, que permita llevar a cabo simulaciones más detallas y complejas con un mayor rendimiento, aumentando la precisión de los resultados. Así, los resultados obtenidos mediante ANSYS se podrían comparar con los obtenidos mediante MATLAB. Esto serviría como método de verificación y validación de los resultados. Además, el empleo de ANSYS facilitaría el estudio de situaciones de convección natural turbulenta. En cuanto al método numérico de resolución, se emplea el método de colocación con interpolantes de Lagrange que aproximan a las funciones incógnitas en los puntos de colocación, tomándose los denominados nodos de Chebyshev para la coordenada r y nodos equiespaciados para la coordenada θ . En el presente trabajo, se ha detallado la importancia que tiene la selección de estos nodos en la calidad de la interpolación de Lagrange de una función al analizar el fenómeno de Runge. No obstante, aunque a priori este método numérico propociona muy buenos resultados, se propone, como línea futura de mejora, analizar cómo varían los 5.2 Líneas futuras 81 resultados de la convección de Rayleigh-Bénard al emplear otro método numérico. Por otro lado, se ha observado que si el número de Rayleigh es muy próximo al crítico, la convección que se produce es muy débil y el tiempo en alcanzar el régimen estacionario es elevado. Por tanto, la precisión para determinar el Rayleigh crítico podría verse afectada. Así, se propone mejorar dicha precisión como línea futura del trabajo. Finalmente, se propone estudiar el, aparentemenete no muy estudiado hasta ahora, efecto sobre la convección natural del movimiento de una de las paredes del recinto, lo que introduce una componente de convección forzada adicional. Índice de Figuras 1.1 Celdas originadas debido a la convección de Rayleigh-Bénard 2 1.2 Ejemplos de convección de Rayleigh-Bénard en la naturaleza 3 1.3 Patrón de movimiento de convección desarrollado en una capa de fluido al calentarla desde abajo 3 1.4 Vista lateral de las celdas de Bénard 3 2.1 Corona circular en coordenadas ryθ8 3.1 Interpolación con nodos equiespaciados y con nodos de Chebyshev para N=720 3.2 Interpolación con nodos equiespaciados y con nodos de Chebyshev para N=15 21 3.3 Representación de un five-point stencil en una y dos dimensiones, respectivamente 22 3.4 Representación del mallado bidimensional en coordenadas xyz29 3.5 Ejemplo de mallado bidimensional con Nx=6yNz=330 3.6 Representación del mallado bidimensional en coordenadas θyr32 3.7 Índices de los nodos donde se aplican las condiciones de contorno del problema 34 3.8 Solución numérica de la ecuación de Poisson en una corona circular 37 3.9 Solución analítica de la ecuación de Poisson en una corona circular 37 3.10 Solución numérica de la ecuación de Laplace no estacionaria para t=0.001 42 3.11 Solución numérica de la ecuación de Laplace no estacionaria para t=0.06 43 3.12 Solución numérica de la ecuación de Laplace no estacionaria para t=0.53 43 3.13 Solución numérica de la ecuación de Laplace no estacionaria para t=1.63 43 3.14 Solución numérica de la ecuación de Laplace no estacionaria para t=3.644 3.15 Solución numérica de la ecuación de Laplace no estacionaria para t=4.444 3.16 Corona circular en coordenadas adimensionales r∗yθ∗45 3.17 Representación del mallado bidimensional en coordenadas adimensionales θ∗yr∗46 3.18 Índices de los nodos donde se aplican las condiciones de contorno 53 3.19 Solución de las ecuaciones de Saltzman adimensionales en una corona circular con λ=2para t∗=3000 58 4.1 Solución para λ=2,Ra =1000 e instante de tiempo t∗=50 60 4.2 Solución para λ=2,Ra =1000 e instante de tiempo t∗=500 60 4.3 Solución para λ=2,Ra =2500 e instante de tiempo t∗=50 61 4.4 Solución para λ=2,Ra =2500 e instante de tiempo t∗=500 62 4.5 Evolución temporal de la función de corrienteΨ∗para Ra =2500 cuando λ=263 4.6 Evolución temporal del campo de temperatura T∗para Ra =2500 cuando λ=263 83 84 Índice de Figuras 4.7 Evolución temporal de la temperatura T∗ en el radio medio de la corona para Ra =2500 cuando λ=2(I) 64 4.8 Evolución temporal de la temperatura T∗ en el radio medio de la corona para Ra =2500 cuando λ=2(II) 64 4.9 Solución para λ=2,Ra =1834 e instante de tiempo t∗=100000 66 4.10 Función de corriente para λ=2cuando Ra =2000 yRa =5000, respectivamente 67 4.11 Geometrías del recinto variando la relación de aspecto λ69 4.12 Variación de RaRicen función de λ70 4.13 Solución para λ=5con RaHc=2304 71 4.14 Solución para λ=3con RaHc=1976 72 4.15 Solución para λ=2.5con RaHc=1934 72 4.16 Solución para λ=1.8con RaHc=1800 72 4.17 Solución para λ=1.7con RaHc=1792 73 4.18 Evolución temporal del Nusselt medio Nu para λ=2yRa =2150 75 4.19 Evolución temporal del Nusselt medio Nu para λ=2yRa =1900 75 4.20 Resultados experimentales obtenidos por Silverstone (1958) 76 4.21 Resultados experimentales obtenidos por Koschmieder y Pallas (1974) 77 4.22 Nusselt medio Nu en función del número de Rayleigh para λ=277 Índice de Tablas 2.1 Aproximación de valores (ρ,β,k) para agua y aire 9 4.1 Búsqueda del Rayleigh crítico (Rac)para λ=2(I) 65 4.2 Búsqueda del Rayleigh crítico (Rac)para λ=2(II) 65 4.3 Número de rollos en régimen estacionario al variar el número de Rayleigh para λ=267 4.4 Números de Rayleigh críticos RaRicen función de la geometría del problema (I) 69 4.5 Números de Rayleigh críticos RaRicen función de la geometría del problema (II) 69 4.6 Números de Rayleigh críticos RaHc y número de rollos en función de la geometría del problema (I) 70 4.7 Números de Rayleigh críticos RaHc y número de rollos en función de la geometría del problema (II) 70 4.8 Tiempos adimensionales y valores de Nu asociados al régimen permanente para λ=2y diferentes Ra >Rac76 4.9 Nusselt medio Nu en función del número de Rayleigh para λ=277 85 Índice de Códigos 3.1 Cálculo de matrices de diferenciación para un stencil de M puntos en un punto ¯xi del mismo 25 3.2 Cálculo de matrices de diferenciación para una malla de N nodos en un punto xi de la misma 27 3.3 Resolución numérica y analítica de la ecuación de Poisson en una corona circular 35 3.4 Resolución numérica de la ecuación de Laplace en régimen no estacionario en una corona circular 39 3.5 Resolución numérica de las ecuaciones de Saltzman adimensionales en una corona circular 53 87 Bibliografía [1] Jasim Abdulateef, Effect of prandtl number and diameter ratio on laminar natural convection from isothermal horizontal cylindrical annuli, Diyala Journal of Engineering Sciences 9(2016), 39–54. [2] Pablo José Ruiz Contreras, El método de colocación para el problema de convección de rayleigh-bénard, Proyecto fin de carrera, Universidad de Sevilla, 2013. [3] Belén Jurado de la Cuadra, Estudio numérico de la convección de rayleigh-bénard y de flujos giratorios en cavidades cilíndricas, Trabajo fin de carrera, Universidad de Sevilla, 2022. [4] Eduardo M. García Juárez, Procesos de convección natural con hipótesis anelástica, Proyecto fin de carrera, Universidad de Sevilla, 2014. [5] T. H. Kuehn and R. J. Goldstein, An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, Journal of Fluid Mechanics 74 (1976), no. 4, 695–719. [6] Jiro Mizushima, Sachiko Hayashi, and Takahiro Adachi, Transitions of natural convection in a horizontal annulus, International Journal of Heat and Mass Transfer 44 (2001), no. 6, 1249–1257. [7] G. Petrone, E. Chénier, and G. Lauriat, Stability of free convection in air-filled horizontal annuli: influence of the radius ratio, International Journal of Heat and Mass Transfer 47 (2004), no. 17, 3889–3907. [8] Miguel Pérez-Saborid, A simple matlab program to compute differentiation matrices for arbitrary meshes via lagrangian interpolation, (2019). [9] Lord Rayleigh, Lix. on convection currents in a horizontal layer of fluid, when the higher temperature is on the under side, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 32 (1916), no. 192, 529–546. [10] Rao Yan-Fei, Miki Yasutomi, Fukuda Kenji, Takata Yasuyuki, and Hasegawa Shu, Flow patterns of natural convection in horizontal cylindrical annuli, International Journal of Heat and Mass Transfer 28 (1985), no. 3, 705–714. [11] Xiufeng Yang and Song-Charng Kong, Numerical study of natural convection in a horizontal concentric annulus using smoothed particle hydrodynamics, Engineering Analysis with Boundary Elements 102 (2019), 11–20. 89