TFM-2418 Torne Alaminos

Vista previa del material en texto

Equation Chapter 1 Section 1 
Trabajo Fin de Máster 
Máster Universitario en Ingeniería Aeronáutica 
 
La hipótesis anelástica en problemas de convección 
de Rayleigh-Bénard para cavidades rectangulares: 
aspectos termo-fluidomecánicos y energéticos 
Autor: Pablo Torné Alaminos 
Tutor: Miguel Pérez-Saborid Sánchez-Pastor 
Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos 
Escuela Técnica Superior de Ingeniería 
Universidad de Sevilla 
 Sevilla, 2022 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
iii 
 
Trabajo Fin de Máster 
Máster Universitario en Ingeniería Aeronáutica 
 
 
 
 
 
La hipótesis anelástica en problemas de convección 
de Rayleigh-Bénard para cavidades rectangulares: 
aspectos termo-fluidomecánicos y energéticos 
 
 
Autor: 
Pablo Torné Alaminos 
 
 
Tutor: 
Miguel Pérez-Saborid Sánchez-Pastor 
Profesor Titular 
 
 
 
Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos 
Escuela Técnica Superior de Ingeniería 
Universidad de Sevilla 
Sevilla, 2022 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
v 
 
 
Trabajo Fin de Máster: La hipótesis anelástica en problemas de convección de Rayleigh-Bénard para 
cavidades rectangulares: aspectos termo-fluidomecánicos y energéticos 
 
 
 
 
 
Autor: Pablo Torné Alaminos 
Tutor: Miguel Pérez-Saborid Sánchez-Pastor 
 
 
El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: 
Presidente: 
 
 
 
Vocales: 
 
 
 
 
Secretario: 
 
 
 
 
Acuerdan otorgarle la calificación de: 
 
 
 
 Sevilla, 2022 
 El Secretario del Tribunal 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
vii 
 
 
Agradecimientos 
Este proyecto es el culmen a mi formación como Ingeniero Aeronáutico. Durante toda esta etapa, compuesta 
por el Grado en Ingeniería Aeroespacial y este máster, se pasan por momentos de suma felicidad, pero también 
de gran agotamiento y dificultad. Sin las personas que me han acompañado durante este camino, sé que todo 
habría sido infinitamente más complicado. Gracias a mis padres, a mi abuela, a mis hermanos, a mis amigos de 
la carrera, a los que he conocido fuera de ella, y a mis amigos que siempre me han acompañado. Cada momento 
compartido a lo largo de estos años con ellos ha sido siempre una ayuda para lograr mis objetivos. 
Agradecer a Miguel, tutor de este proyecto, tanta dedicación, y compromiso durante su realización. La Escuela 
Técnica Superior de Ingeniería de Sevilla, y en especial todos los estudiantes de Ingeniería Aeronáutica, tienen 
la gran suerte de contar con un docente como él. 
La carga de conocimientos adquiridos durante todo este camino es muy extensa, e incluso complicada de 
recordar a veces. Sin embargo, tengo claro lo que caracteriza a cada persona que ha conseguido superar esta 
carrera: se podrá no recordar un concepto concreto de un tema, ni tener la respuesta inmediata a un problema, 
pero mientras exista una solución, un ingeniero nunca bajará los brazos hasta conseguirla. 
 
Pablo Torné Alaminos 
Sevilla, 2022 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
ix 
 
 
Resumen 
Este proyecto consiste en proponer un modelo que aproxime y simule el comportamiento y la dinámica de un 
fluido sometido a un gradiente de temperaturas en un recinto bidimensional. El avance de este proyecto frente a 
otros realizados en el departamento es la inclusión de la densidad como variable del problema; es decir, se 
considera la estratificación de la densidad. Para ello, se emplea la hipótesis anelástica, formulando las ecuaciones 
del modelo en base a sus simplificaciones y exigencias. Además, se consideran polítropos para modelar el estado 
de referencia. 
Las condiciones de contorno más representativas son la fijación de temperaturas en la base y en la superficie 
superior, siendo 𝑇𝑖𝑛𝑓 > 𝑇𝑠𝑢𝑝 . En la resolución numérica, se emplea el método de diferencias finitas para la 
obtención de matrices de derivación, un mallado siguiendo nodos de Chebyshev, y una de discretización que 
permite resolver el problema de manera semi-implícita. 
Para poner a prueba el modelo, y entender la física del problema, se analizan una serie de resultados, comparando 
resultados de la hipótesis anelástica con los de la aproximación de Boussinesq, un modelo que no permite la 
estratificación de la densidad, empleado en proyectos anteriores. A su vez, se estudia la influencia de los 
parámetros característicos sobre el modelo, concluyendo el análisis con un balance energético. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xi 
 
 
Abstract 
This project consists of proposing a model that approximates and simulates the behavior and dynamics of a fluid 
subjected to a temperature gradient in a two-dimensional enclosure. The advance of this project compared to 
others carried out in the department is the inclusion of density as a variable of the problem; that is, density 
stratification is considered. For this purpose, the anelastic hypothesis is used, formulating the model equations 
based on its simplifications and requirements. In addition, polytropes are considered to model the reference state. 
The most representative boundary conditions are the setting of temperatures at the base and at the top surface, 
being 𝑇𝑏𝑜𝑡 > 𝑇𝑡𝑜𝑝. In the numerical resolution, the finite difference method is used to obtain derivation 
matrices, a meshing following Chebyshev nodes, and a discretization method that allows solving the problem in 
a semi-implicit way. 
To test the model, and to understand the physics of the problem, a series of results are analyzed, comparing 
results of the anelastic hypothesis with those of the Boussinesq approximation, an incompressible model used 
in previous projects. In turn, the influence of the characteristic parameters on the model is studied, concluding 
the analysis with an energy balance. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xiii 
 
 
Índice 
Agradecimientos vii 
Resumen ix 
Abstract xi 
Índice xiii 
Índice de Tablas xv 
Índice de Figuras xvii 
1 Introducción 1 
1.1. Convección de Rayleigh-Bénard 1 
1.2. Ejemplos de convección natural 3 
1.3. Objetivo del proyecto. Estructura. 4 
2 Formulación del problema 7 
2.1. Estabilidad a convección de un fluido 7 
2.2. Modelos de convección 10 
2.2.1. Aproximación de Boussinesq 11 
2.2.2. Modelo anelástico 12 
2.3. Resolución del modelo anelástico 22 
2.3.1. Ecuación de la vorticidad. Función de corriente 22 
2.3.2. Estado de referencia: Polítropos 24 
2.3.3. Condiciones de contorno. Condiciones iniciales 27 
2.3.4. Adimensionalización 28 
3 Método numérico de resolución 31 
3.1. Método de diferencias finitas 31 
3.1.1 Mallado de M puntos 31 
3.1.2 Matrices de derivación para una malla arbitraria 35 
3.2. Discretización de las ecuaciones del problema anelástico 37 
3.2.1. Discretización completa de las ecuaciones 39 
3.2.2. Aplicación de las condiciones de contorno sobre el mallado 42 
4 Análisis de los resultados 45 
4.1. Número de Rayleigh crítico (𝐑𝐚𝐜) 46 
4.1.1 Influencia del índice politrópico 𝒏 en el 𝐑𝐚𝐜 50 
4.1.2 Influencia del factor de escala𝑳/𝑯 en el 𝐑𝐚𝐜 54 
4.1.3 Influencia de la relación 𝚫𝑻/𝑻𝟏 en el 𝐑𝐚𝐜 58 
4.1.4 Evolución de las magnitudes de referencia, y los campos 𝑻𝒏, 𝚿𝒏 en función del tiempo. 59 
4.2. Número de celdas convectivas 62 
4.1.1 Influencia del número del factor de escala y del número de Rayleigh en el número de celdas 62 
4.1.2 Situaciones de interés 63 
4.3. Número adimensional de conducción, 𝐍𝐜𝐨𝐧𝐝 65 
4.3.1. Influencia del número de Rayleigh y del índice politrópico en 𝐍𝐜𝐨𝐧𝐝 70 
4.4. Transferencia de energía en el proceso 72 
 
 
 
4.4.1. Obtención de las expresiones energéticas 72 
4.4.2. Evolución de las componentes de la energía y flujos de calor para distintos casos. 77 
5 Conclusiones y líneas de desarrollo 83 
Anexos 85 
A: Código Matlab 86 
Referencias 101 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xv 
 
 
ÍNDICE DE TABLAS 
 
 
Tabla 4.1. Parámetros de simulación de Figura 4.1 . 47 
Tabla 4.2. Parámetros de la simulación de Figura 4.2. 47 
Tabla 4.3. Parámetros de simulación de la Figura 4.3. 49 
Tabla 4.4. Parámetros de simulación de Figura 4.4. 49 
Tabla 4.5. Parámetros de simulación de la Figura 4.5 50 
Tabla 4.6. Parámetros de simulación de la Figura 4.6 50 
Tabla 4.7. Parámetros de simulación de la Figura 4.7. 51 
Tabla 4.8. Parámetros de simulación de la Figura 4.8 52 
Tabla 4.9. Dependencia número de Rayleigh crítico con el índice politrópico para un gas con 𝛾 = 1.4 y 𝐿/𝐻 =
2 53 
Tabla 4.10. Parámetros para la simulación de la Figura 4.11. 54 
Tabla 4.11. Parámetros para la simulación de la Figura 4.12. 54 
Tabla 4.12. Parámetros para la simulación de la Figura 4.13. 55 
Tabla 4.13. Parámetros para la simulación de la Figura 4.14. 56 
Tabla 4.14. Dependencia del número de Rayleigh crítico con el factor de escala. 57 
Tabla 4.15. Valores del Rayleigh crítico en función de la relación Δ𝑇/𝑇1 para distintos índices politrópicos, 
con 𝐿/𝐻 = 2 y Pr = 0.733 58 
Tabla 4.16. Número de celdas convectivas en función del factor de escala, para dos casos distintos de índice 
politrópico, y Ra = Rac 63 
Tabla 4.17. Número de celdas convectivas en función del factor de escala, para dos casos distintos de índice 
politrópico, y Ra = 3Rac 63 
Tabla 4.18. Parámetros escogidos para la simulación de Figura 4.20 . 64 
Tabla 4.19. Parámetros escogidos para la simulación de la Figura 4.21. 64 
Tabla 4.20. Parámetros para la simulación de la Figura 4.23. 67 
Tabla 4.21. Parámetros para la simulación de 70 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
xvii 
 
 
ÍNDICE DE FIGURAS 
 
 
Figura 1.1. Diagrama de celdas hexagonales de Bénard [13] 2 
Figura 1.2. Celdas de Bénard en esperma. Reproducción de una de las fotografías originales de Bénard 2 
Figura 1.3 Esquema de los fenómenos involucrados en la convección de Rayleigh-Bénard. [13] 2 
Figura 1.4. Esquema de los mecanismos de transferencia de calor en el calentamiento de agua de una olla 3 
Figura 1.5. Modelo de circulación atmosférica propuesto por George Hadley [24] 4 
Figura 1.6. Modelo de circulación atmosférica de tres células convectivas [23] 4 
Figura 1.7. Capas del Sol [22] 5 
Figura 1.8. Convección en el manto terrestre [21] 5 
Figura 2.1. Estabilidad a convección térmica: estado subadiabático. Evolución de la temperatura y la densidad 
frente a la altura. 8 
Figura 2.2. Inestabilidad a convección térmica: estado superadiabático. Evolución de la temperatura y la 
densidad frente a la altura. 9 
Figura 2.3. Recinto prismático rectangular de modelado del volumen fluido. 10 
Figura 2.4. Evolución de las magnitudes de referencia adimensionalizadas frente a la altura, con Δ𝑇 =
100𝐾, 𝑇1 = 303 𝐾: temperatura hidrostática (izquierda) y densidad hidrostática (derecha). 26 
Figura 2.5. Recinto rectangular de modelado del fluido 27 
Figura 3.1. Sentido de numeración de los nodos del mallado 38 
Figura 3.2. Ubicación de los nodos del contorno y subcontorno del recinto 43 
Figura 4.1. Función de corriente (izquierda) y campo de temperaturas (derecha) para Ra < Rac, asociados a los 
parámetros de la Tabla 4.1. 47 
Figura 4.2. Función de corriente y campo de temperaturas para Ra > Rac, asociados a los parámetros de la 
Tabla 4.2. 48 
Figura 4.3. Función de corriente y campo de temperaturas para Ra ≈ Rac, asociados a los parámetros de la 
Tabla 4.3.. 49 
Figura 4.4. Función de corriente y campo de temperaturas para Ra ≫ Rac, asociados a los parámetros de la 
Tabla 4.4. 49 
Figura 4.5. Función de corriente y campo de temperaturas para 𝑛 = 1.5, Ra ≈ Rac, asociados a los parámetros 
de la Tabla 4.5. 50 
Figura 4.6. Función de corriente y campo de temperaturas para 𝑛 = 2, Ra ≈ Rac, asociados a los parámetros 
de la Tabla 4.6. 51 
Figura 4.7. Función de corriente y campo de temperaturas para 𝑛 = 2.4, Ra ≈ Rac, asociados a los parámetros 
de la Tabla 4.7. 51 
Figura 4.8. Función de corriente y campo de temperaturas para 𝑛 = 1.5, Ra ≫ Rac, asociados a los parámetros 
de la Tabla 4.8. 52 
Figura 4.9. Dependencia del número de Rayleigh crítico con el índice politrópico para un gas con 𝛾 = 1.4
 52 
 
 
 
Figura 4.10. Evolución del Rac en función de 𝑛, para un gas con 𝛾 = 1.4, obtenida en [5]. 53 
Figura 4.11. Función de corriente y campo de temperaturas para 𝐿/𝐻 = 3, Ra ≈ Rac, asociados a los 
parámetros de la Tabla 4.10. 54 
Figura 4.12. Función de corriente y campo de temperaturas para 𝐿/𝐻 = 5, Ra ≈ Rac, asociados a los 
parámetros de la Tabla 4.11. 55 
Figura 4.13. Función de corriente y campo de temperaturas para 𝐿/𝐻 = 10, Ra ≈ Rac, asociados a los 
parámetros de la Tabla 4.12. 55 
Figura 4.14. Función de corriente y campo de temperaturas para 𝐿/𝐻 = 5, Ra ≫ Rac, asociados a los 
parámetros de la Tabla 4.13. 56 
Figura 4.15. Dependencia del número de Rayleigh crítico con el factor de escala para 𝑛 = 1.5 57 
Figura 4.16. Dependencia del número de Rayleigh crítico con el factor de escala para varios índices politrópicos
 58 
Figura 4.17. Representación del Rayleigh crítico en función de la relación Δ𝑇/𝑇1 para distintos índices 
politrópicos, con 𝐿/𝐻 = 2 y Pr = 0.733 59 
Figura 4.18. Evolución temporal para 𝑛 = 0 del campo de función de circulación, campo de temperaturas, 
perfil de la temperatura total e hidrostática en 𝑥 = 𝐿/2, y densidad total e hidrostática en 𝑥 = 𝐿/2. 60 
Figura 4.19. Evolución temporal para 𝑛 = 1.5 del campo de función de circulación, campo de temperaturas, 
perfil de la temperatura total e hidrostática en 𝑥 = 𝐿/2, y densidad total e hidrostática en 𝑥 = 𝐿/2. 61 
Figura 4.20. Evolución de los campos de función de corriente y temperatura correspondientes a los parámetros 
de la Tabla 4.18. en distintos instantes de tiempo. 64 
Figura 4.21. Evolución de los campos de función de corriente y temperatura correspondientes a los parámetros 
de la Tabla 4.19. en distintos instantes de tiempo 65 
Figura 4.22. Flujos de calor por la superficie horizontal z 66 
Figura 4.23. Representación del número de conducción (izquierda) y del perfil de perturbaciones de 
temperaturas en 𝑥 = 𝐿/2 (derecha) frente a 𝑧 ∗, según los parámetros de la Tabla 4.20. 68 
Figura 4.24. Evolución del número de conducción adimensional frente al tiempo, para los datos de la Tabla 
4.20 (Tabla 4.13 modificada en 𝑁𝑧 y 𝜏 ∗). 69 
Figura 4.25. Representación del número de conducción (izquierda) y evolución temporal del número de 
condición en 𝑧 = 0 y 𝑧 = 𝐻 (derecha), para 𝑛 = 0, 𝐿/𝐻 = 3 y Ra = 1𝑒4. 69 
Figura 4.26. Representación del número de conducción (izquierda) y del perfil de perturbaciones de 
temperaturas en 𝑥 = 𝐿/2 (derecha) frente a 𝑧 ∗ según los parámetros de la Tabla 4.21. 70 
Figura 4.27. Número de conducción en función del número de Rayleigh para distintos índices politrópicos, 
con 𝐿/𝐻 = 3. 71 
Figura 4.28. Número de conducción en función del número de Rayleighpara distintos factores de escala, con 
𝑛 = 1.5 72 
Figura 4.29. Recinto de modelado del fluido. 74 
Figura 4.30. Evolución de las componentes de la energía y de flujos de calor para los parámetros de la Tabla 
4.4. 78 
Figura 4.31. Evolución de las componentes de la energía y de flujos de calor para los parámetros de la Tabla 
4.13. 79 
Figura 4.32. Evolución de las componentes de la energía y de flujos de calor para los parámetros de la Tabla 
4.19. 80 
Figura 4.33. Balance energético correspondiente a la introducción de parámetro de la Tabla 4.19 81 
 
 
xix 
 
 
 
 
 
 
1 INTRODUCCIÓN 
 
 
El fenómeno de la convección natural está implicado en numerosos procesos de la vida que desempeñamos en 
nuestro planeta. Desde hechos cotidianos como el de cocinar un plato de pasta, para el que es necesario llevar a 
ebullición una cantidad de agua depositada en un recipiente, o el de soplar el café para enfriarlo, hasta ser parte 
del origen de fenómenos atmosféricos, como el viento, o geológicos, como el movimiento de las placas 
tectónicas en nuestro planeta, entre otros muchos. 
La convección es uno de los tres mecanismos básicos de transferencia de calor. La transferencia de calor, como 
dice Holman [1], es aquella ciencia que busca predecir la transferencia de energía que puede ocurrir entre cuerpos 
materiales, provocada por una diferencia de temperatura. En termodinámica, a dicha energía se le llama calor; 
en cambio, la temperatura es una propiedad intrínseca a la materia que se está midiendo, asociada a lo caliente 
(mayor temperatura) o fría (menor) que se encuentre. La convección, por tanto, consiste en la transferencia de 
calor desde la superficie de un sólido o de un fluido hacia un fluido (a diferente temperatura) que está en contacto 
a través de la interfase de separación, por la acción combinada de la difusión de calor y el transporte de masa. 
Sin contacto directo entre sólido/fluido y el fluido este proceso no sería posible, ni tampoco sin el transporte de 
masa en el seno del fluido hacia el que se transfiere el calor. Es decir, la convección es el mecanismo de 
transferencia de calor que se da por el transporte de un fluido. 
 Además de la convección, se tienen dos mecanismos más: la conducción y la radiación. La conducción, en 
cambio, es la transferencia de calor por acción de dos cuerpos que se encuentran a diferente temperatura, pero 
que no requiere de movimiento fluido. Por último, la radiación, que es el mecanismo que no requiere de un 
medio material para transferir calor, ya que se efectúa mediante ondas electromagnéticas, por lo que solo requiere 
de diferencia de temperaturas. 
Retomando la convección, existen dos formas de clasificar la convección según el origen del movimiento del 
fluido. Si el movimiento del fluido es causado por las fuerzas de flotabilidad que resultan de las variaciones de 
densidad originadas por variaciones de temperatura en el fluido, estaremos hablando de convección natural o 
libre. Si el movimiento del fluido es causado por la acción de un elemento externo, como una bomba o ventilador, 
entonces hablaremos de convección forzada. A lo largo de esta introducción se habla de la primera, también 
denominada convección de Rayleigh-Bénard, donde se ejemplifican las bases físicas del modelo que se va a 
desarrollar. 
1.1. Convección de Rayleigh-Bénard 
La convección natural se empezó a estudiar de manera sistemática a comienzos del siglo XX, siendo Henri 
Bénard quien, mediante una serie de experimentos, encontró un patrón de comportamiento en una capa de fluido, 
que era calentada desde abajo y por encima estaba expuesto al ambiente. En concreto, el fluido en 
experimentación era esperma de ballena, que goza de gran conductividad térmica, en capas delgadas. Dicho 
patrón consistía en, tras un régimen transitorio inicial, se alcanzaba un estado estable de convección (Figura 1.1), 
y la superficie quedaba dividida en un diagrama de celdas hexagonales, como si de un panal de abeja se tratara 
(Figura 1.2). 
 
 
 
 Introducción 
 
 
2 
 
En 1916, Lord Rayleigh [2] desarrolló una teoría en la que se daba por primera vez sobre papel una explicación 
a los mecanismos que originaban la convección (no así a las celdas que se formaban en el experimento de 
Bénard). Rayleigh postuló que para que se produzca la convección natural, además de tener en cuenta la 
flotación, es también imprescindible considerar los efectos de la fricción, originada por las fuerzas de viscosidad, 
y el de difusividad térmica (Figura 1.3), que tiende a homogeneizar el campo de temperaturas en el fluido y, por 
ende, el de densidades. 
 
 
Una manera de entenderlo es la siguiente. Una porción de fluido, mediante una perturbación, es desplazada hacia 
arriba, teniendo más temperatura que las partículas fluidas circundantes, y por tanto menor densidad. Esta 
diferencia de densidad empuja dicha porción (menos densa) a ascender, por flotabilidad, y a hundir la zona del 
fluido más denso por acción de la gravedad; a la misma vez, la masa que queda inmediatamente por encima 
dificulta su ascenso, por fricción viscosa; y la difusión térmica tiende a disminuir el gradiente térmico que 
produce la flotación. Si el mecanismo de flotación es más fuerte que los otros dos, entonces el fluido tenderá a 
ascender (equilibrio inestable); si es más fuerte la difusividad térmica y la fricción viscosa, entonces la porción 
de fluido permanecerá en reposo (equilibrio estable). 
En su análisis, Rayleigh consideró una capa de fluido entre dos placas infinitas paralelas a distinta temperatura, 
siendo la inferior la de mayor valor. 
La importancia relativa entre estos mecanismos se mide mediante el número de Rayleigh (Ra), que halló Lord 
Rayleigh al linealizar las ecuaciones de Navier-Stokes que gobiernan el dominio fluido. Relaciona las fuerzas 
de flotación con el producto de las fuerzas viscosas y la conductividad térmica. 
 
Figura 1.3 Esquema de los fenómenos involucrados en la convección de 
Rayleigh-Bénard. [13] 
Figura 1.2. Celdas de Bénard en esperma. 
Reproducción de una de las fotografías originales 
de Bénard 
Figura 1.1. Diagrama de celdas hexagonales de 
Bénard [13] 
 
 
3 
Ra =
𝑔 𝛽 Δ𝑇 𝐻3
𝜈 𝛼
 (1.1) 
 
Donde 𝑔 es la aceleración de la gravedad, 𝛽 es el coeficiente de dilatación térmica, Δ𝑇 es la diferencia de 
temperatura entre la superficie inferior y superior del fluido, 𝐻 es el espesor de la capa fluida, 𝜈 es la viscosidad 
cinemática y 𝛼 es la difusividad térmica del fluido. Este número adimensional es de gran relevancia a la hora de 
modelar la convección natural en un dominio fluido. 
En resumen, Bénard fue quien primero estudió experimentalmente el fenómeno de la convección natural, y 
Rayleigh quien dio el fundamento teórico. Por ello, se conoce como convección de Rayleigh- Bénard al proceso 
termoconvectivo de inestabilidad que puede producirse en un medio fluido, sujeto a una perturbación en el 
gradiente de temperaturas, dando lugar en su superficie superior a las ya mencionadas celdas de Bénard. 
1.2. Ejemplos de convección natural 
Este mecanismo es fácil de entender mediante ejemplos, que irán desde lo más cotidiano, hasta algunos bastante 
complejos, pero que son parte fundamental de la vida en nuestro planeta. 
El primer ejemplo, en el que se ven implicados los otros dos mecanismos, es el mencionado lineas atrás. Cuando 
se quiere llevar a ebullición una cantidad de agua depositada en una olla, se transfiere calor por radiación desde 
la vitrocerámica hacia la base de esta. A medida que la base va adquiriendo temperatura, el agua más próxima a 
ella adquiere calor, viendo alterada su densidad, por lo que asciende debido a la flotabilidad, desplazando la 
masa de agua que estaba por encima hacia abajo. Una vez esta masa entra en contacto con la base, se calienta 
mientras que la que antes ascenció, se enfría, repitiéndose el proceso anterior.Es decir, vuelve a subir por 
flotabilidad la masa de agua que está en contacto con la base, desplazando la masa que se encontraba por encima 
hacia abajo. 
 
Este proceso cíclico de desplazamiento de agua desde la base hacia la parte superior de la olla, y posterior bajada 
a la base, es un ejemplo cotidiano de transferencia de calor por convección natural. Además de la convección, 
también intervienen los otros dos mecanismos de transferencia antes mencionados: conducción, en el mango de 
la olla a través del calentamiento del metal; y radiación, por acción de la vitrocerámica 
Otro ejemplo de convección natural, a gran escala, es el que se produce en la circulación de la atmósfera terrestre. 
George Hadley, buscando una explicación a la dirección occidental que siguen los vientos alisios en las 
proximidades del Ecuador, expuso en 1735 [3] un modelo de circulación atmosférica en el que el aire sigue una 
circulación cerrada. Tal y como se observa en la ¡Error! No se encuentra el origen de la referencia., el aire 
húmedo y caliente (4) convergen en el Ecuador, debido a que el efecto Coriolis allí es despreciable comparado 
con latitudes mayores. Tras calentarse, por diferencia de densidad con el aire circundante, se eleva, formando 
tormentas eléctricas (1). Cuando alcanza la tropopausa, no pueden subir más por flotabilidad, ni tampoco 
permanecer allí, ya que se ve empujado por el aire que viene desde abajo, por lo que se dirigen, según el 
Figura 1.4. Esquema de los mecanismos de transferencia de calor 
en el calentamiento de agua de una olla 
 
 Introducción 
 
 
4 
hemisferio, bien hacia el norte (2a), o bien hacia el sur (2b). A medida que el aire va alejándose del ecuador, va 
enfriándose por intercambio de calor con el aire circundante, por lo que comienza a descender (3). También 
aumenta en este desplazamiento del aire el efecto Coriolis, lo que va desviando el aire superior hacia el Este 
(visto desde un observador terrestre). Cuando el aire llega a la superficie terrestre, el anticiclón lo empuja hacia 
el Ecuador (4), viendo un observador terrestre cómo también se dirige hacia el Oeste por el efecto Coriolis, 
dando explicación al movimiento de los vientos alisios hacia el Oeste. 
 
 
 
Hadley creyó que únicamente había una celda de circulación por hemisferio, es decir, que el aire viajaba desde 
el Ecuador hacia los polos directamente. Pero en el siglo XIX se demostró que el efecto Coriolis también 
provocaba que las celdas cerradas de circulación se cerraran antes de llegar a los polos, dando origen a dos más 
por hemisferio (Figura 1.6). Se mantuvo el nombre de célula de Hadley para la que iba desde el Ecuador a 
latitudes aproximadas a los 30º, célula de Ferrel o de latitudes medias para la que se encuentra entre la célula de 
Hadley y la Polar (hasta latitudes cercanas a los 60º), y célula Polar para la que llega hasta cada polo. Este modelo 
de tres celdas convectivas es algo simplificado, ya que la extensión de estas es variable y depende de otros 
factores, como la estación del año o el calentamiento global, pero es útil para prever el movimiento de los vientos 
a gran escala. 
Otros ejemplos donde podemos encontrar la convección en la naturaleza son en el manto superior de la Tierra, 
donde se transporta calor por el movimiento de grandes cantidades de material de gran viscosidad, desde el 
interior (a 3773 K) de esta hacia su superficie (a 873 K), y que da origen al movimiento de las placas tectónicas 
(Figura 1.7); o en el Sol, donde su zona convectiva alcanza una temperatura de 2 millones de Kelvin en su radio 
interno y de unos 5700 K en su radio externo (Figura 1.8). Los datos térmicos han sido extraídos de [4] y Figura 
1.8. 
Como se detallará más adelante, para que un sistema que contenga un fluido sea estable a convección, debe 
cumplirse el criterio de Schwarzschild, un criterio matemático que permite predecir tal mecanismo bajo unas 
determinadas circunstancias. 
1.3. Objetivo del proyecto. Estructura. 
Este proyecto es una reformulación del proyecto realizado en [5], donde se modeló el comportamiento de un 
dominio fluido según el problema de convección de Rayleigh-Bénard en dos recintos: uno rectangular y otro 
anular. En dichos problemas no se tenía en cuenta la radiación, y se empleó la hipótesis anelástica, es decir, se 
Figura 1.6. Modelo de circulación atmosférica de tres 
células convectivas [23] 
Figura 1.5. Modelo de circulación atmosférica 
propuesto por George Hadley [24] 
 
 
5 
calculaban también las variaciones de densidad en la coordenada vertical/radial para unas condiciones de 
contorno en las que se fijaba la temperatura en la base y en la superficie superior del recinto. 
El objetivo de este proyecto es hacer una extensión de las ecuaciones obtenidas en [5] para un recinto rectangular, 
realizando consideraciones termo-fluidomecánicas y aspectos energéticos, y establecer comparaciones. 
Consiste, por tanto, en el estudio de problemas de convección natural con hipótesis anelástica en recintos 
rectangulares, sin tener en cuenta la radiación, para temperaturas impuestas en la base y la superficie superior 
del dominio fluido. 
La estructura del proyecto comienza con una formulación del problema, en la que se habla del modelo de 
Boussinesq, que fue la que Rayleigh tomó de partida para su estudio teórico; y de la hipótesis anelástica, que es 
una generalización del modelo de Boussinesq considerando la estratificación de la densidad. El siguiente 
capítulo versa sobre el método de cálculo que se ha empleado para la resolución de las ecuaciones de Navier-
Stokes, presentando el método de diferencias finitas. En el cuarto capítulo se presentan los resultados obtenidos 
para distintas combinaciones de parámetros, haciendo un profundo análisis fluido-termomecánico, y realizando 
un balance energético del proceso. Por último, se habla de las conclusiones del proyecto, y de las posibles líneas 
de desarrollo a partir del mismo. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 1.8. Convección en el manto terrestre [21] 
Figura 1.7. Capas del Sol [22] 
 
 
 
2 FORMULACIÓN DEL PROBLEMA 
 
 
Este capítulo comienza presentando físicamente la naturaleza del fenómeno de la convección térmica, como 
resultado de una inestabilidad en el campo termofluido. Matemátimente, se verá que dicha inestabilidad viene 
provocada por un gradiente de temperaturas menor al del estado adiabático. Tras ello, se presentan dos modelos 
para estudiar la convección térmica de Rayleigh- Bénard: la aproximación de Boussinesq, y la hipótesis 
anelástica. En el desarrollo de ambos modelos, especialmente en el que se emplea la hipótesis anelástica, se han 
seguido las formulaciones y obtención de las ecuaciones de Juárez [5], Glatzmaier [6], y Randall [7], que son 
ampliamente conocidas y empleadas en el estudio de fenómenos convectivos. Entre otros aspectos de interés, se 
presentará un estado de referencia basado en polítropos. La obtención final de las ecuaciones del modelo con 
hipótesis anelástica conduce a la aplicación, en el capítulo posterior, de un método numérico de resolución no 
tan conocido en la resolución de estos problemas. 
2.1. Estabilidad a convección de un fluido 
La convección natural sucede como consecuencia de una inestabilidad en el campo de temperaturas de una 
atmósfera. Para entender cuál es el origen del fenómeno, se va a dar una explicación física del mismo, y también 
un criterio matemático, que ayuda a predecirlo. 
Para la explicación física, siguiendo un enfoque similar al realizado en [6] y [8], es preciso considerar una 
porción de fluido dentro de una atmósfera (indicados con el subíndice p y Ω, respectivamente), inicialmente en 
equilibrio con su entorno en temperatura, presión y densidad, al se llamará posición (1). 
 
𝑇𝑝1 = 𝑇Ω1 , 𝑝𝑝1 = 𝑝Ω1 , 𝜌𝑝1 = 𝜌Ω1 
 
A continuación, se desplaza dicha porción hacia arriba,a una nueva posición (2), lo suficientemente rápido como 
para que no haya transferencia de calor, pero lo suficientemente lento como para que se mantenga el equilibrio 
de presiones. Es decir, el desplazamiento se da de tal manera que el proceso sea reversible y adiabático, y por 
tanto isentrópico, donde la velocidad del fluido sea muy inferior a la del sonido. 
Una vez llegados a la posición (2), existen dos escenarios posibles: 
 
• Situación subadiabática: la temperatura de la atmósfera circundante es mayor a la de la porción 
fluida, y con ello, la densidad del entorno es inferior a la de la porción, por lo que ésta, si 
desaparece la perturbación inicial, desciende por flotabilidad a su posición original (¡Error! No se 
encuentra el origen de la referencia.). 
 
𝑇𝑝2 < 𝑇Ω2 , 𝜌𝑝2 > 𝜌Ω2 
 
A este proceso se le llama onda de gravedad interna, y en este caso, el fluido es estable a convección. Se debe 
a que la atmósfera circundante se encuentra en estado subadiabático con respecto a la porción fluida. 
 
 
 Formulación del problema 
 
 
8 
• Situación superadiabática: la temperatura de la atmósfera circundante es inferior a la de la porción 
fluida, siendo la densidad del entorno superior a la de la porción, por lo que ésta sigue ascendiendo 
por flotabilidad (Figura 2.2). 
 
𝑇𝑝2 > 𝑇Ω2 , 𝜌𝑝2 < 𝜌Ω2 
 
A este proceso se le llama convección térmica, y en este caso, el fluido es inestable a convección. Se debe a 
que la atmósfera circundante se encuentra en estado superadiabático con respecto a la porción fluida. De esta 
manera, dependiendo del estado en que se encuentre la atmósfera circundante con respecto a la porción 
perturbada, podremos anticipar en si el recinto fluido será estable o no a convección. 
 
 
Otro enfoque, que va a servir como puente para conectar con el criterio matemático para predecir la convección, 
es hacer un balance de energías, tal y como se propone en [9]. 
 Se supone que una partícula, en las mismas condiciones que el enfoque anterior en (1) ≡ (𝑇, 𝑧), se hace 
ascender a una posición (2) ≡ (𝑇 − 𝑑𝑇, 𝑧 + 𝑑𝑧), de tal manera que en el proceso no pierde energía interna. Al 
llegar a (2), la partícula aún conserva su temperatura inicial 𝑇, que difiere con la temperatura de su entorno en 
una cantidad 𝑑𝑇. Si se asume que, mientras la partícula libera el exceso de energía interna que tiene con su 
entorno, la presión local se mantiene constante, la cantidad de energía térmica (por unidad de masa) liberada es 
igual a 𝑐𝑝𝑑𝑇, siendo 𝑐𝑝 el calor específico a presión constante. 
Otra contribución energética al proceso, que se opone a la liberación de energía térmica, es la energía potencial 
gravitatoria de la partícula, definida (por unidad de masa) como 𝑔𝑑𝑧, siendo 𝑔 la aceleración de la gravedad. 
Por tanto, el trabajo total empleado en desplazar la partícula desde (1) hasta (2) es: 
 
Partícula adiabática 
Atmósfera subadiabática 
 
Figura 2.1. Estabilidad a convección térmica: estado subadiabático. Evolución de la 
temperatura y la densidad frente a la altura. 
 
 
9 
ΔW|(1)→(2) = 𝑔𝑑𝑧 + 𝑐𝑝𝑑𝑇 (2.1) 
 
 
Si la energía potencial gravitatoria es, en valor absoluto, superior a la energía interna liberada por la partícula, 
ésta volverá a su posición inicial; esta situación se corresponde con el caso subadiabático antes descrito. Pero si 
la energía interna supera a la energía potencial gravitatoria, la partícula no solo se mantendrá en (2), si no que 
subirá hasta que la energía liberada por la partícula se compense con la gravitatoria; es decir, existirá una 
inestabilidad, dando lugar a la convección térmica. 
Por tanto, el valor límite del gradiente de temperaturas que establece la frontera con la inestabilidad (estado 
adiabático, ad) es aquel que anula la expresión (2.1): 
 
ΔW|(1)→(2) = 0 →
𝑑𝑇
𝑑𝑧
|
𝑎𝑑
= −
𝑔
𝑐𝑃
 (2.2) 
 
Este valor del gradiente de temperaturas es el límite del criterio de Schwarzschild. Este criterio establece que 
para valores del gradiente de temperaturas menores (es decir, más negativos, ya que la temperatura decrece con 
el aumento de la coordenada z) que −𝑔/𝑐𝑃, el medio en estudio es inestable, originándose la convección térmica 
(se tiene un estado superdiabático de la atmósfera circundante). En cambio, para valores mayores del gradiente 
térmico (es decir, más cercanos a 0), el medio es estable a convección (se tiene un estado subadiabático de la 
atmósfera circundante). Otro punto de vista de la obtención este criterio, con idéntica conclusión, se puede 
obtener en [10]. 
Partícula adiabática 
Atmósfera superadiabática 
 
Figura 2.2. Inestabilidad a convección térmica: estado superadiabático. Evolución de 
la temperatura y la densidad frente a la altura. 
 
 Formulación del problema 
 
 
10 
2.2. Modelos de convección 
El interés de este proyecto reside en modelar la zona convectiva de una estrella, como el Sol, o atmósferas o 
zonas interiores de ciertos planetas que tienen un comportamiento convectivamente inestable. En ambos casos, 
dicha zona es, suponiendo que el astro tenga un volumen perfectamente esférico, corona circular. Sin embargo, 
si se descompone tal corona en muchos sectores consecutivos, preservando siempre que la anchura de esos 
sectores sea suficientemente mayor a la altura de la corona, se consigue buenos resultados aproximando tales 
sectores por prismas rectangulares, dada las grandes dimensiones que los astros en cuestión. En las paredes 
laterales se aplican condiciones de contorno que permiten aislar un sector, para contribuir a una mayor 
simplificación del modelo. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Por tanto, se va a modelar un volumen fluido en un recinto prismático rectangular (Figura 2.3), con dimensiones 
de anchura 𝐿, profundidad, 𝐵, y altura, 𝐻 (asociadas a las coordenadas 𝑥, 𝑦, 𝑧, respectivamente). La base del 
recinto está a una temperatura fija 𝑇1, y la superficie superior o tapadera se encuentra a una temperatura también 
fija 𝑇2, donde 𝑇1 > 𝑇2. Las paredes laterales son adiabáticas, es decir, no transfieren calor al exterior. Además, 
se tiene en cuenta la acción de la gravedad sobre el fluido. 
 Para resolver este problema se han empleado las ecuaciones generales de la Mecánica de Fluidos. La resolución 
completa de las mismas para los procesos convectivos es inviable en la práctica, ya que son altamente no lineales 
y con términos diferenciales. Además, incluyen el efecto de las ondas acústicas en las perturbaciones de presión 
y densidad. Tales ondas viajan normalmente a una velocidad mucho mayor que la del fluido en los planetas y 
los interiores de estrella, por lo que requieren de un paso de computación muy pequeño para poder ser resueltos. 
Los modelos de Boussinesq y anelástico filtran este efecto, y sus resultados serán, tanto más válidos, en la 
medida en que se verifique que la velocidad del fluido en estudio sea mucho menor que la del sonido. Por ello, 
han sido los modelos escogidos para profundizar en este proyecto. 
La forma genérica del sistema de ecuaciones en forma diferencial de Navier Stokes del que se parte para el 
desarrollo del proyecto está formada por las ecuaciones (2.3 - (2.5). 
 
𝜕𝜌
𝜕𝑡
= − ∇ ∙ (𝜌�⃗�) (2.3) 
𝐷�⃗�
𝐷𝑡
= −
1
 𝜌
∇𝑝 +
1
𝜌
∇ ∙ �̿� + �⃗� (2.4) 
Figura 2.3. Recinto prismático rectangular de modelado del volumen fluido. 
 
 
11 
𝜌
𝐷𝑒
𝐷𝑡
+ 𝑝 ∇ ∙ �⃗� = ∇ ∙ (𝑘∇𝑇) + 𝑄 (2.5) 
 
La ecuación (2.3) es la de conservación de la masa, que verifica que la variación de densidad 𝜌 por unidad de 
tiempo 𝑡 es igual a la negativa de la divergencia del flujo de masa 𝜌�⃗�, siendo �⃗� la velocidad del fluido. La 
ecuación (2.4) es la Segunda Ley de Newton aplicada a un dominio fluido (también llamada ecuación de 
conservación de cantidad de movimiento). En ella intervienen, además de otras variables ya mencionadas, la 
presión (𝑝), el tensor de viscosidad�̿�, y la gravedad �⃗�. La tercera ecuación del sistema (2.5) es la Primera Ley 
de la Termodinámica para un gas perfecto, donde intervienen entre otras la energía interna específica 𝑒, la 
conductividad térmica 𝑘 (𝑘 = 𝑐𝑝𝜌𝜅, donde 𝑐𝑝 es la capacidad calorífica a presión constante y 𝜅 la difusividad 
térmica), la temperatura 𝑇, y término energético asociado a la viscosidad u otras aportaciones en forma de calor, 
como de tipo nuclear o eléctrico 𝑄. 
 
2.2.1. Aproximación de Boussinesq 
En el artículo [11] es descrita la aproximación atribuida a Boussinesq (1903) basándola en dos hipótesis: a) la 
perturbación en densidad originadas con el movimiento del fluido (gas perfecto) resulta de los efectos de la 
temperatura, y no de la presión; b) en las ecuaciones de continuidad y cantidad de movimiento, la variación de 
densidad es despreciable, salvo cuando estén acopladas a la aceleración de la gravedad en la fuerza de flotación. 
Esto mismo es aplicable siempre y cuando se cumpla que: 
• La dimensión de la extensión vertical donde se encuentra el fluido (gas perfecto) es mucho menor 
que cualquier escala, ya sea de presión, temperatura, o densidad. 
• El movimiento induce perturbaciones en densidad y presión que no debe exceder, en orden de 
magnitud, el valor de referencia estático de estas variables. 
• La velocidad del fluido en estudio sea muy inferior a la velocidad del sonido; es decir, que el Mach 
del fluido sea inferior a 0.1. 
Además, bajo estas condiciones, las ecuaciones son formalmente equivalentes a las de un sistema incompresible 
cuando el gradiente de temperatura se sustituye por su exceso adiabático, y se puede intercambiar el calor 
específico a presión constante, 𝑐𝑝, por el correspondiente a volumen constante, 𝑐𝑉 . 
Una escala de altura en densidad es, por ejemplo, equivalente a −(
1
𝜌
𝑑𝑙𝑛𝜌
𝑑𝑧
)
−1
, siendo z la altura del dominio. 
Esto suele ser así en los océanos, mantos o fluidos de núcleos de planetas terrestres, incluso en el interior de las 
estrellas, donde la escala de altura de densidad es mucho más grande que el radio del cuerpo. Sin embargo, no 
da buenos resultados en las atmósferas de planetas gigantes ni en la zona exterior de las estrellas, ya que la 
estratificación de la densidad es mucho mayor que en los casos anteriores. 
Estas consideraciones de la aproximación de Boussinesq, llevadas al sistema de ecuaciones (2.3) - (2.5), se 
traducen en lo siguiente. Con un desarrollo similar a [6], se tiene que las variaciones en densidad, al asumirse 
muy pequeñas frente a su valor de referencia, la ecuación de conservación de la masa (2.3) queda simplificada: 
 
∇ ∙ �⃗� = 0 (2.6) 
 
También se consideran pequeños los efectos locales de la presión sobre la densidad, dependiendo las 
perturbaciones en densidad únicamente de la temperatura, (𝜌 = −𝜌0𝛼𝑇). Esto, unido a la ecuación de equilibrio 
hidrostático (𝑑𝑝0/𝑑𝑧 = −�̅�𝜌0) y a considerar la viscosidad del fluido, la gravedad y la densidad constantes, da 
lugar a: 
 
 
 Formulación del problema 
 
 
12 
𝜕�⃗�
𝜕𝑡
= −(�⃗� ∙ ∇)�⃗� −
1
𝜌0
∇𝑝 + �̅�𝛽𝑇�⃗⃗� + 𝜈∇2�⃗� (2.7) 
 
Donde 𝜌0 es la densidad media del dominio fluido, 𝑔 ̅es la gravedad, 𝛽 es el coeficiente de dilatación y 𝜈 es la 
difusividad viscosa. La simplificación en la ecuación de conservación de la masa es la que hace que el modelo 
de Boussinesq filtre las ondas acústicas, ya que la divergencia de la ecuación (2.7) elimina las derivadas con 
respecto al tiempo debido a la ecuación (2.6). De esta manera, resulta que la presión en cualquier parte del 
dominio y en cada instante de tiempo está dada por perturbaciones en la temperatura y la velocidad. 
Para un gas perfecto, la densidad de energía interna es igual a 𝑐𝑣𝑇, siendo 𝑐𝑣 la capacidad calorífica a volumen 
constante. Por tanto, asumiendo también que los efectos viscosos son despreciables, la ecuación (2.5) para la 
aproximación de Boussinesq quedaría de la siguiente manera: 
 
𝜕𝑇
𝜕𝑡
= −(�⃗� ∙ ∇)𝑇 + 𝛼∇2𝑇 (2.8) 
 
Estas simplificaciones que permiten pasar del sistema de ecuaciones (2.3) – (2.5) al (2.6) – (2.8) reduce mucho 
el costo computacional requerido para simular la convección térmica frente a otros modelos que permiten una 
mayor estratificación de la densidad; la rápida obtención de resultados es la gran ventaja que aporta este modelo. 
Las ecuaciones (2.6) y (2.8) son idénticas a las empleadas de partida en los proyectos [12] y [13] para los casos 
sin radiación. La ecuación (2.7) difiere únicamente en que aquí se ha asumido que 𝛽 ≈ 1/𝑇0. 
A continuación, se presenta el modelo con hipótesis anelástica, el cual es una generalización de la aproximación 
de Boussinesq que permite la estratificación de la densidad. En la práctica, no se han desarrollado dos códigos 
en paralelo que resuelvan las ecuaciones (2.6) – (2.8) por un lado, y las del modelo anelástico por otro; como se 
verá, existe una manera de particularizar el modelo anelástico y obtener la aproximación de Boussinesq, sin 
necesidad de discretizar ambos sistemas de ecuaciones. 
 
2.2.2. Modelo anelástico 
La densidad en los mantos y núcleos de los planetas del sistema solar no suele variar más de un 20% desde el 
límite inferior al superior de los recintos, siendo válidos (en primera instancia) los resultados que proporciona la 
aproximación de Boussinesq. Sin embargo, esta aproximación no es válida para modelos globales de simulación 
de convección en planetas gigantes (diez veces más másicos que la Tierra) o en la zona externa de las estrellas, 
ya que la estratificación de la densidad en estos casos no es despreciable. 
En estos casos, se requiere un modelo que considere las variaciones de la densidad, midiendo como el fluido se 
expande o se contrae, y cómo esto afecta a la transferencia de calor, cantidad de movimiento y la energía cinética. 
En 1953, Batchelor desarrolla unas ecuaciones de movimiento en derivadas parciales, asumiendo que el orden 
de magnitud de la temperatura potencial era pequeño, y que el tiempo de escala estaba dado por la frecuencia de 
Brunt.Väisälä, omitiendo las ondas acústicas. En 1960, Charney y Ogura emplean para el análisis de la 
convección en atmósferas un sistema de ecuaciones análogo al de Batchelor, para un gas perfecto, asumiendo 
que las distribuciones de densidad y presión son cercanas a las de una atmósfera estratificada adiabáticamente. 
A estas ecuaciones, Ogura y Phillips en 1962 [14] las denominaron ecuaciones anelásticas. 
Para el desarrollo del modelo anelástico hay varias formas de llegar al sistema de ecuaciones objetivo, pero todas 
se basan en las dos mismas hipótesis: la velocidad del fluido es muy inferior a la del sonido y las perturbaciones 
termodinámicas son pequeñas comparadas con sus valores medios. 
Se ha seguido una estructura de explicación del modelo anelástico similar a la realizada en [6]. En primer lugar, 
se va a presentar la ecuación de estado y también las distintas hipótesis empleadas. Posteriormente, se detalla la 
obtención de las ecuaciones de conservación de la masa, cantidad de movimiento y energía interna del modelo. 
 
 
13 
Ecuación de estado 
Las variables termodinámicas se descomponen en la suma de un estado de referencia hidrostático y la 
perturbación, de ahora en adelante con la siguiente notación: 
 
{
�̅�(�⃗�, 𝑡) = 𝑇𝐻(𝑧) + 𝑇(�⃗�, 𝑡)
�̅�(�⃗�, 𝑡) = 𝜌𝐻(𝑧) + 𝜌(�⃗�, 𝑡)
�̅�(�⃗�, 𝑡) = 𝑝𝐻(𝑧) + 𝑝(�⃗�, 𝑡)
 (2.9) 
 
Donde �̅� denota la variable termodinámica total y el subíndice 𝑋𝐻 denota el estado de referencia, por lo que la 
variable 𝑋 denota la perturbación. Tal y como se ha dicho anteriormente, empleando la notación de (2.9), debe 
cumplirse que 𝑋 ≪ 𝑋𝐻. 
 
 El estado de referencia de las variables termodinámicas es solo función de la coordenada vertical, y están en 
equilibrio hidrostático, es decir, cumplen que 
 
𝑑𝑝𝐻
𝑑𝑧
= −�̅�𝜌𝐻 (2.10) 
 
El equilibriohidrostático en una estrella es un estado de estabilidad que se alcanza dentro la misma cuando se 
compensan la presión térmica (que tiende a la expansión) y su peso gravitatorio (que tiende a la compresión). 
Las variables hidrostáticas se escogen de tal manera que tengan una estratificación isentrópica, o cuasi-
isentrópica al menos, dentro de la zona de convección (𝑑𝑆𝐻/𝑑𝑧 ≈ 0). Tal y como se indica en [6], 
indistintamente se hablará de proceso isentrópico y adiabático cuando se haga referencia a procesos de 
estratificación de la temperatura o de la densidad. 
Para modelar el fluido, se utiliza la ecuación de gas ideal, 𝑝 + 𝑝𝐻 = 𝑅𝑔(𝜌 + 𝜌𝐻)(𝑇 + 𝑇𝐻), donde 𝑅𝑔 es la 
relación entre el gas y la masa atómica del mismo. Dicha ecuación es buena aproximación para interiores de 
estrellas. Empleando las variables de estado en su forma hidrostática, se tiene también que 
 
𝑝𝐻 = 𝑅𝑔𝜌𝐻𝑇𝐻 (2.11) 
 
Para obtener el comportamiento del fluido considerando las magnitudes totales no hay más que sumar a cada 
magnitud hidrostática su correspondiente perturbación. En todo momento se asumirá que 𝑅𝑔 , 𝑐𝑣 , 𝑐𝑝 son 
constantes para un mismo fluido, por lo que las relaciones existentes entre ellas son 𝑅𝑔 = 𝑐𝑃 − 𝑐𝑣 , 𝛾 = 𝑐𝑃/𝑐𝑣. 
La linealización de la ecuación de estado para las perturbaciones es: 
 
𝑝
𝑝𝐻
=
𝜌
𝜌𝐻
+
𝑇
𝑇𝐻
 (2.12) 
 
La entropía para un gas perfecto en equilibrio hidrostático queda, considerando una constante arbitraria 𝑆0, como 
 
(𝑆 + 𝑆𝐻) = 𝑐𝑣 ln(𝑝 + 𝑝𝐻) − 𝑐𝑃 ln(𝜌 + 𝜌𝐻) + 𝑆0 
 
 Formulación del problema 
 
 
14 
 
En términos únicamente de variables hidrostáticas, se tiene 
 
𝑆𝐻 = 𝑐𝑣 ln 𝑝𝐻 − 𝑐𝑃 ln 𝜌𝐻 + 𝑆0 (2.13) 
 
Derivando esta última ecuación con respecto a la altura, junto con la ecuación de equilibrio hidrostático (2.10) 
y de la ley de gases ideales (2.11), se obtiene el estado de referencia del gradiente de la entropía, 
 
𝑑𝑆𝐻
𝑑𝑧
=
𝑐𝑃
𝑇𝐻
(
𝑑𝑇𝐻
𝑑𝑧
+
�̅�
𝑐𝑃
) (2.14) 
 
Para que el proceso sea isentrópico (o adiabático), la ecuación (2.14) se tiene que anular. Es decir, el gradiente 
de temperatura hidrostática ha de ser tal que 𝑑𝑆𝐻/𝑑𝑧 sea cero. Por lo tanto, 
 
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
= −
�̅�
𝑐𝑃
 (2.15) 
 
Este resultado vuelve a mostrar el criterio de Schwarzschild, ya obtenido en la expresión (2.2) para un caso 
genérico. En primera aproximación, para modelos de estrellas donde se cumpla el equilibrio hidrostático, se 
cumple que: 
𝑑𝑇𝐻
𝑑𝑧
→
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
 (2.16) 
 
 
Volviendo a la ecuación (2.12), empleando (2.13) se puede reescribir de la siguiente manera, 
 
𝑑𝑆𝐻
𝑑𝑧
=
𝑐𝑃
𝑇𝐻
[ 
𝑑𝑇𝐻
𝑑𝑧
−
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
] (2.14.1) 
 
Las escalas de altura, ya mencionadas con anterioridad, siguen la siguiente formulación: 
 
𝐻𝑇 = (
1
𝑇𝐻
𝑑𝑇𝐻
𝑑𝑧
)
−1
= 1/ℎ𝑇 
 𝐻𝜌 = (
1
𝜌𝐻
𝑑 𝜌𝐻
𝑑𝑧
)
−1
= 1/ℎ𝜌 
𝐻𝑝 = (
1
𝑝𝐻
𝑑𝑝𝐻
𝑑𝑧
)
−1
= 1/ℎ𝑝 
(2.17) 
 
 
15 
𝐻𝜃 = (
1
𝜃𝐻
𝑑𝜃𝐻
𝑑𝑧
)
−1
= 1/ℎ𝜃 
 
Dichas escalas son de importancia, ya que sus valores, en comparación con las dimensiones del recinto que 
alberga el volumen fluido, permiten simplificar o no términos en las ecuaciones generales. Su uso es conveniente 
para hacer un análisis de órdenes de magnitud con los términos de las ecuaciones del modelo. 
Para desarrollos posteriores, es útil saber que el corchete de la ecuación (2.14.1), empleando también (2.10) y 
(2.11), y la definición de inversa de la escala de densidad introducida en (2.17), puede escribirse como sigue, 
 
[ 
𝑑𝑇𝐻
𝑑𝑧
−
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
] =
1
𝛾
 [
𝑑𝑇𝐻
𝑑𝑧
− ℎ𝜌(𝛾 − 1)𝑇𝐻 ] (2.18) 
 
Para la obtención de las ecuaciones de conservación de la masa y de la cantidad de movimiento de nuestro 
modelo anelástico, se ha seguido un análisis de órdenes de magnitud, similar al seguido por Randall [7] y 
también empleado en [5]. Para la obtención de la ecuación de la energía, se han empleado procedimientos tanto 
de [7] como de [6], además de emplear el concepto de temperatura potencial (𝜃), introducido en [13]: 
 
𝜃 = 𝑇 (
𝑝
𝑝0
)
𝑅
𝑐𝑝
 (2.19) 
 
Se define como la temperatura que adquiría un volumen fluido para llevarlo de una presión p a una presión 
media de referencia, 𝑝0. Se conserva en los procesos adiabáticos y sirve de medida de la estabilidad a convección 
Por último, antes de la deducción de las ecuaciones de nuestro modelo anelástico, cabe señalar que están 
particularizadas para un recinto rectangular bidimensional, de altura 𝐻 y anchura 𝐿, por lo que las magnitudes y 
derivadas en la tercera dimensión (según eje 𝑦 en este caso) se omiten. 
Ecuación de conservación de la masa 
Tomando como partida la ecuación (2.3), se reescribe teniendo en cuenta la descomposición de (2.9) como sigue: 
 
𝜕(𝜌 + 𝜌𝐻)
𝜕𝑡
+ ∇ ∙ (𝜌�⃗� + 𝜌𝐻�⃗�) = 0 (2.20) 
 
Las magnitudes hidrostáticas no varían con respecto al tiempo, por lo que se anula su derivada temporal. 
Desarrollando la divergencia del flujo, y teniendo en cuenta la hipótesis de que la densidad hidrostática es mayor 
que la perturbada (𝜌/𝜌𝐻 ≪ 1), se llega a lo siguiente: 
 
𝜕𝜌
𝜕𝑡
+ 𝜌𝐻 (1 +
𝜌
𝜌𝐻
)∇ ∙ �⃗� + 𝑣𝑥 
𝜕𝜌
𝜕𝑥
+ 𝑣𝑧 (
𝜕𝜌
𝜕𝑧
+
𝜕𝜌𝐻
𝜕𝑧
) ≅ 
≅
𝜕𝜌
𝜕𝑡
+ 𝜌𝐻 (
𝜕𝑣𝑥
𝜕𝑥
+
𝜕𝑣𝑧
𝜕𝑧
 ) + 𝑣𝑥 
𝜕𝜌
𝜕𝑥
+ 𝑣𝑧 (
𝜕𝜌
𝜕𝑧
+
𝜕𝜌𝐻
𝜕𝑧
) = 0 
 
(2.21) 
 
 Formulación del problema 
 
 
16 
Se definen los valores característicos del tiempo, 𝜏; longitud vertical, 𝐻; longitud horizontal, 𝐿; velocidad 
horizontal, 𝑉; y velocidad vertical, 𝑉𝑧. A continuación, se realiza un análisis comparativo del término 
𝜌𝐻 𝜕𝑣𝑧/𝜕𝑧 con el resto de la ecuación (2.21). 
 
• |
𝜕𝜌
𝜕𝑡
| 
|
𝜕𝜌
𝜕𝑡
|
|𝜌𝐻
𝜕𝑣𝑧
𝜕𝑧
|
~
|𝜌|
𝜏
𝜌𝐻
|𝑉𝑧|
𝐻
=
|𝜌|
𝜌𝐻
𝐻
|𝑉𝑧|𝜏
=
𝑔
|𝜌|
𝜌𝐻
|𝑉𝑧|
𝜏
𝐻
𝑔 𝜏2
 (2.22) 
 
Si 𝑔|𝜌|/𝜌𝐻 es del orden de |𝑉𝑧/𝜏|, la derivada temporal de la densidad es despreciable cuando se 
verifique que 
𝜏2 ≫
𝐻
𝑔
 (2.23) 
 
• |𝑣𝑧
𝜕𝜌𝐻
𝜕𝑧
| 
|𝑣𝑧
𝜕𝜌𝐻
𝜕𝑧 |
|𝜌𝐻
𝜕𝑣𝑧
𝜕𝑧 |
=
1
𝐻𝜌
𝑣𝑧
𝜕𝑣𝑧
𝜕𝑧
~
1
𝐻𝜌
𝑉𝑧
𝑉𝑧
𝐻
=
𝐻
𝐻𝜌
 
 
 
Si 𝐻𝜌 ≪ 𝐻, se podría despreciar el término |𝜌𝐻𝜕𝑣𝑧/𝜕𝑧|,. Sin embargo, esto sería buena aproximación 
en el caso de movimientos superficiales. En cambio, el interés del modelo anelástico es mejorar los 
resultados de la aproximación de Boussinesq, y al mantener este término se tienen en cuenta los efectos 
de la estratificación de la densidad. 
 
• 𝑣𝑥
𝜕𝜌
𝜕𝑥
 
𝑣𝑥
𝜕𝜌
𝜕𝑥
|𝜌𝐻
𝜕𝑣𝑧
𝜕𝑧 |
~
|𝜌|𝑉
𝐿
𝜌𝐻
|𝑉𝑧|
𝐻
=
|𝜌|
𝜌𝐻
𝑉
|𝑉𝑧|
𝐻
𝐿
 (2.24) 
 
 
Si se cumple que 
𝑉
|𝑉𝑧|
𝐻
𝐿
≤ 1 (2.25) 
El término |𝑣𝑥
𝜕𝜌
𝜕𝑥
| puede ser despreciado, siempre y cuando la relación de aspecto 𝐿/𝐻 sea lo 
suficientemente grande. 
 
 
 
17 
• |𝑣𝑧
𝜕𝜌
𝜕𝑧
| 
|𝑣𝑧
𝜕𝜌
𝜕𝑧
|
|𝜌𝐻
𝜕𝑣𝑧
𝜕𝑧 |
~
|𝜌 𝑉𝑧|
𝐻
|𝜌𝐻
𝑉𝑧
𝐻
|
=
|𝜌|
𝜌𝐻
 (2.26) 
 
Verificándose la hipótesis inicial, en que la contribución hidrostática es mucho mayor a la contribución 
perturbada, se puede despreciar también este término. 
 
Finalmente, se llega a la ecuación de conservación de la masa para nuestro modelo anelástico: 
 
∇ ∙ (𝜌𝐻�⃗�) = 0 (2.27) 
 
 
Otra forma de escribir la ecuación (2.27) es la siguiente, que será de utilidad más adelante: 
 
∇ ∙ �⃗� = −ℎ𝜌𝑣𝑧 (2.28) 
 
Conservación de la cantidad de movimiento 
Se parte de la ecuación (2.4), aplicando nuevamente la descomposición realizada en (2.9) y el desarrollo de 
derivada sustancial. Se supone que la viscosidad dinámica es constante en el espacio (𝜌𝐻𝜈 = 𝑐𝑡𝑒), 
 
𝜕�⃗�
𝜕𝑡
+ �⃗� ∙ ∇�⃗� = −
1
𝜌 + 𝜌𝐻
∇(𝑝 + 𝑝𝐻) − �̅��⃗⃗� + 𝜈 (∇
2�⃗� +
1
3
∇(∇ ∙ �⃗�)) (2.29) 
 
Se van a analizar los dos primeros términos del segundo miembro, descomponiéndolos en componente 
horizontal y componente vertical. La componente horizontal se obtiene demanera inmediata 
 
−
1
𝜌 + 𝜌𝐻
𝜕(𝑝 + 𝑝𝐻)
𝜕𝑥
= −
1
𝜌𝐻 (1 +
𝜌
𝜌𝐻
)
𝜕𝑝
𝜕𝑥
≅ −
1
𝜌𝐻
𝜕𝑝
𝜕𝑥
= −
𝜕
𝜕𝑥
(
𝑝
𝜌𝐻
) (2.30) 
 
La componente vertical, aplicando la ecuación de equilibrio hidrostático (2.10) 
 
 
 Formulación del problema 
 
 
18 
−
1
𝜌 + 𝜌𝐻
(
𝜕𝑝
𝜕𝑧
+
𝜕𝑝𝐻
𝜕𝑧
) − �̅� = −
1
𝜌 + 𝜌𝐻
𝜕𝑝
𝜕𝑧
+ (
𝜌𝐻
𝜌 + 𝜌𝐻
− 1) �̅� = −
1
𝜌 + 𝜌𝐻
𝜕𝑝
𝜕𝑧
−
𝜌
𝜌 + 𝜌𝐻
�̅�
= −
1
𝜌𝐻 (
𝜌
𝜌𝐻
− 1)
( 
𝜕𝑝
𝜕𝑧
− 𝜌�̅�) ≅ −
1
𝜌𝐻
 
𝜕𝑝
𝜕𝑧
−
𝜌
𝜌𝐻
�̅� =
= −
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) −
𝑝
𝜌𝐻
(
1
𝜌𝐻
𝜕𝜌𝐻
𝜕𝑧
) −
𝜌
𝜌𝐻
�̅� 
 
(2.31) 
Ahora, aplicando la definición de temperatura potencial (2.18) sobre la ley de los gases ideales, se llega a la 
siguiente expresión de la densidad hidrostática: 
𝜌𝐻 =
𝑝0
𝑅𝑔𝜃𝐻
(
𝑝𝐻
𝑝0
)
1
 𝛾
 
Derivando la misma con respecto a la coordenada z, se llega a: 
 
1
𝜌𝐻
𝜕𝜌𝐻
𝜕𝑧
=
1
𝛾𝑝𝐻
𝜕𝑝𝐻
𝜕𝑧
−
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
 → 
1
𝐻𝜌
=
1
𝛾
1
𝐻𝑝
−
1
𝐻𝜃
 (2.32) 
 
Considerando nuevamente la hipótesis de pequeñas perturbaciones frente a estado hidrostático en todas las 
variables (𝜌 ≪ 𝜌𝐻 , 𝑝 ≪ 𝑝𝐻 , 𝑇 ≪ 𝑇𝐻 , 𝜃 ≪ 𝜃𝐻), se tiene: 
 
𝜌
𝜌𝐻
=
1
𝛾
𝑝
𝑝𝐻
−
𝜃
𝜃𝐻
 (2.33) 
 
Aplicando las relaciones (2.32) y (2.33) a (2.31), así como la ecuación de equilibrio hidrostático (2.10): 
 
−
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) −
𝑝
𝜌𝐻
(
1
𝜌𝐻
𝜕𝜌𝐻
𝜕𝑧
) −
𝜌
𝜌𝐻
�̅�
= −
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) −
𝑝
𝜌𝐻
(
1
𝛾𝑝𝐻
𝜕𝑝𝐻
𝜕𝑧
−
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
) − (
1
𝛾
𝑝
𝑝𝐻
−
𝜃
𝜃𝐻
)𝑔 ̅ = 
= −
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) −
𝑝
𝜌𝐻
(
1
𝛾𝑝𝐻
(−𝜌𝐻�̅�) −
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
) − (
1
𝛾
𝑝
𝑝𝐻
−
𝜃
𝜃𝐻
) �̅� =
= −
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) +
𝑝
𝜌𝐻
(
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
) + (
𝜃
𝜃𝐻
) �̅� 
(2.34) 
Para este caso, se vuelve a hacer un análisis de órdenes de magnitud para simplificar la ecuación con respecto al 
término | 
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) |. Antes, es preciso definir la frecuencia de Brunt-Väislälä, 𝑁, que aparece de forma natural en 
la descripción de ondas gravitacionales: 
 
 
19 
 
𝑁2 =
�̅�
𝜃𝐻
𝑑𝜃𝐻
𝑑𝑧
=
�̅�
𝐻𝜃
 (2.35) 
 
• |
𝑝
𝜌𝐻
(
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
) | 
|
𝑝
𝜌𝐻
(
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
) |
|
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) |
~
𝑝
𝜌𝐻
𝑁2
�̅�
 
𝑝
𝜌𝐻
1
𝐻
=
𝑁2𝐻
𝑔
=
𝐻
𝐻𝜃
 (2.36) 
 
Si se cumple que 𝐻 ≪ 𝐻𝜃 , podemos despreciar el término 
𝑝
𝜌𝐻
(
1
𝜃𝐻
𝜕𝜃𝐻
𝜕𝑧
); es decir, siempre que la escala 
de temperatura potencial sea mucho mayor a la altura del recinto del volumen fluido. 
 
• |
𝜕
𝜕𝑥
(
𝑝
𝜌𝐻
) | 
|
𝜕
𝜕𝑥
(
𝑝
𝜌𝐻
) |
|
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) |
~
𝑝
𝜌𝐻
1
𝐿
 
𝑝
𝜌𝐻
1
𝐻
=
𝐻
𝐿
 
 
(2.37) 
Si la anchura del recinto en estudio es lo suficientemente mayor con respecto a la altura, podemos 
despreciar el término |
𝜕
𝜕𝑥
(
𝑝
𝜌𝐻
) |, obtenido en (2.30). 
 
Se quiere comprobar ahora si 𝜃/𝜃𝐻 ≅ 𝑇/𝑇𝐻. Para ello, recordando la expresión (2.32) y asumiendo que las 
escalas de temperatura y temperatura potencial son suficientemente mayores a las de densidad y presión, se 
deduce que 𝛾𝐻𝑝 ≅ 𝐻𝜌 , por lo que (2.33) se puede reescribir, despejando el cociente de temperaturas 
potenciales, de la siguiente manera: 
 
𝜃
𝜃𝐻
=
1
𝛾
𝑝
𝑝𝐻
𝑑𝑝𝐻
𝑑𝑧
𝑑𝑝𝐻
𝑑𝑧
 −
𝜌
𝜌𝐻
= −
𝑝
𝛾𝐻𝑝𝜌𝐻�̅�
−
𝜌
𝜌𝐻
≅ −
𝑝
𝐻𝜌𝜌𝐻�̅�
−
𝜌
𝜌𝐻
 (2.38) 
 
Por otro lado, tomando logaritmos en la ecuación de gas perfecto, y derivando con respecto a la coordenada z, 
se llega fácilmente a 
 
1
𝜌𝐻
𝑑𝜌𝐻
𝑑𝑧
=
1
𝑝𝐻
𝑑𝑝𝐻
𝑑𝑧
−
1
𝑅𝑔𝑇𝐻
𝑑𝑇𝐻
𝑑𝑧
→
1
𝐻𝜌
=
1
𝐻𝑝
−
1
𝑅𝑔𝐻𝑇
 (2.39) 
 
 
 Formulación del problema 
 
 
20 
Donde volviendo a aplicar la hipótesis 𝐻𝜌, 𝐻𝑝 ≪ 𝐻𝑇 , se obtiene que 𝐻𝑝 ≅ 𝐻𝜌. Con esta última deducción, 
retomando la expresión (2.34) y volviendo a considerar la ecuación de equilibrio hidrostático en el último paso, 
 
𝜃
𝜃𝐻
≅ −
𝑝
𝐻𝜌𝜌𝐻�̅�
−
𝜌
𝜌𝐻
≅ −
𝑝
𝐻𝑝𝜌𝐻�̅�
−
𝜌
𝜌𝐻
=
𝑝
𝑝𝐻
𝐻𝑝
𝐻𝑝
−
𝜌
𝜌𝐻
=
𝑝
𝑝𝐻
−
𝜌
𝜌𝐻
 
(2.40) 
 
Comparando con la ecuación de estado para perturbaciones linealizada (2.12), se tiene finalmente que 𝜃/𝜃𝐻 ≅
𝑇/𝑇𝐻. Es importante hacer énfasis de que este resultado es válido siempre y cuando las perturbaciones de las 
escalas de altura de densidad y presión sean suficientemente pequeñas comparadas con las escalas de altura de 
temperatura y temperatura potencial. 
También se puede comprobar, derivando logarítmicamente la expresión de temperatura potencial (2.18), que las 
escalas de altura de temperatura y temperatura potencial son similares en primera aproximación, bajo la hipótesis 
mencionada. 
 
Empleando las expresiones obtenidas en (2.30) y (2.34) y la aproximación 𝜃/𝜃𝐻 ≅ 𝑇/𝑇𝐻, la ecuación de 
conservación de cantidad de movimiento (2.29) puede reescribirse como sigue, 
 
𝜕�⃗�
𝜕𝑡
+ �⃗� ∙ ∇�⃗� = −∇(
𝑝
𝜌𝐻
) +
𝑝
𝜌𝐻
(
1
𝑇𝐻
𝜕𝑇𝐻
𝜕𝑧
) +
𝑇
𝑇𝐻
�̅� + 𝜈 (∇2�⃗� +
1
3
∇(∇ ∙ �⃗�)) (2.41) 
 
Ecuación de conservación de la energía interna 
Tomando como partida la ecuación (2.5), se desarrollan sus términos recordando nuevamente la descomposición 
(2.9), y considerando gas perfecto para la energía interna (𝑒 = 𝑐𝑣𝑇): 
 
𝑐𝑣 𝜌𝐻
𝜕(𝑇 + 𝑇𝐻)
𝜕𝑡
+ 𝑐𝑣(𝜌 + 𝜌𝐻)�⃗� ∙ ∇(𝑇 + 𝑇𝐻) + (𝑝 + 𝑝𝐻) ∇ ∙ �⃗� = ∇ ∙ (𝑘∇(𝑇 + 𝑇𝐻)) + �̅�𝑣 (2.42) 
 
Se ha asumido en el término de la derivada temporal que 𝜌 ≪ 𝜌𝐻. El término �̅�𝑣 representa la contribución al 
calor de los efectos viscosos. En [14] se describe como la ratio de disipación de energía mecánica originada por 
la viscosidad. Su expresión viene dada por: 
 
�̅�𝑣 = 2𝜇 (𝑒𝑖𝑗𝑒𝑖𝑗 −
1
3
(∇ ∙ �⃗�)2), 
 𝑒𝑖𝑗 =
1
2
(
𝜕𝑣𝑖
𝜕𝑥𝑗
+
𝜕𝑣𝑗
𝜕𝑥𝑖
) , 𝑖, 𝑗 = 𝑥, 𝑧 
(2.43) 
 
Se va a estudiar cada miembro de la ecuación (2.42) por separado. Empezando por el de la izquierda, 
desarrollando sus términos y teniendo presente la ecuación de conservación de la masa en su forma (2.28): 
 
 
 
21 
𝑐𝑣 𝜌𝐻
𝜕𝑇
𝜕𝑡
+ 𝑐𝑣(𝜌𝐻 + 𝜌)(�⃗� ∙ ∇𝑇 + 𝑣𝑧
𝑑𝑇𝐻
𝑑𝑧
) − (𝑝 + 𝑝𝐻)ℎ𝜌𝑣𝑧 ≅
≅ 𝑐𝑣 𝜌𝐻
𝜕𝑇
𝜕𝑡
+ 𝑐𝑣𝜌𝐻�⃗� ∙ ∇𝑇 + 𝑐𝑣𝜌𝐻𝑣𝑧 (
𝑑𝑇𝐻
𝑑𝑧
−
ℎ𝜌𝑝𝐻
𝑐𝑣𝜌𝐻
)+ (𝑐𝑣𝜌𝑣𝑧
𝑑𝑇𝐻
𝑑𝑧
− ℎ𝜌𝑝𝑣𝑧) 
(2.44) 
 
En la segunda igualdad de la ecuación (2.44) se han eliminado los términos de orden superior. Para un gas 
perfecto en equilibrio hidrostático, se puede emplear la siguiente relación: 
 
ℎ𝜌𝑝𝐻
𝑐𝑣𝜌𝐻
= (𝛾 − 1)ℎ𝜌𝑇𝐻 (2.45) 
 
Además, reescribiendo la ecuación de estado linealizada utilizando la ley de gases ideales (2.11), tenemos 
 
𝑝
𝑅𝑔𝜌𝐻𝑇𝐻
=
𝜌
𝜌𝐻
+
𝑇
𝑇𝐻
→
𝑝
𝜌𝐻
=
𝜌𝑅𝑔𝑇𝐻
𝜌𝐻
+ 𝑅𝑔𝑇 (2.46) 
 
Por tanto, empleando las relaciones (2.45) y (2.46) en los dos últimos términos de (2.44) se tiene, 
respectivamente: 
 
𝑐𝑣𝜌𝐻𝑣𝑧 (
𝑑𝑇𝐻
𝑑𝑧
−
ℎ𝜌𝑝𝐻
𝑐𝑣𝜌𝐻
) = 𝑐𝑝𝜌𝐻𝑣𝑧 ( 
𝑑𝑇𝐻
𝑑𝑧
−
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
) (2.47) 
𝑐𝑣𝜌𝐻𝑣𝑧 (
𝜌
𝜌𝐻
𝑑𝑇𝐻
𝑑𝑧
−
ℎ𝜌𝑝
𝑐𝑣𝜌𝐻
) = −𝑐𝑣𝜌𝐻𝑣𝑧((𝛾 − 1)ℎ𝜌𝑇− 𝛾(
𝑑𝑇𝐻
𝑑𝑧
−
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
)
𝜌
𝜌
𝐻
 ) (2.48) 
 
En (2.47) y (2.48) se ha empleado la equivalencia (2.18). 
 
El término proporcional a 𝜌/𝜌𝐻 en (2.48) es de orden inferior al resto, por lo que será despreciado en lo que 
sigue. Teniendo en cuenta todo lo anterior, aplicando las expresiones (2.47) y (2.48) sobre (2.42), la ecuación de 
conservación de la energía interna queda de la siguiente manera: 
 
 
𝜕𝑇
𝜕𝑡
= −�⃗� ∙ ∇𝑇+ (𝛾 − 1)ℎ𝜌𝑇𝑣𝑧− 𝛾(
𝑑𝑇𝐻
𝑑𝑧
− (
𝑑𝑇𝐻
𝑑𝑧
)
𝑎𝑑
)𝑣𝑧 +
1
𝜌
𝐻
𝑐𝑣
 (∇2𝑇+
𝑑2𝑇𝐻
𝑑𝑧2
++�̅�𝑣) (2.49) 
 
 Observando los términos primero y segundo del segundo miembro, y usando (2.28), se puede reescribir la 
ecuación (2.49): 
 
 
 Formulación del problema 
 
 
22 
𝜕𝑇
𝜕𝑡
= −∇ ∙ (𝑇�⃗�) + (𝛾 − 2)ℎ𝜌𝑇𝑣𝑧
− 𝛾(
𝑑𝑇𝐻
𝑑𝑧
− (
𝑑𝑇𝐻
𝑑𝑧
)
𝑎𝑑
)𝑣𝑧 +
1
𝜌
𝐻
𝑐𝑣
 (∇2𝑇+
𝑑
2
𝑇𝐻
𝑑𝑧2
+ �̅�𝑣) (2.50) 
 
En resumen, para obtener el modeloanelástico se han considerado las siguientes hipótesis: 
 
• Las variables termodinámicas se desvían muy poco del estado de equilibrio. 
• La velocidad del fluido es mucho menor que la del sonido. 
• Las escalas de temperatura de presión y densidad son mucho menores que las de temperatura y 
temperatura potencial. 
 
El sistema de ecuaciones que finalmente da forma a nuestro modelo anelástico es el conformado por el sistema 
(2.51) - (2.53) 
 
∇ ∙ (𝜌𝐻�⃗�) = 0 
(2.51) 
 
𝜕�⃗�
𝜕𝑡
+ �⃗� ∙ ∇�⃗� = −
𝜕
𝜕𝑧
(
𝑝
𝜌𝐻
) +
𝑝
𝜌𝐻
(
1
𝑇𝐻
𝜕𝑇𝐻
𝜕𝑧
) + (
𝑇
𝑇𝐻
) �̅� + 𝜈 (∇2�⃗� +
1
3
∇(∇ ∙ �⃗�)) 
(2.52) 
 
𝜕𝑇
𝜕𝑡
= −∇ ∙ (𝑇�⃗�) + (𝛾 − 2)ℎ𝜌𝑇𝑣𝑧− 𝛾(
𝑑𝑇𝐻
𝑑𝑧
− (
𝑑𝑇𝐻
𝑑𝑧
)
𝑎𝑑
)𝑣𝑧 +
1
𝜌
𝐻
𝑐𝑣
 (∇2𝑇+
𝑑
2
𝑇𝐻
𝑑𝑧2
+ �̅�𝑣) (2.53) 
 
2.3. Resolución del modelo anelástico 
El sistema de ecuaciones anelástico (2.51) – (2.53) lo conforman 𝑖 + 3 ecuaciones, siendo n la dimensión 
espacial. Al tratarse de un problema bidimensional, 𝑖 = 2. El número de incógnitas es también 𝑖 + 3, es decir, 
las componentes de la velocidad (𝑣𝑥 , 𝑣𝑧) y las perturbaciones de las tres magnitudes de estados (𝜌, 𝑝, 𝑇). Las 
magnitudes de referencia serán conocidas (𝜌𝐻 , 𝑝𝐻 , 𝑇𝐻), una vez se imponga un estado de referencia que cumpla 
el equilibrio hidrostático. La aceleración de la gravedad se asume constante y conocida (�̅�), así como las 
propiedades del fluido (𝛾, 𝑐𝑣 , 𝑐𝑃). 
Para reducir el costo computacional en la resolución del modelo, interesa reducir al mínimo el número de 
ecuaciones, obteniendo el campo de temperaturas y velocidades en todo el espacio para cada instante de tiempo. 
También se requiere imponer unas condiciones de contorno, y realizar la adimensionalización del sistema 
completo, para reducir el número de parámetros. Una vez se haya llegado a ese punto, en el capítulo posterior 
se explica el método numérico, y se discretiza el sistema de ecuaciones que permite implementar 
computacionalmente el modelo para obtener los resultados. 
 
2.3.1. Ecuación de la vorticidad. Función de corriente 
Una forma de reducir el número de ecuaciones en este modelo es empleando una función cuyo comportamiento 
cumpla con la ecuación de conservación de la masa, y que pueda suplir a alguna de las incógnitas de las 
ecuaciones restantes. Se define para ello la función de corriente, Ψ, cuya relación con las componentes (𝑥, 𝑧) de 
la velocidad tiene la siguiente formulación: 
 
 
23 
 
{
 
 𝑣𝑥 = −
1
𝜌𝐻
𝜕Ψ
𝜕𝑧
𝑣𝑧 =
1
𝜌𝐻
𝜕Ψ
𝜕𝑥
 (2.54) 
 
De manera inmediata, se puede comprobar como dicha función cumple la ecuación de conservación de la masa 
para el modelo anelástico (2.51). Se define la vorticidad, 𝜔, que sigue la siguiente expresión: 
 
𝜔 = ∇ × �⃗� = −
1
𝜌𝐻
(∇2Ψ− ℎ𝜌
𝜕Ψ
𝜕𝑧
) (2.55) 
 
La vorticidad no se obtendrá de manera explícita en el problema, pero sí será útil para simplificar la escritura de 
las ecuaciones a lo largo del proyecto. Cuando se vayan a obtener las ecuaciones definitivas a implementar en 
el código numérico, éstas se expresarán en términos de la función de corriente. Esto se debe a que, la 
implementación de las condiciones de contorno es más intuitiva y sencilla de hacer sobre la función de corriente 
que sobre la vorticidad. Para simplificar la ecuación de conservación de la cantidad de movimiento conviene 
recordar algunas relaciones matemáticas, relacionadas con el rotacional y la divergencia de una función. Así, 
considerando las magnitudes 𝑓, �⃗�, vectoriales ambas: 
 
(𝑓 ∙ ∇)𝑓 = ∇(
1
2
 |𝑓|
2
) − 𝑓 × (∇ × 𝑓) (2.56) 
∇ × (𝑓 × �⃗�) = (∇ ∙ �⃗�)𝑓 + (�⃗� ∙ ∇)𝑓 − (∇ ∙ 𝑓)�⃗� − (𝑓 ∙ ∇)�⃗� (2.57) 
 
Además, se emplea que el rotacional de un gradiente es nulo, al igual que la divergencia de un rotacional. 
Tomando rotacional en ambos miembros de la ecuación de conservación de cantidad de movimiento (2.52), y 
empleando la definición de vorticidad, junto con las expresiones (2.56) - (2.57): 
 
𝜕𝜔
𝜕𝑡
= −�⃗� ∙ ∇𝜔 + ℎ𝜌𝜔𝑣𝑧 − (
�̅�
𝑇𝐻
𝜕𝑇
𝜕𝑥
+
1
𝜌𝐻𝑇𝐻
𝜕𝑝
𝜕𝑥
𝑑𝑇𝐻
𝑑𝑧
) +
𝜇
𝜌𝐻
∇2𝜔+ 𝜇
𝑑
𝑑𝑧
(
1
𝜌𝐻
)(
𝜕𝜔
𝜕𝑧
−
4
3
ℎ𝜌
𝜌𝐻
𝜕2Ψ
𝜕𝑥2
) (2.58) 
 
Tras esta transformación, las incógnitas que quedan en la ecuación son 𝑝, 𝑇,Ψ. El término que incluye la 
derivada horizontal de la perturbación en presión se puede calcular con la componente horizontal de la ecuación 
de cantidad de movimiento: 
 
1
𝜌𝐻
𝜕𝑝
𝜕𝑥
= −
𝜕𝑣𝑥
𝜕𝑡
− ((�⃗� ∙ ∇)�⃗�)
𝑥
+
𝜇
𝜌𝐻
(∇2𝑣𝑥 −
1
3
ℎ𝜌
𝜕𝑣𝑧
𝜕𝑥
) (2.59) 
 
Sin embargo, tal y como se comenta en [6], el término que incluye el gradiente horizontal de la presión, dentro 
de la aproximación anelástica, empeora la conservación de la energía. Además, no calcular este término implica 
 
 Formulación del problema 
 
 
24 
simplificar el código, ya que no se requiere calcular la componente horizontal de la ecuación, por lo que de ahora 
en adelante se suprime de la ecuación (2.58). Esta simplificación será tanto más buena en función de si se cumple, 
recordando el análisis (2.36), que 𝐻 ≪ 𝐻𝜃 . De esta manera, se elimina la presión como incógnita dentro de la 
ecuación. 
Teniendo en consideración todo lo anterior, la ecuación de conservación de cantidad de movimiento quedaría 
reescrita de la siguiente manera: 
 
𝜕𝜔
𝜕𝑡
= −�⃗� ∙ ∇𝜔 + ℎ𝜌𝜔𝑣𝑧 −
�̅�
𝑇𝐻
𝜕𝑇
𝜕𝑥
+
𝜇
𝜌𝐻
∇2𝜔 + 𝜇
𝑑
𝑑𝑧
(
1
𝜌𝐻
)(
𝜕𝜔
𝜕𝑧
−
4
3
ℎ𝜌
𝜌𝐻
𝜕2Ψ
𝜕𝑥2
) (2.60) 
 
Los dos primeros términos del segundo miembro pueden agruparse en uno solo, como se muestra en [6]: 
 
−�⃗� ∙ ∇𝜔 + ℎ𝜌𝜔𝑣𝑧 = −[𝑣𝑥
𝜕𝜔
𝜕𝑥
+ 𝑣𝑧 (
𝜕𝜔
𝜕𝑧
− ℎ𝜌𝜔)] = −∇ ∙ (𝜔�⃗�) (2.61) 
 
Este cálculo es más preciso numéricamente que calculando cada término por separado. 
Finalmente, el modelo anelástico propuesto para calcular el campo de temperaturas y el campo de velocidades 
(por medio del campo de función de corriente) en cada instante y en cada punto del espacio, está compuesto por 
dos siguientes ecuaciones: 
 
𝜕𝜔
𝜕𝑡
= −∇ ∙ (𝜔�⃗�) −
�̅�
𝑇𝐻
𝜕𝑇
𝜕𝑥
+
𝜇
𝜌𝐻
∇2𝜔 + 𝜇
𝑑
𝑑𝑧
(
1
𝜌𝐻
)(
𝜕𝜔
𝜕𝑧
−
4
3
ℎ𝜌
𝜌𝐻
𝜕2Ψ
𝜕𝑥2
) (2.62) 
𝜕𝑇
𝜕𝑡
= −∇ ∙ (𝑇�⃗�) + (𝛾 − 2)ℎ𝜌𝑇𝐻𝑣𝑧
− 𝛾(
𝑑𝑇𝐻
𝑑𝑧
− (
𝑑𝑇𝐻
𝑑𝑧
)
𝑎𝑑
)𝑣𝑧 +
1
𝜌
𝐻
𝑐𝑣
 (∇2𝑇+
𝑑2𝑇𝐻
𝑑𝑧2
+ �̅�𝑣) (2.63) 
 
En lo referente a este sistema de ecuaciones cabe indicar que, aunque aparezcan términos que incluyan a la 
velocidad (o a alguna de sus componentes) y/o a la vorticidad, todas ellas dependen directamente de la función 
de corriente, que es sobre la que se aplicarán las condiciones de contorno. Es decir, las incógnitas en este sistema 
son el campo de temperaturas y el campo de funciones de corriente. 
 
2.3.2. Estado de referencia: Polítropos 
Para la resolución del modelo, es necesario establecer una relación adicional entre las variables de estado 𝑝, 𝜌, 𝑇 
que nos permita despejar una de ellas en función de las demás, y así hacer el sistema de ecuaciones del modelo 
anelástico compatible determinado. Una elección apropiada sería adoptar un modelo de polítropos, tal y como 
se propone en [5] y se desarrolla en [6] donde la presión y la densidad satisfacen la condición de equilibrio 
hidrostático (2.10), y la presión es proporcional a la densidad, elevada a potencia de (𝑛 + 1)/𝑛 : 
 
𝑝𝐻 = 𝑝1 (
𝜌𝐻
𝜌1
)
𝑛+1
𝑛
 (2.64) 
 
 
 
25 
Donde 𝑛 es el índice politrópico. Las magnitudes con subíndice 1 se corresponden con valores del estado de 
referencia para la presión y la densidad en la superficie inferior del recinto. Este modelo de estado de referencia 
es propuesto por ajustarse bastante bien al comportamiento de temperatura, presión y densidad en el interior de 
las estrellas y planetas. Por ejemplo, según [6], para 𝑛 = 0 (sin estratificación de la densidad), se tiene una buena 
aproximación en el comportamiento del manto y núcleos de planetas terrestres; para polítropos con 𝑛 =1.5, 
considerando gas perfecto, se tiene una estratificación de la temperatura adiabática, lo cual es apropiado para 
una zona de convección; para polítropos con 𝑛 = 3 se aproxima el comportamiento interior de las estrellas que 
están dominadas por la presión de radiación, como el núcleo del Sol. 
Considerándose el caso en estudio, de un recinto rectangular con altura 𝐻 y longitud 𝐿, en que la temperatura en 
la superficie inferior es mayor que la superficie superior (𝑇1 > 𝑇2), y la aceleración de la gravedad es constante 
y con valor −�̅��⃗⃗�, la relación politrópica se escribe de la siguiente manera: 
 
𝑝𝐻(𝑧) = 𝑝1Θ
𝑛+1(𝑧), 𝜌𝐻(𝑧) = 𝜌1Θ
𝑛(𝑧) (2.65) 
 
Donde Θ es el polítropo adimensional. Su expresión se deduce sustituyendo (2.65) en la ecuación de equilibrio 
hidrostático (2.10) con la condición Θ(0) = 1: 
 
Θ(𝑧) = 1 −
�̅�𝜌1
(𝑛 + 1) 𝑝1
𝑧 = 1 −
𝑧
𝑧1
 (2.66) 
 
Donde sea ha definido 𝑧1 = (𝑛 + 1)𝑝1/(𝜌1�̅�). La presión y temperaturas del estado de referencia, usando la 
expresión (2.66), quedan de la siguiente manera: 
 
𝑝𝐻(𝑧) = 𝑝1 (1 −
�̅�𝜌1
(𝑛 + 1) 𝑝1
𝑧)
𝑛+1
, 𝜌𝐻(𝑧) = 𝜌1 (1 −
�̅�𝜌1
(𝑛 + 1) 𝑝1
𝑧)
𝑛
 (2.67) 
 
Para hallar la expresión de la temperatura en función del polítropo adimensional, se emplea la ley de los gases 
ideales (2.11), 
 
𝑇𝐻(𝑧) =
𝑝𝐻(𝑧)
𝑅𝑔 𝜌𝐻(𝑧)
= 𝑇1𝛩(𝑧) = 𝑇1 (1 −
�̅�
(𝑛 + 1) 𝑅𝑔𝑇1
𝑧) = 𝑇1 −
�̅�
(𝑛 + 1) 𝑅𝑔
𝑧 = 𝑇1 −
Δ𝑇
𝐻
𝑧 (2.68) 
 
A la igualdad final se llega imponiendo la condición de contorno en la superficie superior del recinto. Para que 
se cumpla que 𝑇𝐻(𝑧 = 𝐻) = 𝑇2, Δ𝑇/𝐻 = (𝑇1 − 𝑇2)/𝐻 =
�̅� 
(𝑛+1) 𝑅𝑔
 . 
Por la expresión (2.68) se pueden sacar dos conclusiones inmediatas: a) la evolución de la temperatura 
hidrostática es lineal decreciente (Figura 2.4), y b) es independiente del índice politrópico. La expresión del 
estado de referencia del gradiente de la temperatura hidrostática quedaría como 
 
𝑑𝑇𝐻
𝑑𝑧
= −
𝑇1
𝑧1
= −
�̅�
(𝑛 + 1)𝑅𝑔
 (2.69) 
 
 Formulación del problema 
 
 
26 
 
 
Comparando este último resultado con el criterio de Schwarzschild (2.15), es inmediato obtener el índice 
politrópico que nos lleva a tal estado: 
 
𝑑𝑇𝐻
𝑑𝑧
|
𝑎𝑑
 = −
�̅�
𝑐𝑃
= −
�̅�
(𝑛 + 1)𝑅𝑔
→ 𝑛 =
𝑐𝑃
𝑅𝑔
− 1 =
𝛾
𝛾 − 1
− 1 → 𝑛𝑎𝑑 =
1
𝛾 − 1
 (2.70) 
 
 Este valor del índice politrópico es de suma importancia ya que marca la frontera de la estabilidad a convección. 
De este modo, el proceso será subadiabático (estable a convección) si 𝑛 > 1/(𝛾 − 1); en cambio, será 
superadiabático (convección térmica) si 𝑛 < 1/(𝛾 − 1). Dicho de otra manera, el criterio de Schwarzschild 
para nuestro polítropo predice que para un índice politrópico superior a 1/(𝛾 − 1), el medio estelar será estable 
a convección. 
Otro aspecto relevante es que, para zonas de estrellas donde se cumpla el equilibrio hidrostático, se caractericen 
por números de Rayleigh muy altos, y el índice politrópico sea cercano a 𝑛𝑎𝑑, la relación (2.16) es más próxima; 
es decir, la variación de la temperatura hidrostática con la altura tenderá a la variación de la temperatura 
hidrostático para el caso adiabático. En el presente proyecto, al ser lineal la expresión de la temperatura 
hidrostática (2.68), la relación (2.16) se cumplirá en todo momento. 
Por último, cabe hacer un apunte adicional respecto a la utilidad de los polítropos. Aunque se ha planteado el 
sistema de ecuaciones del modelo de Boussinesq, este no es más que un caso particular del modelo anelástico. 
En concreto, considerando el fluido sin estratificación de la densidad, es decir, 𝑛 = 0, el sistema de ecuaciones 
del modelo anelástico (2.51) – (2.53) equivaldría al de Boussinesq. Por tanto, como ya se indicó, para mayor 
agilidad y simplicidad en la obtención y comparación de los resultados de ambos métodos, se impondrá 𝑛 = 0 
para sacar resultados de la aproximación de Boussinesq, en vez de crear un nuevo código que resuelva el sistema 
(2.6) – (2.8). 
 
Figura 2.4. Evolución de las magnitudes de referencia adimensionalizadas frente a la altura, con Δ𝑇 =
100𝐾, 𝑇1 = 303 𝐾: temperatura hidrostática (izquierda) y densidad hidrostática (derecha). 
 
 
27 
2.3.3. Condiciones de contorno. Condiciones iniciales 
El recinto en que se va a modelar el fluido, como ha sido mencionado con anterioridad, es bidimensional. Consta 
de una cavidad rectangular, con altura 𝐻 y anchura 𝐿. Las temperaturas impuestas en toda la superficie inferior 
y superior son, respectivamente, 𝑇1 y 𝑇2, cumpliéndose siempre que 𝑇1 > 𝑇2. El sistema de referencia a emplear 
es cartesiano, tal que (𝑥, 𝑧) ∈ [0, 𝐿] × [0, 𝐻]. El intervalo de tiempo en que se evaluará el modelo será tal que 
𝑡 ∈ [0, τ], siendo 𝜏 el instante final de la simulación. La gravedad se asume constante en todo el recinto, y se 
expresa de manera que �⃗� = −�̅��⃗⃗�. Las paredes verticales se consideran adiabáticas y sin fricción. 
 
 
 
 
 
 
 
 
 
 
 
 
 
El campo de temperaturas y de velocidades estará condicionado por estas características impuestas en el recinto. 
En forma matemática, tales condiciones de contorno se expresan de la siguiente manera: 
 
 
{ 
(𝑇 + 𝑇𝐻)|𝑧=0 = 𝑇1
(𝑇 + 𝑇𝐻)|𝑧=𝐻 = 𝑇2
, ∀ 𝑡 ∈ (0, 𝜏) → {
 𝑇|𝑧=0 = 0
 𝑇|𝑧=𝐻 = 0
, ∀ 𝑡 ∈ (0, 𝜏) 
 
{ 
𝜕
𝜕𝑥
(𝑇 + 𝑇𝐻)|𝑥=0 = 0
𝜕
𝜕𝑥
(𝑇 + 𝑇𝐻)|𝑥=𝐿 = 0
, ∀ 𝑡 ∈ (0, 𝜏) → 
{
 
 
𝜕𝑇
𝜕𝑥
|
𝑥=0
 = 0
 
𝜕𝑇
𝜕𝑥
|
𝑥=𝐿
 = 0
, ∀ 𝑡 ∈ (0, 𝜏) 
(2.71) 
 
{ 
�⃗�|𝑧=0 = �⃗�|𝑧=𝐻 = 0
�⃗�|𝑥=0 = �⃗�|𝑥=𝐿 = 0
, → 
{
 
 Ψ|𝑧=0 = Ψ|𝑧=𝐻 =
𝜕Ψ
𝜕𝑧
|
𝑧=0
=
𝜕Ψ
𝜕𝑧
|
𝑧=𝐻
= 0
 Ψ|𝑥=0 = Ψ|𝑥=𝐿 =
𝜕Ψ
𝜕𝑥
|
𝑥=0
=
𝜕Ψ
𝜕𝑥
|
𝑥=𝐿
= 0
, ∀ 𝑡 ∈ (0, 𝜏) (2.72) 
 
Las condiciones iniciales corresponden a una perturbación inicial en la función de corriente: 
 
Figura 2.5. Recinto rectangular de modelado del fluido 
 
 Formulación del problema 
 
 
28 
{
𝑇|𝑡=0 = 𝑇0(𝑥, 𝑧)
Ψ|𝑡=0 = Ψ0(𝑥, 𝑧)
 (𝑥, 𝑧) ∈ [0, 𝐿] × [0, 𝐻] (2.73) 
 
Las propiedades del fluido (𝛾, 𝑐𝑃 , 𝑐𝑣) , como ya se mencionó con anterioridad, se asumen constantes. 
 
2.3.4. Adimensionalización 
Con el objetivo de reducir el número de parámetros dentro de las ecuaciones (2.62) - (2.63), se propone una 
adimensionalización similar a la adoptada en [5]. La variable adimensional sigue la notación 𝑋∗. 
 
𝑥 = 𝐻𝑥∗, 𝑧 = 𝐻𝑧∗, 𝑡 =
𝐻2
𝛼
𝑡∗ 
 �⃗� =
𝛼
𝐻
�⃗�∗, Ψ = 𝜌1𝛼Ψ
∗, ω =
𝛼
𝐻2
ω∗ 
 𝜌1 = 𝜌1𝜌
∗, 𝜌𝐻 = 𝜌1𝜌𝐻
∗ , 𝑝 = 𝑝1𝑝
∗ =
𝛼2𝜌1
𝐻2
𝑝∗, 𝑝𝐻 = 𝑝1𝑝𝐻
∗ =
𝛼2𝜌1
𝐻2
𝑝𝐻
∗ 
 𝑇 = Δ𝑇 𝑇∗, 𝑇𝐻 = Δ𝑇 𝑇𝐻
∗ 
(2.74) 
 
Donde 𝜌1 = 𝜌𝐻(𝑧 = 0), 𝛼 = 𝑘/𝑐𝑃𝜌1 es el coeficiente de difusividad térmica. Se introducen también los 
números adimensionales de Rayleigh (Ra) y Prandtl (Pr), 
 
Ra =
�̅�Δ𝑇 𝐻3
𝑇𝑚𝑒𝑑𝜈𝛼
, Pr =
𝜈
𝛼
 (2.75) 
 
 Donde 𝑇𝑚𝑒𝑑 = (𝑇1 + 𝑇2)/2, y 𝜈 = 𝜇/𝜌1, coeficiente de viscosidad cinemática. El número de Rayleigh, como 
ya se describió en el capítulo 1, mide la importancia relativa entre las fuerzas de flotación y el producto de las 
fuerzas viscosas y la conductividad térmica. La diferencia entre la expresión (1.1) y la escrita en (2.75) está en 
haber utilizado la aproximación 𝛽 ≅ 1/𝑇𝑚𝑒𝑑. En cambio, el número de Prandtl compara los efectos de la 
velocidad de difusión de cantidad de movimiento con los de velocidad de difusión de calor. 
Por último, antes de presentar las ecuaciones adimensionalizadas, se muestra la notación de los distintos 
operadores que intervienen en ella, también sin dimensiones: 
 
∂
∂𝑥
=
1
𝐻
∂
∂x∗
,
∂
∂z
=
1
𝐻
∂
∂z∗
, ∇=
1
𝐻
 ∇∗ , 
𝜕
𝜕𝑡
=
𝛼
𝐻2
𝜕
𝜕𝑡∗
 
(2.76) 
 
Adimensionalización de las ecuaciones 
Empleando los cambios de variable propuestos en (2.74) en las ecuaciones (2.62) – (2.63), y teniendo en cuenta 
(2.75) y (2.76), se tiene el siguiente