1. Educación

UNIVERSIDAD NACIONAL AUTONOMA DE MEXICO
PROGRAMA DE MAESTRÍA Y DOCTORADO EN
INGENIERÍA
FACULTAD DE INGENIERÍA
TRANSFERENCIA DE CALOR EN UN FLUIDO CONFINADO
ENTRE ESFERAS CONCÉNTRICAS CON ROTACIÓN
TESIS
QUE PARA OPTAR POR EL GRADO DE:
MAESTRO EN INGENIERÍA
INGENIERÍA MECÁNICA - TERMOFLUIDOS
P R E S E N T A:
ING. ARES CABELLO GONZALEZ
TUTOR:
DR. RUBEN AVILA RODRIGUEZ
2011
JURADO ASIGNADO
Presidente:
Dr. Cervantes De Gortari Jaime
Secretario:
Dr. Solorio Ordaz Francisco Javier
Vocal:
Dr. Ávila Rodríguez Rubén
1er. Suplente: Dr. Ramos Mora Eduardo
2do. Suplente: Dr. Cuevas García Sergio
Lugar donde se realizó la tesis:
DEPARTAMENTO DE TERMOFLUIDOS – FACULTAD DE INGENIERÍA
TUTOR DE TESIS:
Dr. RUBÉN ÁVILA RODRÍGUEZ
______________________________
FIRMA
´Indice general
´
Indice general
´
Indice de figuras
Agradecimientos
Nomenclatura
I
II
V
VII
Resumen
IX
Abstract
XI
1. Introducci´
on
1
2. Modelo f´ısico y modelo matem´
atico
9
2.1. Modelo F´ısico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.1.1. Estructura planetaria . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.1.2. Modelo propuesto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.2. Modelo Matem´atico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.2.1. Ecuaci´on de conservaci´on de masa . . . . . . . . . . . . . . . . . . . . . .
13
2.2.2. Ecuaci´on de conservaci´on de cantidad de movimiento . . . . . . . . . . . .
13
2.2.3. Ecuaci´on de conservaci´on de energ´ıa . . . . . . . . . . . . . . . . . . . . .
16
2.2.4. Par´
ametros adimensionales . . . . . . . . . . . . . . . . . . . . . . . . . .
17
3. Algoritmo num´
erico
21
3.1. M´etodo de los elementos espectrales . . . . . . . . . . . . . . . . . . . . . . . . .
21
3.1.1. Funciones de expansi´
on . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
3.1.2. Definici´
on de los puntos Gauss - Lobatto - Legendre en el dominio f´ısico .
23
3.1.3. M´etodo de Galerkin (Bubnov-Galerkin) . . . . . . . . . . . . . . . . . . .
24
´
INDICE GENERAL
II
3.1.4. Cuadratura de Gauss - Lobato - Legendre . . . . . . . . . . . . . . . . . .
25
3.2. Esfera c´
ubica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
4. Resultados
31
4.1. Conducci´on de calor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
4.2. Convecci´
on natural . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.2.1. N´
umero de Rayleigh cr´ıtico . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.2.2. Patrones de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
4.3. Convecci´
on natural y rotaci´
on . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
5. Conclusiones
53
A. Adimensionalizaci´
on de las ecuaciones
55
A.1. Ecuaci´on de continuidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
A.2. Ecuaci´on de cantidad de movimiento . . . . . . . . . . . . . . . . . . . . . . . . .
56
A.3. Ecuaci´on de la energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
B. Simplificaci´
on de los t´
erminos de rotaci´
on
59
C. C´
alculo del n´
umero de Nusselt para el caso difusivo.
63
Bibliograf´ıa
68
´Indice de figuras
2.1. Modelo f´ısico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
3.1. Ejemplo de la definici´on de los puntos dentro del dominio est´
andar . . . . . . . .
25
3.2. La esfera se divide en seis secciones las cuales corresponden a cada una de las
caras del cubo inscrito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
3.3. Representaci´
on de un elemento y su posici´
on dentro del macro elemento. . . . . .
28
3.4. Cuatro elementos del primer cuadrante de la cara interna del sector
. . . . . . .
28
3.5. Explosi´
on de los 6 (seis) elementos que forman la esfera . . . . . . . . . . . . . .
30
4.1. Comparaci´on entre la soluci´
on anal´ıtica y la soluci´
on num´erica . . . . . . . . . .
32
4.2. L´ıneas isotermas en los planos meridional y ecuatorial en el estado estacionario
del proceso de conducci´on de calor. . . . . . . . . . . . . . . . . . . . . . . . . . .
32
4.3. Isosuperficie de la temperatura promedio en el estado estacionario . . . . . . . .
33
4.4. Campos de velocidad para el n´
umero de Rayleigh cr´ıtico. . . . . . . . . . . . . .
34
4.5. L´ıneas isotermas en el plano meridional y ecuatorial (Ra(a) = 1606) . . . . . . .
35
4.6. Isosuperficies de temperatura (Ra(a) = 1606) . . . . . . . . . . . . . . . . . . . .
36
4.7. Patrones de convecci´
on en ´
anulos esf´ericos en el caso sin rotaci´
on T a = 0 y
Ra = 5 × 103 , visualizaci´
on del campo de temperatura en la direcci´
on radial. Las
sombras obscuras corresponden al flujo caliente ascendente y los colores brillantes
corresponden al flujo fr´ıo [Futterer et al., 2010]. . . . . . . . . . . . . . . . . . . .
37
4.8. Campo de temperatura en el hemisferio sur. Convecci´on para el caso sin rotaci´
on
T a = 0, Ra = 4 × 103 [Futterer et al., 2008]. . . . . . . . . . . . . . . . . . . . . .
38
4.9. N´
umeros de Nusselt en la convecci´on natural (Ra(a)=1606). . . . . . . . . . . . .
39
4.10. Isosuperficie de la temperatura en el inicio de la convecci´on (Ra(a) =1606, t∗ =0.2). 40
4.11. Isosuperficie de temperatura promedio en el pico del Nu externo (Ra(a) =1606,
t∗ =0.378).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
IV
´
INDICE DE FIGURAS
4.12. Isosuperficie de temperatura promedio en la convecci´
on natural (Ra(a) =1606,
t∗ =0.6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.13. Isosuperficie de temperatura promedio despues de modificarse el patr´
on convectivo (Ra(a) =1606, t∗ =1.183). . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.14. Isosuperficie de temperatura promedio una vez que los n´
umeros de Nusselt vuelven
a coincidir (Ra(a) =1606, t∗ =1.538). . . . . . . . . . . . . . . . . . . . . . . . . .
43
4.15. N´
umeros de Nusselt en la convecci´on natural (Ra(a) = 16500) . . . . . . . . . . .
43
4.16. Velocidad y temperatura promedio del caso 4 . . . . . . . . . . . . . . . . . . . .
45
4.17. Velocidad y temperatura promedio del caso 5 . . . . . . . . . . . . . . . . . . . .
46
4.18. Velocidad y temperatura promedio del caso 6 . . . . . . . . . . . . . . . . . . . .
48
4.19. Velocidad y temperatura promedio del caso 7 . . . . . . . . . . . . . . . . . . . .
49
4.20. Velocidad y temperatura promedio del caso 8 . . . . . . . . . . . . . . . . . . . .
50
4.21. Velocidad y temperatura promedio del caso 9 . . . . . . . . . . . . . . . . . . . .
51
Agradecimientos
Mi mayor agradecimiento es a mis padres V´ıctor y Rosa, quienes con su amor y su apoyo me
han permitido llegar hasta donde estoy. El significado de este escrito es un logro m´
as de su
trabajo.
A mis hermanos Ayax y Apolo, mis compa˜
neros y amigos de siempre, con quienes he
aprendido a compartir la vida misma.
A Selene, mi compa˜
nera incansable de incontables aventuras y de algunas otras desventuras,
con sus consejos y observaciones he aprendido a tomar decisiones en los momentos importantes.
´
Al doctor Rub´en Avila,
director de este trabajo, por su tiempo, dedicaci´
on y paciencia durante
mi formaci´
on profesional. Porque de ´el he aprendido otra forma de concebir la Ingenier´ıa.
A los doctores Jaime Cervantes, Sergio Cuevas, Eduardo Ramos y Francisco Solorio, que con
sus comentarios y ense˜
nanzas me han llevado a ser un mejor estudiante.
´
A la familia Avila
Mart´ınez, por su generosidad y su hospitalidad incondicional durante mi
estancia en la Universidad de California.
A los profes: Daniel, Ponce, Joaqu´ın y Conchita por su invaluable amistad, y su apoyo
desinteresado en los momentos dif´ıciles.
Al Consejo Nacional de Ciencia y Tecnolog´ıa por el apoyo econ´
omico brindado durante los
estudios de maestr´ıa.
Al personal del observatorio de visualizaci´on de la DGSCA, por su apoyo y dedicaci´
on en todo
momento.
A todos ellos mi m´
as sincero agradecimiento.
Nomenclatura
Cp
Calor espec´ıfico a presi´on constante [kJ/kg K]
Cg
Constante gravitacional [1/s2 ]
k
Conductividad t´ermica [W/m K]
∆T
Diferencial o incremento de temperatura [K]
E
Q˙
Energ´ıa total [J]
D
Longitud caracter´ıstica Re − Ri [m]
m
Masa [kg]
P
Presi´on [ Pa ]
P0
Presi´on est´
atica [Pa]
P∗
Presi´on din´
amica [Pa]
T
Temperatura [K]
T0
Temperatura inicial. Temperatura de referencia [K]
t
Tiempo [s]
W
Trabajo [W]
Re
Radio exterior [m]
Ri
Radio interior [m]
M
Vector cantidad de movimiento [kg/m s]
r
Vector de posici´
on de la part´ıcula [m]
F
Vector fuerza [N]
g
Vector gravedad [m/s2 ]
v
Vector velocidad [m/s]
V
Volumen [m3 ]
Flujo de calor por unidad de tiempo [W]
Nomenclatura
VIII
φ
´
Angulo
polar de las coordenadas esf´ericas [rad]
´
Angulo azimutal de las coordenadas esf´ericas [rad]
β
Coeficiente de expansi´
on t´ermica [1/K]
ω
Componente del vector velocidad angular [rad/s]
ρ
Densidad [kg/m3 ]
ρ0
Densidad de referencia [kg/m3 ]
α
Difusividad t´ermica [m2 /s]
δρ
Variaci´
on de la densidad [kg/m3 ]
Ω
Vector velocidad angular [rad/s]
ν
Viscosidad cinem´
atica [m2 /s]
µ
Viscosidad din´
amica [kg/m s]
D
Dt
d
dt
Derivada material o total
∇
Vector nabla
θ
Derivada o variaci´
on en el tiempo
Resumen
Estudios geol´
ogicos han permitido determinar que las capas internas de la Tierra est´
an formadas por una mezcla de materiales en estado l´ıquido o gaseoso a altas temperaturas, por lo
que su comportamiento en conjunto puede considerarse como el de un fluido viscoso. Adem´
as
los estudios y experimentos realizados han establecido que nuestro planeta en sus capas m´
as
internas est´
a formado en su mayor´ıa por metales pesados, y que en la capa m´
as profunda del
planeta la temperatura y la presi´on son tales que los materiales se encuentran en estado s´
olido.
Con esa idea se ha planteado el desarrollo de un modelo computacional que permita simular
la din´
amica interna del planeta bas´
andose en el estudio del material contenido en las capas
internas, haciendo uso para ello de las ecuaciones de la din´
amica de fluidos, lo cual permite
predecir no u
´nicamente el movimiento, sino tambien la distribuci´on de temperatura del fluido, y
con ello calcular la transferencia de calor desde el interior del fluido hacia los alrededores. Para
este modelo se considera una cavidad anular esf´erica con una relaci´
on de radios de 0.35, con
fronteras r´ıgidas a temperaturas constantes y diferentes, sin fuente interna de calor, un fluido
incompresible, viscoso, no conductor y sin deslizamiento en las fronteras.
En esta etapa de desarrollo del modelo se ha introducido el an´
alisis de la transferencia de
calor a trav´es de las fronteras, como funci´on de las tasas de rotaci´
on y de convecci´on natural
dentro de la cavidad. Considerando que son estos dos movimientos, junto con el campo magn´etico (no incluido en esta etapa), las fuerzas dominantes que provocan el movimiento del material
interno.
Las ecuaciones de Navier-Stokes tridimensionales para flujo incompresible se resuelven num´ericamente utilizando el m´etodo de elementos espectrales con polinomios Gauss-Lobatto-Legendre,
en una malla generada por medio del algoritmo de la esf´era c´
ubica que reproduce una geometr´ıa
Resumen
X
esf´erica en coordenadas cartesianas.
Se reportan tres casos de convecci´
on natural sin rotaci´
on, adem´
as para cada una de los casos
convectivos se simulan dos velocidades de rotaci´
on. Se presentan los patrones de velocidad y
temperatura y resultados de la transferencia de calor por ambas fronteras en los nueve casos.
De acuerdo con los resultados obtenidos en este trabajo de investigaci´
on es posible concluir
que el m´etodo num´erico permite determinar de manera satisfactoria el n´
umero de Rayleigh cr´ıtico en un sistema en donde la distribuci´on de la gravedad var´ıa inversamente con el radio de la
forma 1/r 5 . Se ha observado que la convecci´on natural no es un proceso estacionario y que la
rotaci´
on tiene un efecto estabilizante en el flujo que modifica la transferencia de calor tal y como
ha sido reportado en la literatura.
Una contribuci´on importante de esta tesis es el acoplamiento del m´etodo de elementos espectrales con el algoritmo de la esfera c´
ubica lo que ha permitido obtener resultados satifactorios,
mismos que pueden ser comparados con resultados previamente publicados.
Abstract
Several geological studies have shown the composition that the inner layers of the Earth’s,
this materials can be regarded as a high temperature viscous fluid. Furthermore, the studies and
experiments have established that our planet in its inner layers is formed mainly by heavy metals, and in the deepest layer, the temperature and pressure are such that the materials are solids.
Therefore in this research a computational numerical model to simulate the internal dynamics of the planet has been developed based on the study of the material in the most inner layers,
using the equations of fluid dynamics, which can predict the movement and the temperature
distribution of the fluid, also thus calculate the heat transfer from inside to the exterior of the
fluid. This model consists of a spherical annular cavity with a radius ratio of 0.35, with rigid
boundaries at constant and different temperatures, without internal heat source, an incompressible, viscous and non-conductive fluid with non-slip condition at the boundaries.
The three-dimensional Navier-Stokes equations are solved numerically using the spectral
element method (SEM) with Gauss-Lobatto-Legendre polynomials, in a mesh created by the
cubed-sphere algorithm in Cartesian coordinates.
We report three cases of natural convection without rotation, in addition for each convective
cases are simulated two rotation speeds. We present the velocity and temperature patterns and
results of heat transfer at the boundaries for each nine cases.
According to the results obtained in this research we conclude that the numerical method
allows satisfactorily determine the critical Rayleigh number in a system where the distribution
of gravity varies inversely with the radius of the form 1/r 5 . Also it was observed that natural
convection is not a stationary process and that the rotation has a stabilizing effect on the flow
Abstract
XII
changing the heat transfer as has been reported in the literature.
An important contribution of this work is the coupling of the spectral element method with
cubed sphere algorithm which allowed provide satisfactory results, they can be compared with
previously published results.
Cap´ıtulo 1
Introducci´
on
El fen´omeno de transferencia de calor a trav´es de un medio fluido est´
a presente en gran parte
de los procesos ingenieriles, desde la refrigeraci´
on de equipo electr´
onico por microcanales, hasta
los procesos de siderurgia y de fundici´
on de aleaciones. En la naturaleza tambi´en existen una
gran diversidad de fen´omenos de transferencia de calor en fluidos, desde el flujo capilar en los
alveolos pulmonares, hasta el efecto convectivo del sol. Es por ello que para los ingenieros es
muy importante conocer, estudiar y dominar los mecanismos por medio de los cuales se lleva a
cabo este transporte de energ´ıa.
En el caso particular de este trabajo de tesis, el fen´omeno que se estudia es el que se lleva
a cabo en las capas internas del planeta, con la finalidad de entender cuales son los precursores de los efectos que se observan en la superficie, como el campo magn´etico terrestre o la
disposici´
on geogr´
afica de los ejes volc´
anicos, para ello se propone un modelo num´erico que simula la capa exterior del n´
ucleo terrestre. Ese modelo toma en cuenta tanto las caracter´ısticas
f´ısicas, la geometr´ıa y las fuerzas involucradas, de las cuales se tiene evidencia como resultado
de estudios geol´
ogicos y sismol´ogicos [Jordan, 1979], [Stixrude et al., 1997], [Alfe et al., 2002],
[Monnereau et al., 2010].
Diversos trabajos anal´ıticos, experimentales y num´ericos se han publicado sobre la convecci´
on
natural y forzada de un fluido confinado entre esferas conc´entricas, en los cuales se considera
una gravedad que act´
ua sobre el sistema en direcci´
on radial, una rotaci´
on de cuerpo r´ıgido
(fuerzas de Coriolis y centr´ıfuga), una rotaci´
on diferencial (flujo de Couette esf´erico) y un campo magn´etico. [Roberts, 1968], [Busse, 1970], [Wimmer, 1978] [Marcus and Tuckerman, 1987],
[Dumas and Leonard, 1994], [Zhang, 1994], [Ardes et al., 1997], [Lorenzani and Tilgner, 2001],
[Travnikov et al., 2001], [Noir et al., 2001], [Simitev and Busse, 2003], [Zhang and Liao, 2004],
2
Introducci´
on
[Busse et al., 2003], [Hollerbach et al., 2006], [Evonuk and Glatzmaier, 2006], [Kelley et al., 2007],
[Tilgner, 2007], [Futterer et al., 2008], [Futterer et al., 2008], [Futterer et al., 2010], [Olson, 2011],
[Aubert et al., 2008], [Hernlund and Tackley, 2008], [Wu and Roberts, 2009], [King et al., 2009],
[Zhang et al., 2010].
Roberts(1968) desarroll´o curvas con los estados de convecci´on marginales para varias relaciones de n´
umeros de Taylor y Prandtl, incrementando el n´
umero de Rayleigh hasta obtener
movimientos convectivos. Encontr´
o que en una esfera de fluido con gravedad radial, rotaci´
on,
densidad y una distribuci´
on de las fuentes de calor, todas ellas uniformes, existen modos de
convecci´on que son asim´etricos con respecto al eje de rotaci´
on y adem´
as demostr´o que estos
modos de convecci´
on son los m´
as inestables aquellos donde el n´
umero de Taylor es peque˜
no.
Zhang(1994) propone que existen dos modos alternativos de convecci´on que son los m´
as
com´
unes a peque˜
nos n´
umeros de Prandtl, ´estos son escencialmente los modos de oscilaci´
on inercial de la ecuaci´
on de Poincar´e con la estructura m´
as simple a lo largo del eje de rotaci´
on y
simetr´ıa ecuatorial. Cada uno se propaga en direcci´
on opuesta al otro (este-oeste) y quedan
confinados en la regi´
on ecuatorial. Se plantea que el estudio de la frecuencia de estas oscilaciones
mediante soluciones de la ecuaci´
on de Poincar´e puede ser la base para desarrollar m´etodos de
estabilidad para fluidos en rotaci´
on. Sin embargo estas oscilaciones convectivas son estudiadas
con condiciones de frontera libres de esfuerzo.
Basados en la teor´ıa asint´
otica de Roberts(1968), y la teor´ıa de ondas de Zhang, Zhang &
Liao(2004) describen una nueva teor´ıa para la convecci´on en sistemas esf´ericos que giran r´
apidamente, n´
umero de Ekman muy peque˜
no (Ek << 1), basada en tres hip´
otesis. Primero, para
un n´
umero de Ekman peque˜
no, arbitrario pero fijo, la convecci´on localizada se extiende r´
apidamente en direcci´
on radial cilindrica y polar cil´ındrica con un n´
umero de Prandtl decreciente.
Segundo, se asume que para esas condiciones de n´
umero de Ekman, existe un flujo diferente de
cero en la capa l´ımite de Ekman, y tercero, que la velocidad dominate de la convecci´on se puede
expresar en funci´
on de la frecuencia media de la convecci´on t´ermica, del flujo en la capa l´ımite
de Ekman y de los modos de onda geostr´
oficos, para los cuales existen expresiones anal´ıticas.
Esta teor´ıa permitir´ıa generar expresiones anal´ıticas para la rotaci´
on diferencial generada por la
interacci´
on no-lineal de las ondas.
Aubert et. al.(2008) generan un modelo de la convecci´on termoqu´ımica, con la acci´on del
d´ınamo para encontrar explicaci´on a porqu´e el hemisferio este de la Tierra (40o − 180o E) es
3
s´ısmicamente m´
as r´
apido, m´
as isotr´
opico y m´
as atenuante que el hemisferio oeste. Este modelo
reproduce un n´
ucleo externo terrestre a gran escala y temporalmente muy largo considerando la
heterogeneidad del n´
ucleo interno y del bajo manto terrestre. Con este modelo logran reproducir
un ”viento termomagn´etico” cuya principal caracter´ıstica es una circulaci´
on cicl´onica por debajo
de Asia que concentra campo magn´etico en la frontera entre el manto y el n´
ucleo, lo que provoca
un despendimiento de material ligero en el hemisferio este de la frontera del n´
ucleo interno, lo
cual terminar´ıa generando la anomalia s´ısmica.
King et al.(2009) realizaron experimentos de laboratorio y simulaciones num´ericas de cavidades esf´ericas con rotaci´
on y convecci´
on, determinaron que los efectos de convecci´on- rotaci´
on
no est´
an definidos por el equilibrio global de fuerzas, concluyeron que existe una transici´on entre
los reg´ımenes de flujo dominados por la rotaci´
on y los dominados por la convecci´on y esta transici´on est´
a controlada por el espesor relativo entre la capa l´ımite de Ekman y la capa l´ımite t´ermica.
Dentro de los trabajos num´ericos, Busse(2003) y Simitev & Busse(2003) reportan los resultados de la simulaci´
on de un sistema de esferas conc´entricas con rotaci´
on y sometidas a un campo
gravitacional radial, y en el cu´
al la esfera interna tiene una fuente de calor. Con esta simulaci´
on
ellos encuentran los patrones convectivos dentro del ´anulo esf´erico y los comparan contra resultados experimentales. Debido a la forma de estos patrones Busse(1970) los llam´
o Banana cells
´o celdas o formaciones de banana .
Glatzmaier y Evonuk(2006) realizaron una simulaci´
on bidimensional (en el plano ecuatorial)
de la din´
amica interna de un planeta gaseoso con un n´
ucleo s´
olido, y analizaron la relaci´
on que
existe entre la velocidad angular con la que gira el planeta y la forma de los patrones convectivos formados debido a ese movimiento. Compararon las diferencias entre tener y no tener un
peque˜
no n´
ucleo s´
olido a diferentes velocidades angulares.
Una de las grandes complicaciones de la simulaci´
on num´erica de estos flujos, son las caracter´ısticas geom´etricas de la cavidad, que en principio pareciera ser totalmente adecuada para el
uso de coordenadas esf´ericas, sin embargo estas coordenadas tienen la particularidad de estar
num´ericamente indeterminadas en los polos, se han hecho varios intentos por resolver este inconveniente como el c´
alculo de aproximaciones sucesivas o la interpolaci´on entre los datos obtenidos
de los puntos cercanos a la indeterminaci´on. Una de las propuestas para resolver este problema
fu´e desarrollada por Hernlund(2008), con un modelo bidimensional en el cual se recogen todas
las caracter´ısticas tridimensionales del flujo mediante planos bisectores que atraviesan la esfera
4
Introducci´
on
y en los cuales queda retenida la informaci´
on de la tridimensionalidad, sin embargo este modelo
supone la existencia de una simetr´ıa de los resultados, la cual en general, no est´
a garantizada.
Experimentalmente, Futterer(2008) cre´
o un modelo en el cual se incorpora un sistema de
gravedad radial, inducida a trav´es de un campo el´ectrico de alto voltaje y rotaci´
on, con una diferencia de temperatura entre la frontera interior y la exterior. De este modelo obtienen campos
de velocidad y temperatura que utilizan para validar sus resultados obtenidos mediante simulaci´on num´erica, . De este estudio concluyen que existe una relaci´
on entre los estados convectivos
y los par´
ametros adimensionales propios del sistema. Adem´as determinan que la estabilidad depende mayormente del n´
umero de Taylor, es decir, la rotaci´
on genera un efecto estabilizante en
el flujo.
Otro trabajo experimental fu´e realizado por Egbers et. al.(2009) en la estaci´
on espacial
internacional, con un dispositivo autom´
atico capaz de reproducir diferentes condiciones de temperatura, gravedad y velocidad de rotaci´
on. Este experimento tiene la ventaja que ofrece la
condici´
on de micro-gravedad que existe en la estaci´
on espacial lo cual significa eliminar la gravedad laboratorio. Con los resultados obtenidos por interferometr´ıa se identifican los patrones
de flujo subcr´ıtico y supercr´ıtico, con ellos se elabora un mapa en el que se muestran los rangos
de los n´
umeros de Taylor y de Rayleigh en los que se presenta la primera inestabilidad del flujo.
Se han realizado trabajos experimentales para investigar la din´
amica de los fluidos del n´
ucleo
terrestre utilizando metales l´ıquidos y no metales, considerando la turbulencia, el arrrastre del
n´
ucleo interno, los efectos de la rotaci´
on (constante y variable), la convecci´on t´ermica, los efectos
qu´ımicos, la geometr´ıa de la cavidad, los campos magn´eticos y el cambio de fase, con el fin de
simular las condiciones de la din´
amica del n´
ucleo de la Tierra, sin embargo estos experimentos
presentan diversas dificultades, entre otras, la gran diferencia de escalas f´ısica y temporal entre
un experimento de laboratorio y el n´
ucleo real, la cantidad de fuerzas que est´
an interactuando
entre s´ı, como la flotaci´
on, viscosidad, electromagn´etica, interfacial, adem´
as de estar expuesto
a muchas aceleraciones, como la inercia, centr´ıfuga y Coriolis, aunado adem´
as a la restringida
cantidad de fluidos disponibles para hacer experimentos.
En una revisi´
on de los experimentos realizados para modelar la din´
amica del n´
ucleo planetario, Olson(2001) plantea que debido a lo anterior es que los experimentos de laboratorio
se han enfocado en tres aspectos principalmente. Primero la caracterizaci´
on de las propiedades de los materiales, que mezclados semejan las caracter´ısticas del n´
ucleo, como viscosidad,
5
difusividad t´ermica, difusividad magn´etica y propiedades qu´ımicas, ´esto permite tener una idea
m´
as aproximada de la composicion real del material interno. En segundo lugar, generar nuevos
descubrimientos, pues los resultados a menudo inesperados de experimentos de laboratorio proporcionan nuevas interpretaciones de las actuales observaciones geof´ısicas del n´
ucleo, y motivan
nuevos experimentos y observaciones. Por u
´ltimo, en los experimentos de laboratorio se pueden
proyectar las condiciones b´
asicas proporcionando datos a trav´es de un amplio rango de par´
ametros, que sirven como punto de referencia para el modelado num´erico.
El flujo de Couette esf´erico tambi´en se ha estudiado y documentado, [Hollerbach et al., 2006],
[Kelley et al., 2007], [Marcus and Tuckerman, 1987] este flujo se logra mediante la rotaci´
on diferencial de una de las esferas con respecto a la otra, y se ha propuesto que existen ”ondas”
generadas por el movimiento de rotaci´
on y que estas ondas est´
an fuertemente ligadas a las inestabilidades en el flujo.
Numericamente Hollerbach et. al.(2006) consideran un ´anulo esf´erico con la esfera exterior
fija y la interior rotando, con una relaci´
on de radios β = (ro − ri )/ri , en un rango entre 0.1 y
10. Encuentran un flujo base formado por un chorro desde la esfera interna hacia la externa y
que no depende de la relaci´
on de aspecto, este flujo se modifica a partir del n´
umero de Reynolds
(Rec ) que cambia en cada caso de acuerdo a la relaci´
on de aspecto. Estudian la inestabilidad
del flujo en todo el rango de β. Para 0.1 ≤ β ≤ 3.8 las inestabilidades tienen simetr´ıa opuesta
como el flujo base, y consisten de una serie de ondas en el chorro, el n´
umero de ondas decrece
desde 12 hasta 2, conforme la cavidad esf´erica se hace mas grande. Por arriba de esta relaci´
on
de radios, es decir ri /ro <0.2, la inestabilidad cambia de una perturbaci´on en el chorro y se
convierte en una variaci´
on en la intensidad del flujo de retorno despu´es de que ´este choca contra
la esfera exterior. Estos patrones de flujo coinciden con las simulaciones de Marcus & Tuckerman(1987) quienes adem´
as midieron el tama˜
no de los v´ortices de Taylor como funci´on del n´
umero
de Reynolds, cuyas mediciones concuerdan con las mediciones experimentales de Wimmer(1976).
Experimentalmente Keley et. al.(2007) generan un flujo de Couette en una cavidad esf´erica
con una relaci´
on de radios de 0.33 utilizando sodio l´ıquido, para medir la presencia de ondas
inerciales en el fluido en rotaci´
on por medio de inducci´on magn´etica. Las relaciones encontradas
entre los n´
umeros de onda del flujo, el campo magn´etico externo y el campo magn´etico inducido
concuerdan con lo propuesto por la teor´ıa, sin embargo las reglas que rigen los modos de inercia
tanto en el experimento como en los n´
ucleos planetarios siguen siendo desconocidas.
6
Introducci´
on
Finalmente otro aspecto que se ha estudiado en este tipo de modelos de la din´
amica interna es
el movimiento de precesi´on de la cavidad anular, [Lorenzani and Tilgner, 2001], [Noir et al., 2001],
[Tilgner, 2007], [Wu and Roberts, 2009], [Zhang et al., 2010] en diversos trabajos se ha encontrado que a´
un con una velocidad de precesi´on muy lenta, P o = Ωp /ω ∼ 10−2 existen efectos
notables sobre la capa de Ekman (Noir 2001), lo cual es visible al analizar los campos de vorticidad que se ven modificados por el efecto de la precesi´on (Lorenzani 2001), adem´
as el movimiento
de precesi´on genera una onda inercial con movimiento retrogrado, es decir en sentido contrario
al movimiento de precesi´on (Zhang 2001), (Zhang 2010). Todos estos movimientos tendr´ıan repercusiones en el flujo del n´
ucleo, y como consecuencia en el campo magn´etico inducido (Wu
2009), por lo que es conveniente que en el modelo que incorpora las ecuaciones de magnetohidrodin´
amica sea considerado el movimiento de precesi´on.
Adem´as de los trabajos de estudios de la estabilidad citados anteriormente, [Roberts, 1968],
[Zhang, 1994] y [Zhang and Liao, 2004], Travnikov(2001) hace un analisis de estabilidad de
energ´ıa, plantea que en el caso cuando la diferencia entre los radios de las esferas es muy
peque˜
no, (radio interno entre 0.9 y 1 veces el radio exterior) existe un par´
ametro definido como la raiz cuadrada del n´
umero de Grasoff, que permite determinar en funci´on del n´
umero de
Prandtl el punto cr´ıtico en el cual el flujo se vuelve inestable, sin embargo y de acuerdo a los
resultados planteados, ese par´
ametro de inestabilidad es v´alido u
´nicamente para n´
umeros de
Prandtl muy grandes ya que para el caso de n´
umeros de Prandtl de orden unidad el par´
ametro
de estabilidad difiere del analisis de estabilidad lineal.
En este trabajo de tesis se plantea un primer modelo muy simplificado de la din´
amica interna
del n´
ucleo de la Tierra, considerando una cavidad anular esf´erica con una relaci´
on de radios de
0.35, con fronteras r´ıgidas a temperaturas constantes y diferentes, sin fuente interna de calor. El
fluido utilizado es incompresible, viscoso, no conductor y con condiciones de no deslizamiento
en las fronteras. Se simulan tres condiciones de convecci´on; subcr´ıtico, marginalmente cr´ıtico y
supercr´ıtico y en cada una de ellas se simulan tres condiciones de rotaci´
on, incluyendo la velocidad de rotaci´
on igual a cero. Se resuelven las ecuaciones de la din´
amica de fluidos mediante el
m´etodo de elementos espectrales utilizando una malla generada con el algor´ıtmo de esfera c´
ubica.
En el cap´ıtulo 2 se presenten el modelo f´ısico y el modelo matem´
atico que rige al caso en
estudio.
En el cap´ıtulo 3 se presenta una breve descripci´on del m´etodo num´erico utilizado, asi como
7
el proceso de generaci´
on de la malla mediante el algor´ıtmo de la esfera c´
ubica.
En el cap´ıtulo 4 se presentan los resultados obtenidos para los nueve casos planteados, as´ı como los resultados de la validaci´
on del m´etodo.
Por u
´ltimo en el cap´ıtulo 5 se presentan las conclusiones del trabajo.
8
Introducci´
on
Objetivos
Para este trabajo de tesis se han planteado los siguientes objetivos:
Determinar el n´
umero de Rayleigh cr´ıtico considarando un campo de gravedad radial que
varia de la forma 1/r 5 .
Estudiar los patrones convectivos producidos por un campo gravitacional radial.
Estudiar el efecto de la rotaci´
on sobre los patrones y la transferencia de calor de la convecci´on natural.
Calcular la transferencia de calor a trav´es de las fronteras de un ´anulo esf´erico sometido a
diferentes relaciones de convecci´
on y rotaci´
on.
Visualizar y animar en realidad virtual los campos de velocidad, presi´on y temperatura.
Cap´ıtulo 2
Modelo f´ısico y modelo matem´
atico
Cualquier fen´omemo que pretenda ser estudiado y analizado, primeramente debe ser comprendido tomando en cuenta la f´ısica involucrada, ´esto con la finalidad de conocer la complejidad
y la viabilidad del estudio. El modelo f´ısico adem´
as permite conocer y tomar en cuenta todos los
principios f´ısicos involucrados, para que ´estos sean inclu´ıdos al momento de plantear los modelos
matem´
aticos que representan esa f´ısica.
En este cap´ıtulo se presenta la f´ısica involucrada en el fen´omeno natural, y posteriormente
las consideraciones tomadas en cuenta para el planteamieto de un modelo que reproduzca lo
m´
as posible al fen´omeno. Posteriormente se hace el desarrollo de cada una de estas ecuaciones
que constituyen el modelo matem´
atico a resolver y se plantean los par´
ametros adimensionales
para el an´
alisis de los flujos investigados en este trabajo.
2.1.
Modelo F´ısico
Como ya se ha dicho, es muy importante tener en cuenta la f´ısica del fen´omeno bajo estudio,
para ello es necesario situar de manera amplia el contexto en el cual se est´
a trabajando. Para
este trabajo se debe tener en cuenta las caracter´ısticas de la estructura interna de la Tierra, y
particularmente la f´ısica de la capa exterior del n´
ucleo.
2.1.1.
Estructura planetaria
Las diversas investigaciones desarrolladas para plantear modelos de la estructura interna
de la Tierra, confluyen en que nuestro planeta esta formado de diversas capas conc´entricas de
10
Modelo f´ısico y modelo matem´
atico
caracter´ısticas muy particulares cada una de ellas. Dichos modelos sugieren que la capa exterior (corteza terrestre) tiene un espesor aproximado de 350 km y esta dividida en dos capas,
la m´
as superficial llamada litosfera, de unos 100 km de espesor y con una composici´
on r´ıgida y
fr´agil. La capa profunda de la corteza se extiende entre los 100 y los 350 km, esta zona llamada
asten´
osfera es de una composici´
on pl´
astica y d´ebil. Por debajo de la corteza se encuentra la
mes´
osfera, tambien llamado manto terrestre, esta capa se extiende entre los 350 km y los 2800
km de profundidad [Jordan, 1979],[Monnereau et al., 2010].
M´as profundo se encuentra el n´
ucleo, constitu´ıdo por dos capas, una interna en estado s´
olido
y una externa en estado l´ıquido. La teor´ıa de la composici´
on interna de la Tierra propone que
el material que forma el n´
ucleo del planeta estar´ıa dado de la siguiente forma:
El n´
ucleo exterior l´ıquido formado de un material constitu´ıdo en su mayor parte por metales que por su densidad menor a la del n´
ucleo s´
olido permanecen en una capa de menor
presi´on que la interior, por lo cual poseen caracter´ısticas que lo asemejan a un l´ıquido
viscoso, este n´
ucleo exterior tiene un radio aproximado de 3488 km.
El n´
ucleo s´
olido formado por una mezcla de metales pesados que debido a su densidad
y al campo gravitacional planetario est´
an sometidos a un campo de presi´on tan grande
que los mantiene en estado s´
olido, este n´
ucleo s´
olido tiene un radio aproximado de 1231 km.
En el caso particular de este trabajo, la parte de la estructura interna de la Tierra a considerar es el n´
ucleo, por lo cual es necesario conocer algunas caracter´ısticas particulares de la
composici´
on del material y sus propiedades.
La sismolog´ıa tambi´en nos permite conocer los aspectos m´as sobresalientes de la estructura
interna de la Tierra, por medio del estudio de la velocidad de propagaci´on de las ondas s´ısmicas,
lo cual est´
a documentado para la mayoria de las capas externas, sin embargo el efecto de las
capas m´
as internas es todav´ıa una suposici´
on [Alfe et al., 2002].
Las t´ecnicas ab initio basadas en la teor´ıa funcional de la densidad permiten calcular los
potenciales qu´ımicos de los principales componentes de las impurezas tanto del n´
ucleo interno
como del n´
ucleo externo (S, O y Si), adem´
as de poder obtener estimaciones de las concentraciones presentes de cada uno de estos elementos. Los resultados muestran que el elemento con
2.1 Modelo F´ısico
mayor concentraci´
on es el ox´ıgeno, por encima del azufre y el silicio, lo cual sugiere que si la
convecci´on del n´
ucleo es la responsable del campo magn´etico, el ox´ıgeno es parte muy importante
de esta convecci´
on.[Alfe et al., 2002] [Stixrude et al., 1997]
Se ha propuesto tambi´en una teor´ıa que sostiene que no existe transferencia de masa entre el n´
ucleo externo y el manto inferior sino u
´nicamente un intercambio de energ´ıa, y que el
n´
ucleo s´
olido de la Tierra est´
a creciendo debido al aporte de material proveniente del n´
ucleo
exterior l´ıquido lo cual supondr´ıa modificar la raz´
on de aspecto entre las esferas que contienen
al fluido en estudio y tomar en cuenta tambi´en al material que compone el manto [Jordan, 1979].
Basado en lo anterior, es muy complicado conocer todas las caracter´ısticas de las capas profundas del planeta por lo que es necesario hacer suposiciones e hip´
otesis al respecto y trabajar
con ellas para obtener resultados.
2.1.2.
Modelo propuesto
Para poder simular un fen´omeno natural es necesario hacer algunas acotaciones, principalmente porque si se pretende simular un sistema tan complejo como el de la din´
amica interna
de un planeta, surgir´an muchas variables que no se pueden parametrizar y que, si se dejaran libres, la simulaci´
on requerir´ıa tantos recursos computacionales que ser´ıa imposible llevarla a cabo.
Es por eso que en este caso particular el estudio se centra sobre la din´
amica del n´
ucleo
externo, considerando su interacci´
on con los alrededores como condiciones de frontera sencillas
que permitan llevar a cabo simulaciones en tiempos de c´omputo considerables y con resultados
significativos.
Las condiciones tomadas en cuenta para este estudio est´
an basadas en la f´ısica del problema
anteriormente descrita, pero tambien tomando en cuenta las condiciones con las que se han
llevado a cabo experimentos y simulaciones anteriormente publicados, con la finalidad de tener
puntos de referencia en la literatura para poder comparar los resultados obtenidos.
En el modelo propuesto las condiciones de frontera son las siguientes: un n´
ucleo interno s´
olido, con una temperatura constante en su superficie, la frontera exterior r´ıgida a una temperatura
constante menor a la interior, sin transferencia de masa con los alrededores y con condiciones
de no deslizamiento para el fluido.
11
12
Modelo f´ısico y modelo matem´
atico
La cavidad bajo estudio, equivalente al n´
ucleo externo, se encuentre completamente llena de
un fluido newtoniano, viscoso, y con una densidad sujeta a la aproximaci´on de Boussinesq .
Figura 2.1: Modelo f´ısico
Es importante se˜
nalar que en el modelo propuesto no son considerados ni el manto, ni las
capas externas del planeta, por lo que la transferencia de calor calculada sobre la frontera exterior es equivalente a la energ´ıa que estar´ıa recibiendo el manto, por lo cual ser´ıa necesario,
posteriormente considerar, tanto la energ´ıa absorbida por el manto, la energ´ıa transmitida a la
corteza en forma de calor latente de cambio de fase, as´ı como la energ´ıa que se transmite a la
atm´
osfera, y en un estudio m´
as profundo y detallado considerar la energ´ıa radiante proveniente
del sol.
Los casos analizados en el presente trabajo de tesis se enlistan a continuaci´
on:
a) Conducci´on de calor
b) Convecci´
on natural
c) Convecci´
on natural con rotaci´
on
2.2 Modelo Matem´
atico
2.2.
13
Modelo Matem´
atico
El fluido confinado en la cavidad exsitente entre las esferas conc´entricas, se supone como un
fluido: newtoniano, incompresible, cuyas propiedades (densidad (ρ), viscosidad (µ), conductividad t´ermica (k), difusividad t´ermica (α), coeficiente de expansi´
on t´ermica (β)) son constantes.
Toda la cavidad gira con una velocidad angular constante (ω) alrededor del eje vertical (i3 ).
Las ecuaciones de la din´
amica de los fluidos est´
an basadas en el cumplimiento de tres principios fundamentales, conservaci´
on de masa, conservaci´on de cantidad de movimiento y conservaci´on de energ´ıa.
2.2.1.
Ecuaci´
on de conservaci´
on de masa
La ecuaci´
on de conservaci´
on de masa para un fluido incompresible esta dada como:
∇ · v = 0.
2.2.2.
(2.1)
Ecuaci´
on de conservaci´
on de cantidad de movimiento
La ecuaci´
on que describe el movimiento del fluido es la de conservaci´on de cantidad de movimiento que se deriva de la segunda ley de Newton,
F =
d(m · v)
,
dt
(2.2)
donde m es la masa y v es la velocidad.
Tomado en cuenta todas las fuerzas presentes en el movimiento del fluido, e igualandolas al
producto de la masa por la variaci´
on de la velocidad se tiene:
Fpresion + FCoriolis + Fcentrif uga + Fviscosa + Ff uente =
donde M es el vector de cantidad de movimiento.
DM
,
Dt
(2.3)
14
Modelo f´ısico y modelo matem´
atico
Puede demostrarse que la ecuaci´
on de cantidad de movimiento de un fluido incompresible
sujeto a las fuerzas mencionadas puede expresarse como:
ρ
donde
Dv
Dt
Dv
= −∇P − 2ρΩ × v − ρΩ × (Ω × r) − µ∇ × (∇ × v) + ρg,
Dt
(2.4)
es la derivada material de la velocidad, P es la presi´on, 2Ω × v es la aceleraci´
on de
Coriolis, Ω×(Ω×r) es la aceleraci´
on centr´ıfuga siendo Ω la velocidad angular, ∇×(∇×v) es la aceleraci´
on debido al efecto de la viscosidad y g es la aceleraci´
on de la gravedad [Greenspan, 1968].
Debido a que la variaci´
on de la densidad en el t´ermino ρg no puede ser despreciada, se utiliza
la aproximaci´on de Boussinesq, la cual dice que en todos los t´erminos de la ecuaci´
on de cantidad de movimiento la densidad es constante excepto en el t´ermino gravitacional, en el cual se
define como una funci´
on de la temperatura, lo cual es v´alido siempre y cuando la variaci´
on de la
densidad sea de alrededor del 1 %, para lo cual es necesario que la variaci´
on de la temperatura
no sea mayor a 10 grados. [Arpaci and Larsen, 1984]
En esta aproximaci´on tenemos entonces
ρ = ρ0 [1 − β(T − T0 )],
(2.5)
donde ρ0 es la densidad de referencia a la temperatura T0 y β es el coeficiente de expansi´
on
t´ermica.
La densidad entonces se escribe como una constante m´
as una variaci´
on, es decir: ρ0 + δρ, al
introducirlo al t´ermino fuerza de cuerpo de la ecuaci´
on de cantidad de movimiento (ec. (2.3)),
se escribe como:
1
[ρg] = 1/ρ0 [ρ0 + δρ]g,
ρ0
(2.6)
1
δρ
[ρg] = 1 +
g.
ρ0
ρ0
(2.7)
o bien
De la ecuaci´
on de estado (2.5) se tiene:
2.2 Modelo Matem´
atico
15
ρ − ρ0 = −ρ0 β(T − T0 ),
(2.8)
ρ − ρ0 = δρ,
(2.9)
δρ = −ρ0 β(T − T0 ).
(2.10)
donde
entonces:
Sustituyendo la nueva expresi´
on de la ecuaci´
on de estado en la ecuaci´
on de cantidad de movimiento, el t´ermino se reescribe como:
1
ρ0 β(T − T0
[ρg] = 1 −
g,
ρ0
ρ0
(2.11)
1
[ρg] = g − β(T − T0 )g,
ρ0
(2.12)
´o
por lo que la ecuaci´
on de cantidad de movimiento (ec. (2.3)) queda escrita entonces como:
ρ0
Dv
= −∇P + ρ0 g − ρ0 β(T − T0 )g − 2ρ0 Ω × v − ρ0 Ω × (Ω × r) − µ∇ × (∇ × v).
Dt
(2.13)
Si el t´ermino de presi´on se divide en dos, una presi´on est´
atica P0 m´
as una presi´on din´
amica
P∗
y se considera que en el equilibrio est´
atico (sin flotabilidad), las velocidades del fluido se
hacen cero, la ecuaci´
on de cantidad de movimiento (2.13) se reduce a:
−∇P0 + ρ0 g = 0.
(2.14)
16
Modelo f´ısico y modelo matem´
atico
Considerando esto, se puede reescribir la ecuaci´
on de cantidad de movimiento de la siguiente
manera:
ρ0
Dv
= −∇P ∗ + ρ0 β(T0 − T )g − 2ρ0 Ω × v − ρ0 Ω × (Ω × r) − µ∇ × (∇ × v).
Dt
(2.15)
Suponiendo que la gravedad est´
a dirigida hacia el centro de la esfera, el t´ermino de gravedad se
escribe entonces como:
F = −ρ0 β(T0 − T )gr,
(2.16)
donde T¯ es la temperatura promedio del fluido, g es la constante gravitacional y r es el vector
de posici´
on unitario en cada punto, definido por:
r=
r
.
|r|
(2.17)
En el caso particular que se estudia en este trabajo de tesis, el sistema gira alrededor del eje
vertical por lo que los t´erminos de rotaci´
on de la ecuaci´
on (2.15) se simplifican a una m´ınima
expresi´
on como se muestra en el ap´endice B.
2.2.3.
Ecuaci´
on de conservaci´
on de energ´ıa
Basado en el tercer principio fundamental de conservaci´on, la otra ecuaci´
on a tomar en cuenta es la de la energ´ıa que se que se deriva a partir de la primera ley de la termodin´
amica.
˙ = DE
Q˙ + W
Dt
(2.18)
˙ es el trabajo
donde Q˙ es la tasa de transferencia de calor que ingresa o se extrae al sistema, W
que se realiza sobre o por el sistema y E es la energ´ıa total del sistema.
En la ecuaci´
on de balance de energ´ıa (2.18) se deben considerar el transporte de energ´ıa hacia
y desde el volumen, la transferencia de calor por difusi´
on, el trabajo de las fuerzas de presi´on, el
2.2 Modelo Matem´
atico
17
trabajo de las fuerzas de cuerpo, el trabajo de las fuerzas viscosas y el t´ermino fuente de energ´ıa.
Tomando en cuenta lo anterior, y considerando que no existe fuente da calor interna, que la
conductividad t´ermica es constante y que el trabajo de las fuerzas viscosas es despreciable, la
ecuaci´
on de la energ´ıa se simplifica como:
DT
= α∇2 T
Dt
donde α =
2.2.4.
k
ρCp
(2.19)
es la difusividad t´ermica del fluido.
Par´
ametros adimensionales
En la din´
amica de los fluidos existen una serie de par´
ametros adimensionales que permiten
caracterizar a los flujos, estos par´
ametros son de suma importancia ya que gracias a ellos se
pueden hacer comparaciones o escalamientos de los modelos f´ısicos. Estos par´
ametros aparecen
a partir del proceso de adimensionalizaci´on de las ecuaciones como se muestra en el ap´endice A.
N´
umero de Rayleigh.
En el caso particular de este trabajo, es el par´
amertro que permite cuantificar la magnitud
del efecto convectivo sobre el fluido.
Debido a la gran variedad de definiciones de n´
umero de Rayleigh encontradas en la literatura,
en este trabajo se definen dos formas diferentes, una siguiendo el an´
alisis de Chandrasekhar(1961)
tomando como valor de referencia la gravedad en la esfera externa, y se describe a continuaci´
on:
Considerando que en el caso que se estudia la fuerza de gravedad est´
a orientada en la direcci´
on
radial, se puede escribir en coordenadas esf´ericas como:
fs = fr ir + 0iθ + 0iφ
(2.20)
donde fs es la fuerza en coordenadas esf´ericas, y fr es la componente de esa fuerza en la direcci´
on
radial.
Esa fuerza en coordenadas esf´ericas se puede transformar a sus componentes en coordenadas
cartesianas mediante una matriz de cambio de coordenadas A, mediante la ecuaci´
on fc = Afs ,
18
Modelo f´ısico y modelo matem´
atico
por lo tanto, la fuerza en coordenadas cartesianas se puede escribir como:

 

fr
senθcosφ cosθcosφ −senφ
fr senθcosφ

 






fc =  senθsenφ cosθsenφ cosφ   0  =  fr senθsenφ
cosθ
−senφ
0
0
fr cosθ


.

(2.21)
Por otro lado, de la geometr´ıa se pueden definir las funciones trigonom´etricas para los dos
´angulos, en funci´
on de las coordenadas cartersianas:
1
senθ =
(x21 +x22 ) 2
|r|
cosθ =
x3
|r|
senφ =
x2
1
(x21 +x22 ) 2
cosφ =
x1
(2.22)
1
(x21 +x22 ) 2
Sustituyendo las ecuaciones (2.22) en la matriz (2.21) se puede reescribir la fuerza en coordenadas cartesianas como:
fc = fr
✘
1
✘2✘
2+
✘x
2
(x✘
✘
1
2)
|r|
·
x1
✘
✘2✘21
2+
✘x
(x✘
✘
1
2)
o bien:
fc = fr
i1 + fr
✘
1
✘2✘
2+
✘x
2
(x✘
✘
1
2)
|r|
·
x2
✘
1
✘2✘
2
✘
2
(x✘
✘
1 + x2 )
i2 + fr
x2
x3
x1
i1 + fr
i2 + fr
i3 .
|r|
|r|
|r|
x3
,
|r|
(2.23)
(2.24)
Se tiene ya la formulaci´
on de la fuerza en coordenadas cartesianas, y esa fuerza de gravedad
m
depende del radio por lo que se puede escribir como fr = g(R0 ) | r |, cuyas unidades son 2 ,
s
y donde g(RO ) puede tomar cualquier forma.
Se puede entonces escribir la fuerza de gravedad de la forma:
x1
x2
x3
i1 + g(R0 )✚
| r✚|
i2 + g(R0 )✚
| r✚|
i3
| r✚|
| r✚|
| r✚|
✚
✚
✚
fc = g(R0 )✚
| r✚|
(2.25)
´o
fc = g(R0 )x1 i1 + g(R0 )x2 i2 + g(R0 )x3 i3
(2.26)
Cg
, donde Cg es una constante aleatoria con
r5
1
de tal manera que las unidades de g(R0 ) sigan siendo 2 .
s
Se propone entonces que g(R0 ) sea de la forma
unidades
m5
s2
Con esta formulaci´
on propuesta fr toma la forma:
fr =
Cg
·|r|
r5
(2.27)
como el valor m´
aximo de la fuerza de gravedad est´
a en el radio exterior, se puede evaluar
fr (R0 ) =
Cg
Cg
m5
m
·
R
=
≡
≡ 2 .
0
5
4
2
4
(R0 )
(R0 )
s m
s
(2.28)
2.2 Modelo Matem´
atico
19
Por u
´ltimo, sustituyendo esta formulaci´
on de la gravedad en la definici´on del n´
umero de
Rayleigh, ´este puede reescribirse como.
Ra =
Ra ≡
β∆T Cg D 3
ν 2 R04
1
m5 3
K K s2 m
2
2
m4 ms ms
(2.29)
≡[]
A esta definici´on del n´
umero de Rayleigh se le ha llamado Ra(a) para efectos de identificaci´on en el texto.
´
El otro an´
alisis, discutido ampliamente con el Dr. Rub´en Avila,
toma en cuenta la variaci´
on
de la gravedad respecto al radio, para posteriormente hacer una an´
alisis de estabilidad y se
obtiene una expresi´
on para el n´
umero de Rayleigh de la siguiente forma:
Ra =
ββCg
ναD
(2.30)
donde:
β=
Ti − To
1
1
ri − r0
A esta segunda definici´on se le ha llamado Ra(b) para su identificaci´
on.
N´
umero de Nusselt.
Es el par´
ametro que permite cuantificar la transferencia de calor a trav´es de las fronteras, y
se define como la relaci´
on entre el calor de convecci´on entre el calor por conducci´on:
Nu =
hD
k
(2.31)
donde h es el coeficiente convectivo, D es la distancia entre las dos esferas R0 − Ri , y k es la
conductividad t´ermica del fluido.
N´
umero de Taylor.
Es el par´
ametro que permite evaluar el efecto de las fuerzas de Coriolis y centr´ıfuga sobre las
fuerzas viscosas, y en este caso permite cuantificar el efecto de la rotaci´
on sobre la convecci´on.
Se define como:
Ta =
4ω 2 D 4
ν2
(2.32)
20
Modelo f´ısico y modelo matem´
atico
donde ω es la velocidad angular de giro, y ν es la viscosidad cinem´
atica del fluido.
N´
umero de Prandtl.
Es el par´
ametro que permite cuantificar la difusi´
on de cantidad de movimiento sobre la
difusi´
on de calor, y se define como:
Pr =
ν
α
(2.33)
donde ν es la viscosidad cinem´
atica, y α es la difusividad t´ermica del fluido.
De la adimensionalizaci´
on de las ecuaciones (2.1), (2.15) y (2.19) que se realiza detalladamente en el ap´endice A resulta el siguiente sistema de ecuaciones:
∇∗ · v ∗ = 0
Ta ∗
∂P ∗
Ta ∗
Dv ∗
1/2 ∗
1/2 ∗
i
+
T
a
v
−
i
=
−
−
T
a
v
+
r
r
+ RaT ∗ + ∇∗2 v.
1
2
1
2
Dt∗
4 1
4 2
∂x
1
Dθ ∗
=
∇∗2 T ∗ .
Dt∗
Pr
Este sistema se resuelve con el m´etodo descrito en el cap´ıtulo 3.
(2.34)
(2.35)
(2.36)
Cap´ıtulo 3
Algoritmo num´
erico
3.1.
M´
etodo de los elementos espectrales
El m´etodo num´erico utilizado para resolver las ecuaciones planteadas en el modelo matem´
atico es el de los elementos espectrales. Este m´etodo combina la generalidad del m´etodo
de los elementos finitos con la exactitud de las t´ecnicas espectrales. En la discretizaci´on de
los elementos espectrales, el dominio computacional es dividido en una serie de elementos, y
la velocidad en cada elemento es representada por un interpolante Lagrangiano de alto orden.
[Patera, 1983]
Los m´etodos espectrales consideran la expansi´
on de la soluci´
on de la ecuaci´
on diferencial en
una expansi´
on ortogonal de alto orden, donde los coeficientes de ´esta est´
an determinados por
una t´ecnica de proyecci´
on de residuos pesados. Las aproximaciones se vuelven de orden infinito
si las funciones de expansi´
on son seleccionadas apropiadamente.
El m´etodo de elemento finito es, en el sentido m´
as general, una t´ecnica de residuos pesados
aplicada a las series de expansi´
on, cada una con dominio sobre una peque˜
na regi´
on del espacio,
(elemento). Cuando la t´ecnica de los residuos pesados depende directamente del principio de
asociaci´on variacional, la continuidad de las condiciones de frontera natural est´
a impl´ıcitamente
satisfecha en las fronteras del elemento como parte del proceso de convergencia.
Al ser un m´etodo h´ıbrido recopila las ventajas de los dos que lo componen, por lo que se
puede tener una expansi´
on polinomial de la soluci´
on en cada uno de los elementos de la malla
asegurando la continuidad de la soluci´
on en las fronteras de cada elemento.
22
Algoritmo num´
erico
Una de las ventajas significativas de este m´etodo es la posibilidad de reducir considerablemente el n´
umero de elementos que componen la malla, debido a que la precisi´on del m´etodo se
puede modificar tambi´en incrementando el grado de los polinomios de interpolaci´on y viceversa,
esto se traduce directamente en un menor tiempo de c´omputo adem´
as de permitir la implementaci´on directa de algoritmos como el de la esf´era c´
ubica (cubed sphere) y su aplicaci´
on en
geometr´ıas complicadas.
3.1.1.
Funciones de expansi´
on
Al resolver ecuaciones diferenciales con m´etodos num´ericos, no es posible hallar la soluci´
on
exacta de la ecuaci´
on a resolver
L(u) = 0,
(3.1)
por lo que se debe hacer una aproximaci´on polinomial de la funci´on soluci´
on que se est´
a buscando
NT P
uδ = u0 +
ui (t)Φi (x),
(3.2)
i=1
y al hacer dicha aproximaci´on y sustituirla en la ecuaci´
on se genera un residual
L(uδ ) = R(uδ ),
(3.3)
que es equivalente al error en la soluci´
on. Por lo tanto el objetivo del m´etodo num´erico es hacer
que ese residual sea lo m´
as cercano a cero, es decir aproximar lo m´
as posible la soluci´
on num´erica
a la soluci´
on real.
L(uδ ) = R(uδ )
(3.4)
El m´etodo de los residuales pesados, define el producto interno entre funciones
f (x)g(x)dΩ = 0,
(f, g) =
(3.5)
Ω
con la finalidad de encontrar una funci´on cuyo producto interno con el residual sea igual a cero
(funci´on de peso), es decir, que el residual, pesado por la funci´on de peso sea igual a cero.
Entre mayor sea el n´
umero de puntos con los que se aproxima la soluci´
on, m´
as cercana est´
aa
la soluci´
on real, por lo que es v´alido decir que la soluci´
on aproximada es muy cercana a la soluci´on real conforme N T P → ∞ .
3.1 M´
etodo de los elementos espectrales
23
Existen diversos m´etodos para definir las funciones de peso. En el m´etodo de colocaci´
on
se definen las funciones de peso como funciones δ de Dirac (las funciones δ Dirac tienen las
siguientes propiedades
δ(x − a) = 0
x=a
∞
δ(x − a)dx = 1
∞
∞
δ(x − a)f (x)dx = f (a)),
(3.6)
∞
con lo que se garantiza que el residual sea cero en cada punto, es decir, que la funci´on es exacta
en cada punto de colocaci´
on.
Entonces el residual pesado se escribe como:
(δ(x − xj ), R(uδ ))dx = R(uδ (xj , t)) = 0.
(v(x)j , R(uδ )) =
Ω
En el m´etodo de volumenes finitos el dominio Ω se divide en subdominios y las funciones de
peso se definen como 1 dentro de su subdominio y 0 fuera de ´este. En este trabajo se utiliza el
m´etodo de Galerkin (secc. 3.1.3).
3.1.2.
Definici´
on de los puntos Gauss - Lobatto - Legendre en el dominio
f´ısico
Para este trabajo se definieron 384 (trescientos ochenta y cuatro) elementos de ´anulo esf´erico,
pero u
´nicamente los 4 (cuatro) elementos base de la malla fueron utilizados en el proceso de creacion de los puntos Gauss - Lobatto - Legendre, el resto se cre´o mediante el uso de coordenadas
esf´ericas (secci´
on 3.2). Cada uno de estos cuatro elementos de cre´
o a partir de 7 (siete) cubos,
proyectados hacia siete secciones de esfera conc´entricas, de las cuales se obtuvieron 343 puntos
dentro del dominio esf´erico que fueron transportados a un dominio est´
andar con un sistema de
coordenadas local (r, s, t) definido de -1 a 1 en cada uno de sus ejes.
Una vez en el dominio est´
andar se cre´
o una expansi´
on polinomial isoparam´etrica para la
posici´
on de cada punto definida como:
343
φi xi
x = φ1 x1 + φ2 x2 + φ3 x3 + ... =
i=1
24
Algoritmo num´
erico
donde x es la posici´
on del punto donde se est´
a calculando, φi es la funci´on de expansi´
on en cada
punto i, y xi es la posici´
on de cada uno de los 343 puntos de cada elemento definidos en el dominio
est´
andar 1 . La definici´on de las funciones de expansi´
on que se utilizan se trata en la secci´
on 3.1.3.
Lo anterior nos permite conocer la posici´
on de cualquier punto dentro del dominio est´
andar
siempre y cuando se conozca la posici´
on de los puntos base definidos, y como esos puntos est´
an
definidos en el dominio f´ısico, se conoce la posici´
on de cada uno. Como se trata de una expansi´
on
isoparam´etrica, se puede entonces conocer cualquier propiedad en cualquier punto del dominio.
Para acoplar esto al m´etodo de los elementos espectrales se buscan los puntos Gauss - Lobato
- Legendre (GLL) (seccion 3.1.4) dentro del dominio est´
andar, y se calculan con la expansi´
on y
asi se pueden mapear de regreso al dominio f´ısico de la esfera. En este caso se definen 63 = 216
puntos GLL para cada uno de los 384 elementos que componen la esfera.
3.1.3.
M´
etodo de Galerkin (Bubnov-Galerkin)
El m´etodo de los residuos pesados no define el tipo de funci´on de expansi´
on por lo que existen
diversas maneras de elegir estas funciones.
Una variaci´
on del m´etodo de Galerkin es el m´etodo tau, en el que v(x)j = φj (x), adem´
as
las funciones de expansi´
on no satisfacen las condiciones de frontera, ´estas son forzadas por una
serie de ecuaciones adicionales.
Otra variaci´
on del m´etodo es la aproximaci´on de Petrov - Galerkin, en la que las funciones de
´
peso v(x)j son diferentes de las funciones de expansi´
on φj (x). Este
es utilizado en los m´etodos
libres de mallas donde las funciones de expansi´
on (MLS), est´
an definidas por los puntos vecinos
mientras que las funciones de peso son splines de cuarto orden.
En los elementos espectrales se utiliza el m´etodo de Bubnov - Galerkin en el que las funciones
de peso son exactamente iguales a las funciones de expansi´
on φ, con lo que se tiene:
(φj (x), R(uδ ))dx = 0
(v(x)j , R(uδ )) = (φj (x), R(uδ )) =
Ω
1
En la figura 3.1 se muestra un ejemplo con 33 puntos, en el caso de este trabajo se definen 73 puntos para
cada elemento.
3.1 M´
etodo de los elementos espectrales
25
Figura 3.1: Ejemplo de la definici´on de los puntos dentro del dominio est´
andar
estas funciones est´
an definidas a partir de polinomios de Lagrange, lo cual garantiza la ortogonalidad entre ellas.
Adicionalmente se tiene que cada funci´on de expansi´
on satisface las condiciones de frontera
Bφj (x) = 0.
3.1.4.
Cuadratura de Gauss - Lobato - Legendre
Por lo descrito en la secci´
on 3.1.1, se tienen funciones que es necesario integrar sobre cada
dominio elemental, y se debe hacer de forma num´erica a trav´es de una suma finita de la forma:
Q−1
1
wi u(ξi )
u(ξ) dξ ≈
−1
i=0
26
Algoritmo num´
erico
donde wi son los pesos y ξi representan la absisa de Q diferentes puntos en el intervalo −1 ≤
ξi ≤ 1.
Dentro de las m´
ultiples formas de integraci´
on num´erica se encuentra la Cuadratura Gaussiana,
la cual define a las funciones u(ξ) como polinomios de Lagrange y a los pesos wi como integrales
de polinomios de Lagrange. En este caso en particular, se utiliza la cuadratura Gauss - Lobatto
- Legendre. Esta cuadratura emplea las absisas de los puntos que incluyan ambos extremos del
intervalo, es decir ξ = ±1. Entonces los pesos y las absisas se definen de la siguiente forma:
w10,0 =


i=0

 −1
1,1
ξi =
ξi−1,Q−2 i = 1, 2, 3, 4, ..., Q − 2


 1
i=Q−1
2
Q(Q − 1)[LQ−1 (ξi )]2
i = 0, 1, 2, ..., Q − 1
donde LQ (ξ) es el polinomio de Legendre LQ (ξ) = PQ0,0 (ξ).
3.2 Esfera c´
ubica
3.2.
Esfera c´
ubica
Al trabajar con dominios en geometr´ıas c´
ubicas surgen muchos problemas con el uso de
coordenadas esf´ericas ya que se tiene la comunmente llamada singularidad de los polos donde la
coordenada de la longitud puede tomar m´
ultiples valores, o se pueden presentar inconvenientes
con la transformaci´on de las ecuaciones a estas coordenadas. Una soluci´
on pr´
actica a este problema es el transformar el dominio esf´erico en un dominio c´
ubico en el cual se puedan resolver sin
mayor inconveniente las ecuaciones, eliminando el problema de las singularidades en los polos.
[Ronchi et al., 1996] [Nair et al., 2005]
El algoritmo de la esfera c´
ubica (cubed - sphere) es un m´etodo para mapear los puntos de
la superficie de una esfera sobre las caras de un cubo inscrito en la misma (ver Fig. 3.2 2 ). Para
lograr esto se definen una serie de transformaciones para cada uno de los sectores de la esfera el
cual se transladar´a a una de las caras del cubo.
Figura 3.2: La esfera se divide en seis secciones las cuales corresponden a cada una de las caras
del cubo inscrito
Una de las aportaciones m´
as importantes de este trabajo es el desarrollo de la malla num´erica, para ello se tom´
o el algoritmo original de seis sectores y se realiz´
o una subdivision de cada
uno de los macro sectores de la esfera, cada sector se dividi´o en 64 (sesenta y cuatro) elementos
cuya proyecci´
on corresponde a una secci´
on no sim´etrica de la pir´
amide trunca que representa a
cada sector (Fig. 3.3).
Para crear las coordenadas de los puntos de cada subsector, se utiliz´
o como punto de referencia la esquina inferior interior del primer macro sector y a partir de ahi se crearon 7 (siete)
2
Tomada de Ronchi et. al. 1996
27
28
Algoritmo num´
erico
Figura 3.3: Representaci´
on de un elemento y su posici´
on dentro del macro elemento.
sectores de esfera conc´entricos equidistantes, cada uno formado por una red de 625 (seiscientos
veinticinco) puntos, veinticinco en cada direcci´
on. Estos puntos se agruparon en cuatro zonas de
343 (trescientos cuarenta y tres) puntos, siete en cada direcci´
on, compartiendo los puntos de las
caras colindantes, cada unos de estos grupos corresponden a los cuatro elementos que forman el
primer cuadrante de la cara interna del primer sector (Fig. 3.4).
Figura 3.4: Cuatro elementos del primer cuadrante de la cara interna del sector
Cada uno de estos cuatro elementos se transladaron a un dominio estandard [-1,1], en el
cual se definen las raices del polinomio Gauss-Lobatto-Legendre en las tres direcciones, esto con
3.2 Esfera c´
ubica
la finalidad de calcular las coordenadas en el dominio f´ısico, de cada una de las raices de los
polinomios proyectados sobre el dominio estandard (secci´
on 3.1.2).
Teniendo las coordenadas de las raices de los polinimios (de orden 6), de cada uno de estos 4
elementos, se ordenan de manera que cada grupo de 216 puntos correspondan a cada elemento en
forma subsecuente. Estos 864 (ochocientos sesenta y cuatro) puntos, ya distribuidos correctamente y formando el primer cuadrante de la capa m´
as interna del primer elemento, se transformaron
a coordenadas esf´ericas, con la finalidad de aprovechar la simetr´ıa que estas coordenadas ofrecen.
Haciendo uso de las coordenadas esf´ericas se duplican los 4 elementos creados anteriormente,
en la direcci´
on colatitudinal o polar, posteriormente estos ocho elementos se duplican en direcci´
on
azimutal, cubriendo asi la capa m´
as interna del primer macro elemento, con dieciseis elementos
formados por 3456 (tres mil cuatrocientos cincuenta y seis) puntos. Una vez m´
as haciendo uso
de las coordenadas esf´ericas, esta primera capa se multiplica en la direcci´
on radial para obtener
4 capas conc´entricas. En este punto se ha construido el primer macro sector formado por 64
(sesenta y cuatro) elementos de 216 (doscientos dieciseis) puntos cada uno.
Este primer macro sector se multiplica tres veces en la direcci´
on colatitudinal hasta formar
el anillo de los cuatro macro elementos ecuatoriales de la esfera. Para los macro elementos de
los polos se debe hacer un tratamiento especial, ya que la posici´
on del primer elemento dentro
del macroelemento cambia para cada uno de los casquetes polares. Una vez ubicado el primer
elemento dentro de cada macroelemento polar, ´este se duplica primero en una direcci´
on, y posteriormente se duplica en la direcci´
on perpendicular, teniendo en consideraci´on que la cordenada
azimutal para esta replicaci´
on no crece, sino justo despu´es del polo vuelve a decrecer pero la
coordenada polar de esos puntos es una imagen a espejo con signos contrarios a los puntos que
la est´
an generando.
Una vez obtenidos los puntos que conforman a los seis macroelementos se deben corregir los
errores num´ericos en los puntos frontera de cada uno de los 384 (trescientos ochenta y cuatro)
elementos, para evitar que existan diferencias entre las coordenadas de los puntos que comparten
los elementos, posteriormente se crean los archivos de conectividad entre elementos, en los cuales
se especifica la vecindad de cada una de las seis caras de todos los elementos, es decir, se registra
qu´e cara y qu´e elemento comparten cada una de las seis fronteras de cada elemento. Una vez
completa la esfera, se transforman los 82944 (ochenta y dos mil novecientos cuarenta y cuatro)
puntos a sus coordenadas cartesianas, para crear el archivo de datos que ser´
a le´ıdo posteriormente
29
30
Algoritmo num´
erico
Figura 3.5: Explosi´
on de los 6 (seis) elementos que forman la esfera
por el preprocesador del programa de elementos espectrales para generar la malla computacional.
Cap´ıtulo 4
Resultados
En este trabajo de tesis se estudia primero, el proceso de conducci´on de calor entre dos
esferas conc´entricas, despu´es, el fen´omeno de convecci´on natural entre esas esferas conc´entricas
provocado por un campo de gravedad radial, posteriormente, se estudia el efecto de diferentes
intensidades de rotaci´
on sobre los patrones convectivos, haciendo en cada caso un an´
alisis de los
efectos de la convecci´
on y la rotaci´
on sobre la transferencia de calor.
Los resultados obtenidos se comparan con la soluci´
on anal´ıtica en el caso de la conducci´on
de calor, y con los resultados publicados del proyecto GEOFLOW en el caso de la convecci´
on
natural.
4.1.
Conducci´
on de calor
Cuando se realizan simulaciones num´ericas es necesario llevar a cabo una validaci´
on de los
datos obtenidos, ya sea con resultados experimentales o con soluciones anal´ıticas de las ecuaciones.
Para hacer la validaci´
on del c´
odigo computacional se realiz´
o una simulaci´
on del proceso puramente difusivo, del cual se conoce la soluci´
on anal´ıtica (ap´endice C), y se compararon tanto
la distribuci´
on de temperaturas como los coeficientes convectivos en las dos fronteras.
Al considerar el fen´omeno de flujo de calor en un fluido confinado entre dos esf´eras conc´entricas, con temperaturas uniformes constantes y diferentes en las fronteras, sin gravedad, no existe
convecci´on, es decir, no hay movimiento del fluido, por lo tanto la distribuci´on de temperatura
32
Resultados
en el estado estacionario depente u
´nicamente del radio, por lo que el proceso de conducci´on de
calor se vuelve unidimensional, sujeto u
´nicamente a las condiciones de frontera.
Figura 4.1: Comparaci´on entre la soluci´
on anal´ıtica y la soluci´
on num´erica
En la figura 4.1 se muestra la comparaci´on de la distribuci´on de temperatura adimensional
a lo largo del radio en el estado estacionario, donde T ∗ =
T −Te
Ti −Te
y r∗ =
r−Ri
Re −Ri .
En esta gr´
afica
se puede ver que la distribuci´
on de temperaturas obtenidas num´ericamente, corresponden con
la curva que representa la soluci´
on anal´ıtica de la temperatura.
(a) Isotermas meridionales
(b) Isotermas ecuatoriales
Figura 4.2: L´ıneas isotermas en los planos meridional y ecuatorial en el estado estacionario del
proceso de conducci´on de calor.
4.1 Conducci´
on de calor
Se realiz´
o tambien la validaci´
on del c´alculo num´erico del n´
umero de Nusselt (transferencia
de calor), tanto en la frontera interna como en la frontera externa, compar´andolos con los calculados anal´ıticamente a trav´es de las ecuaciones (C.26) y (C.33), desarrolladas en el ap´endice
C. Los resultados num´ericos reproducen el Nusselt unitario en cada frontera.
Se calcularon num´ericamente las componentes del vector flujo de calor a trav´es de ambas
fronteras y se transformaron a coordenadas esf´ericas mediante la matriz de cambio de coordenadas (ver ecuaci´
on (2.21)), con ´esto se demuestra que, en el caso difusivo, u
´nicamente se tiene
flujo de calor en la direcci´
on radial, y las componentes polar y azimutal son cero, lo cual es
f´ısicamente congruente con las isosuperficies de temperatura esf´ericas conc´entricas.
En la figura 4.2 se muestran las isol´ıneas de temperatura en los planos ecuatorial y meridional, y en la figura 4.3 se muestra la isosuperficie de la temperatura media del fluido, obtenida
mediante el promediado de las temperaturas de todos los puntos Gauss-Lobato-Legendre del
dominio computacional.
Figura 4.3: Isosuperficie de la temperatura promedio en el estado estacionario
Los resultados obtenidos en este proceso de validaci´
on muestran satisfactoriamente el acoplamiento del algoritmo de la esfera c´
ubica con el m´etodo de elementos espectrales.
33
34
Resultados
4.2.
4.2.1.
Convecci´
on natural
N´
umero de Rayleigh cr´ıtico
Para el estudio del fen´omeno de convecci´on natural entre esf´eras conc´entricas, se propone
un campo de gravedad radial que decrece con el radio en funci´
on de 1/r 5 , es decir, el valor de
la gravedad es menor en el radio exterior y es m´
as grande en el radio interno. Se utiliza esta
variaci´
on para hacerlo coincidir con el modelo utilizado en el experimento de micro gravedad
llevado a cabo en la estaci´
on espacial internacional [Futterer et al., 2010].
(a) Campo de velocidad del plano ecuatorial (b) Campo de velocidad del plano meridional
(Ra(a) = 1606)
(Ra(a) = 1606)
Figura 4.4: Campos de velocidad para el n´
umero de Rayleigh cr´ıtico.
Para poder observar los patrones convectivos, se busc´o computacionalmente el valor de la
constante gravitacional que iniciara la convecci´on natural en el fluido. Con esa constante se calcul´
o el n´
umero de Rayleigh cr´ıtico, y de acuerdo a los desarrollos planteados en la secci´
on 2.2.4,
los valores resultantes son Rac (a) = 1606 y Rac (b) = 7455.6.
No existen experimentos o simulaciones num´ericas con las caracter´ısticas establecidas en este
caso, temperaturas fijas en las fronteras, relaci´
on de radios 0.35 y gravedad con una variaci´
on
de la forma 1/r 5 .
Haciendo una revisi´
on de los casos con caracter´ısticas similares publicados en la literatura,
4.2 Convecci´
on natural
35
se encontr´
o que para el caso de una relaci´
on de radios 0.4 con una gravedad variando linealmente con el radio, el n´
umero de Rayleigh cr´ıtico definido de la misma forma que Ra(a) es
2847 [Tilgner, 1996]. Compar´andolo con los resultados obtenidos, el n´
umero de Rayleigh cr´ıtico
Rac (a) =1606 es cualitativamente correcto, ya que en este caso la fuerza de flotaci´
on es mayor,
lo cual implica que el n´
umero de Rayleigh cr´ıtico debe ser menor que en el caso de la variaci´
on
lineal de la fuerza.
Para el caso del experimento GEOFLOW con una variaci´
on de la gravedad de la forma 1/r 5
y una relaci´
on de radios ri /re = 0.5, el n´
umero de Rayleigh cr´ıtico es 2491 [Futterer et al., 2010].
Cabe mencionar que en la definici´on del n´
umero de Rayleigh utilizada en el experimento GEOFLOW se toman en cuenta los par´
ametros el´ectricos que generan la fuerza de gravedad. Compar´
andolo con este experimento, el resultado obtenido tambi´en es cualitativamente correcto dado
que en la simulaci´
on el radio interno es m´
as peque˜
no, por lo que la fuerza de flotaci´
on es mayor,
y el n´
umero de Rayleigh cr´ıtico debe ser menor. Este mismo comportamiento ha sido reportado
por Chandrasekhar(1961), es decir, entre mayor sea la relaci´
on de aspecto el n´
umero de Rayleigh
cr´ıtico es mayor [Chandrasekhar(1961), p. 244].
4.2.2.
Patrones de flujo
(a) Isotermas meridionales
(b) Isotermas ecuatoriales
Figura 4.5: L´ıneas isotermas en el plano meridional y ecuatorial (Ra(a) = 1606)
Al observar la figura 4.4(a) en la que se presenta un corte del campo de velocidades en el
plano meridional, para el n´
umero de Rayleigh cr´ıtico (Ra(a) =1606), y compararla con la figura
36
Resultados
4.4(b) se puede observar que existe una clara similitud en el campo de velocidades de los dos
planos, lo cual sugiere que el campo de velocidades es sim´etrico en la geometr´ıa esf´erica.
(a) Isosuperficie T ∗ = 0,25
(b) Isosuperficie T ∗ = 0,50
(c) Isosuperficie T ∗ = 0,65
Figura 4.6: Isosuperficies de temperatura (Ra(a) = 1606)
El efecto de la convecci´
on tambi´en se puede observar claramente en las l´ıneas isotermas (ver
fig. 4.5), las cuales tambi´en muestran simetr´ıa, tanto en el plano meridional como en el plano
ecuatorial, lo anterior tambi´en se puede corroborar al observar la figura 4.6 que presenta una
vista isom´etrica de tres isosuperficies de temperatura.
4.2 Convecci´
on natural
Figura 4.7: Patrones de convecci´
on en ´anulos esf´ericos en el caso sin rotaci´
on T a = 0 y
Ra = 5 × 103 , visualizaci´on del campo de temperatura en la direcci´
on radial. Las sombras
obscuras corresponden al flujo caliente ascendente y los colores brillantes corresponden al flujo
fr´ıo [Futterer et al., 2010].
Los patrones obtenidos para el rango de n´
umeros de Rayleigh 1606 < Ra <1.2×104 coinciden
cualitativamente con los mostrados en la figura 4.7 presentados por Futterer(2009). El n´
umero
de celdas obtenidas en este patr´
on de flujo, tambi´en corresponde con los patrones reportados
por Futterer(2008), (ver figura 4.8).
Por otro lado en diversos art´ıculos publicados en la literatura se hace menci´
on de la inestabilidad de la convecci´
on natural, por lo que se dej´o correr el tiempo de c´alculo sin modificar
los par´
ametros de la convecci´
on y se observ´
o que el estado convectivo encontrado permanece
durante un tiempo pero posteriormente se modifica, es decir no es un estado permamente.
Tambien se observ´
o que el tiempo que el patr´
on de flujo se mantiene, depende directamente
del n´
umero de Rayleigh, es decir entre m´
as grande sea la constante gravitacional menor es el
´
tiempo que se mantiene el patr´
on convectivo. Esto
sucede hasta un n´
umero Ra(a) de 1.2×104 ,
despu´es de este n´
umero ya no se obtiene el mismo patr´
on convectivo.
De los resultados obtenidos se puede ver tambi´en una clara dependencia de la transferencia
de calor a los efectos de la convecci´
on, esta dependencia puede ser observada tanto en el campo
de temperatura como en los n´
umeros de Nusselt tanto interno como externo.
37
38
Resultados
Figura 4.8: Campo de temperatura en el hemisferio sur. Convecci´on para el caso sin rotaci´
on
T a = 0, Ra = 4 × 103 [Futterer et al., 2008].
Una vez que se ha llegado al n´
umero de Rayleigh cr´ıtico, el fluido comienza a moverse, sin
embargo por efecto de la viscosidad, al inicio este movimiento se ve atenuado, hasta que por
efecto de la no-linealidad de las ecuaciones, la magnitud de la velocidad es tal que sobrepasa
al efecto viscoso, y el proceso convectivo se hace presente en la temperatura. Es decir, el inicio
del estado convectivo en la regi´
on cercana al n´
umero de Rayleigh cr´ıtico no es violento, es un
proceso lento debido a los efectos de la viscosidad.
´
Esto
se puede ver en la gr´
afica 4.9 donde, para el n´
umero de Rayleigh cr´ıtico ((Ra(a) =1606),
se observa que antes del tiempo adimensional 0.284 los n´
umeros de Nusselt son igual a 1, a partir
de este punto el fluido llega a un estado de movimiento tal que, la transferencia de calor cambia
y los n´
umeros de Nusselt crecen y oscilan. A partir del tiempo 0.473 la transferencia de calor se
estabiliza y se mantiene hasta el tiempo 0.994, despu´es del cual el patr´
on de flujo cambia y los
n´
umeros de Nusselt tambi´en se modifican.
Estos cambios en el n´
umero de Nusselt corresponden perfectamente con los cambios en el
patr´
on de flujo, y por consiguiente con las isotermas. Las figuras 4.10, 4.11, 4.12 y 4.13 muestran
las isosuperficies de la temperatura promedio para el n´
umero de Rayleigh cr´ıtico (Ra(a) =1606)
en los distintos puntos de la gr´
afica. La figura 4.10 corresponde a la zona en la cual el movi-
4.2 Convecci´
on natural
Figura 4.9: N´
umeros de Nusselt en la convecci´on natural (Ra(a)=1606).
miento del fluido a´
un no es suficiente como para modificar la transferencia de calor. La figura
4.11 corresponde al punto donde el Nusselt externo es m´
aximo, en t∗ = 0.378. La figura 4.12
representa el rango (0.48 < t∗ < 0.95). La figura 4.13 est´
a tomada en t∗ = 1.183 una vez que el
patr´
on convectivo se ha modificado. Finalmente la figura 4.14 muestra la temperatura promedio
en el tiempo t∗ = 1.538 una vez que los n´
umeros de Nusselt se han vuelto a igualar, sin embargo
despu´es de este punto el patr´
on convectivo no vuelve a ser estacionario.
La variaci´
on en la transferencia de calor se puede explicar observando la gr´
afica de los n´
umeros de Nusselt en el tiempo (fig. 4.9). En la primera parte el movimiento del fluido a´
un no es lo
suficientemente intenso como para modificar la transferencia de calor por conducci´on, raz´
on por
la cual el n´
umero de Nusselt se mantiene en la unidad. En el momento en que la velocidad del
fluido logra vencer los efectos viscosos, las capas de fluido comienzan a acelerarse y este movimiento provoca que se transfiera energ´ıa m´
as r´
apidamente, es por ´esto que el n´
umero de Nusselt
aumenta hasta un m´
aximo que se logra cuando la velocidad del fluido tambi´en es m´
axima. Despu´es de ´esto se llega al estado cuasi-estacionario del flujo, en este momento la transferencia de
calor se estabiliza, sin embargo se observa que los n´
umeros de Nusselt no son iguales, ´esto quiere
decir que por la esfera exterior se esta transmitiendo m´
as calor del que se est´
a transfiriendo
desde la esfera interna, ´esto se puede comprobar al observar que la temperatura media es menor
que la temperatura media en el caso difusivo.
39
40
Resultados
Figura 4.10: Isosuperficie de la temperatura en el inicio de la convecci´on (Ra(a) =1606, t∗ =0.2).
Posteriormente vuelve a existir un cambio de los patrones convectivos, y este cambio se ve
reflejado en el n´
umero de Nusselt, los cuales ahora tienden a coincidir nuevamente lo que significa
que las tasas de transferencia de calor tanto en la esfera interna como en la esfera externa vuelven
a ser las mismas, sin embargo no regresan hasta el valor que ten´ıan en el caso difusivo, ´esto es
debido a que el fluido se mantiene en movimiento y ´este aumenta la tasa de transferencia de calor.
La raz´
on por la que no se mantiene el estado cuasi-estacionario del flujo despu´es de despertar
la convecci´on, puede deberse a la diferencia en las tasas de transferencia de calor en las esferas
´
interna y externa. Esto
puede estar generando un desbalance t´ermico que al cabo del tiempo
provoca un nuevo brinco del patr´
on de flujo que buscar´a equilibrar la transferencia de calor por
las dos fronteras, sin embargo este nuevo patr´
on de flujo de ”equilibrio t´ermico” tampoco es un
patr´
on estacionario.
A diferencia de lo que ocurre en el estado marginal, cuando el n´
umero de Rayleigh se incrementa lo suficiente Ra >1.2×104 el inicio de la convecci´on es s´
ubito, y el cambio en los n´
umeros
de Nusselt est´
a dado por un escal´
on desde el inicio de la simulaci´
on, como se puede ver en la
gr´
afica 4.15, en la que se muestran los valores de los n´
umeros de Nusselt para un n´
umero de
Rayleigh de 16500.
4.2 Convecci´
on natural
Figura 4.11: Isosuperficie de temperatura promedio en el pico del Nu externo (Ra(a) =1606,
t∗ =0.378).
Note que para un n´
umero de Rayleigh de 16500, el n´
umero de Nusselt externo crece hasta
un valor cercano a 6, a diferencia del Rayleigh cr´ıtico, donde el valor no llega a 3. Adem´as para
el caso del Rayleigh supercr´ıtico, los n´
umeros de Nusselt oscilan y se estabilizan en un periodo
de tiempo muy corto.
41
42
Resultados
Figura 4.12: Isosuperficie de temperatura promedio en la convecci´on natural (Ra(a) =1606,
t∗ =0.6).
Figura 4.13: Isosuperficie de temperatura promedio despues de modificarse el patr´
on convectivo
(Ra(a) =1606, t∗ =1.183).
4.2 Convecci´
on natural
Figura 4.14: Isosuperficie de temperatura promedio una vez que los n´
umeros de Nusselt vuelven
a coincidir (Ra(a) =1606, t∗ =1.538).
Figura 4.15: N´
umeros de Nusselt en la convecci´on natural (Ra(a) = 16500)
43
44
Resultados
4.3.
Convecci´
on natural y rotaci´
on
Una vez obtenidos los patrones de la convecci´on natural, se estudi´o el efecto de la rotaci´
on
(fuerza de Coriolis) sobre estos patrones convectivos, para ello se tom´
o como punto de partida
el estado cuasi-estacionario de la convecci´on (Ra(a) =1606) y se consideraron dos diferentes
valores de velocidad angular con el fin de poder observar y evaluar los efectos de la intensidad
de esta nueva fuerza. En el primer caso, se aplic´
o una velocidad de rotaci´
on moderada (T a =
7.7×104 ) de tal manera que los efectos de la rotaci´
on y la convecci´on fueran equiparables, es
decir que los dos efectos compitieran entre s´ı. Despu´es se aplic´
o una velocidad de rotaci´
on grande
(T a = 1 × 106 ) dejando que el efecto de la rotaci´
on fuera dominante sobre la convecci´on natural.
Tambi´en se aplicaron estas mismas dos rotaciones a dos n´
umeros de Rayleigh m´
as, uno por
debajo del n´
umero de Rayleigh cr´ıtico (Ra(a) =1098) y otro supercr´ıtico ((Ra(a) =3570). Con
estos casos se tienen nueve condiciones, considerando los casos donde T a = 0. En estos casos se
calcula el n´
umero de Nusselt para ambas fronteras y se observan los patrones convectivos y las
superficies isotermas de cada caso. En la tabla 4.3 se muestran los valores del n´
umero de Nusselt
en cada frontera para cada uno de los casos.
Caso
Ta
Ra(a)
Ra(b)
N ui
N ue
1
0
1098
5097.8
1.003
0.9981
2
0
1606
7455.6
1.73
2.3551
3
0
3570
16568.
2.26
3.1214
4
77757
1098
5097.8
1.003
0.9974
5
77757
1606
7455.6
1.1133
1.1506
6
77757
3570
16568.
1.7096
1.5700
7
1000000
1098
5097.8
1.0037
1.0212
8
1000000
1606
7455.6
1.0058
1.0180
9
1000000
3570
16568.
1.058
1.0580
En la secci´
on 4.2.2 se describieron los casos sin rotaci´
on T a = 0, en esta secci´
on se muestran
algunos aspectos importantes de los seis casos con rotaci´
on.
En la figura 4.16 se observan los vectores de velocidad y la superficie de temperatura promedio del caso 4, se puede observar que a pesar de que el efecto de la rotaci´
on ya est´
a siendo presente
(plano meridional) a´
un no es suficiente para modificar el campo de temperaturas, adem´
as de
´
que al tener un n´
umero de Rayleigh subcr´ıtico la convecci´
on natural a´
un no est´
a presente. Esto
tambien puede observarse en los n´
umeros de Nusselt que permanecen muy cercanos a la unidad.
4.3 Convecci´
on natural y rotaci´
on
(a) Campo de velocidad del plano ecuatorial
45
(b) Campo de velocidad del plano meridional
(c) Isoterma de temperatura media
Figura 4.16: Velocidad y temperatura promedio del caso 4
En la figura 4.17 se observan los campos de velocidad y la temperatura promedio del caso 5,
en este caso ya existe una notoria presencia de las dos fuerzas, la de gravedad y la de rotaci´
on.
A pesar de que en el plano ecuatorial pareciera que persisten las celdas convectivas iniciales, al
observar el plano meridional y la isosuperficie de temperatura es clara la influencia de la rotaci´
on
y la modificaci´
on de las celdas convectivas. En cuanto a la transferencia de calor la influencia
de la rotaci´
on tambi´en es clara pues disminuy´o considerablemente el flujo de calor en direcci´
on
radial representado por los n´
umeros de Nusselt (ver tabla 4.3).
46
Resultados
(a) Campo de velocidad del plano ecuatorial
(b) Campo de velocidad del plano meridional
(c) Isoterma de temperatura media
Figura 4.17: Velocidad y temperatura promedio del caso 5
En la figura 4.18 se presentan los campos de velocidad y la superficie de temperatura promedio del caso 6. Aqu´ı se observa que, como la fuerza de gravedad es mayor, las celdas convectivas
tienden a ser m´
as marcadas, sin embargo por efecto de la rotaci´
on ´estas se ven deformadas y
son arrastradas en la direcci´
on de la rotaci´
on. Esto tambi´en se puede apreciar en los vectores del
plano meridional, los cuales est´
an dirigidos hacia la esfera interna en la zona cercana al ecuador,
por efecto de la fuerza de Coriolis. Al comparar los n´
umeros de Nusselt (tabla 4.3) es claro el
incremento de la convecci´
on con respecto al caso anterior (caso 5) pues la transferencia de calor
aumenta en la direcci´
on radial.
4.3 Convecci´
on natural y rotaci´
on
En las figuras 4.19, 4.20 y 4.21, se presentan los tres u
´ltimos casos (7, 8 y 9), y se observa
que la rotaci´
on es totalmente dominate, las celdas convectivas iniciales han desaparecido y se
´
observan los flujos hacia los polos encontrados en el trabajo previo [Cabello, 2009]. Unicamente
en el caso 9 se observa una ligera tendencia de la isosuperficie a recordar las celdas convectivas
por la magnitud del n´
umero de Rayleigh, sin embargo el efecto de la rotaci´
on es mucho m´
as
fuerte. La influencia de la rotaci´
on puede ser cuantificada en la transferencia de calor, que disminuye pr´
acticamente hasta el valor de difusi´
on (ver tabla 4.3).
Se observa que la rotaci´
on tiene un efecto estabilizante en el flujo pero a su vez disminuye la
tasa de transferencia de calor (ver tabla 4.3), esto puede deberse a que el flujo comienza a tener
componentes en las direcciones polar y azimutal, lo cual disminuye el transporte de energ´ıa en
direcci´
on radial.
Una vez incorporada la rotaci´
on se espera que no exista una nueva inestabilidad al incrementar la velocidad de giro, por el contrario, se espera que el flujo sea cada vez m´
as estable.
Es necesario realizar simulaciones con la fuerza de rotaci´on en una escala temporal muy larga,
con el fin de evaluar si la fuerza de Coriolis va a conducir a un estado de rotaci´
on de cuerpo r´ıgido.
De la figura 4.17 se puede observar que de los casos presentados, ´esta relaci´
on de n´
umeros
de Rayleigh y de Taylor es en la que las fuerzas est´
an m´
as equilibradas, es decir ninguna de
las dos es predominante sobre la otra, las dos aportan al movimiento final del fluido y como
consecuencia a la transferencia de calor. Es factible pensar que muy cercana a esta relaci´
on se
encuentra la raz´
on T a/Ra en la cu´
al los dos efectos est´
an plenamente acoplados, y que en ese
punto los patrones tengan una geometr´ıa particular distinta a las encontradas hasta ahora.
47
48
Resultados
(a) Campo de velocidad del plano ecuatorial
(b) Campo de velocidad del plano meridional
(c) Isoterma de temperatura media
Figura 4.18: Velocidad y temperatura promedio del caso 6
4.3 Convecci´
on natural y rotaci´
on
(a) Campo de velocidad del plano ecuatorial
49
(b) Campo de velocidad del plano meridional
(c) Isoterma de temperatura media
Figura 4.19: Velocidad y temperatura promedio del caso 7
50
Resultados
(a) Campo de velocidad del plano ecuatorial
(b) Campo de velocidad del plano meridional
(c) Isoterma de temperatura media
Figura 4.20: Velocidad y temperatura promedio del caso 8
4.3 Convecci´
on natural y rotaci´
on
(a) Campo de velocidad del plano ecuatorial
51
(b) Campo de velocidad del plano meridional
(c) Isoterma de temperatura media
Figura 4.21: Velocidad y temperatura promedio del caso 9
Cap´ıtulo 5
Conclusiones
Se logr´
o determinar de manera num´erica directa el n´
umero de Rayleigh cr´ıtico (Ra(a)c =
1606 (Ra(b)c = 7455)) para una distribuci´on de gravedad que varia inversamente con el radio de
la forma 1/r 5 , y que coincide cualitativamente con los n´
umeros de Rayleigh cr´ıticos publicados
en la literatura.
Se determin´o el rango del n´
umero de Rayleigh (1606 < Ra(a) < 12000 (7455 < Ra(b) <
56000)) en el cual la primera inestabilidad de la convecci´on natural es modo 2, sin embargo
los patrones convectivos permanecen durante un tiempo y despu´es se presenta otra instabilidad
cuyos patrones convectivos cambian, por lo que se concluye que no existe un estado estacionario.
Para n´
umeros de Rayleigh (Ra(a) > 12000 (Ra(b) > 56000)) existe otra bifurcaci´
on en la
cual el patr´
on convectivo y la transferencia de calor son diferentes al modo anterior.
Se encontraron bifurcaciones en los patrones de flujo que influyen directamente en la transferencia de calor, estas bifurcaciones dependen de la condici´
on inicial de la simulaci´
on. En ning´
un
caso se llega a un estado estacionario en los patrones de flujo, sin embargo la transferencia de
calor evaluada con el n´
umero de Nusselt tanto en la esfera interna como en la esfera externa
tiende a un estado estacionario.
Se observ´
o tambi´en que la rotaci´
on tiene un efecto estabilizante en el flujo, sin embargo decrece la tasa de transferencia de calor hasta el punto de regresar pr´
acticamente al estado difusivo.
Se desarroll´o una nueva malla computacional basada en el algoritmo de la esf´era c´
ubica que
permite obtener resultados con confiables, en un corto tiempo, y que son comparables satisfac-
54
Conclusiones
toriamente con resultados num´ericos y experimentales publicados.
En base a los resultados, podemos decir que el algoritmo utilizado reproduce satisfactoriamente las condiciones de movimiento y transferencia de calor de un fluido confinado en una
´
geometr´ıa esf´erica. Esto
nos puede permitir desarrollar modelos m´
as complicados, en los cuales
se puedan incluir los t´erminos magnetohidrodin´
amicos de las ecuaciones, as´ı como incluir mas
capas externas, para simular las condiciones del manto y de la atm´
osfera terrestres.
Ap´
endice A
Adimensionalizaci´
on de las
ecuaciones
A continuaci´
on se describe el procedimiento para adimensionalizar las ecuaciones
∇ · v = 0,
(A.1)
1
∂v
+ v · ∇v + 2Ω × v + Ω × (Ω × r) = − ∇P + g + ν∇2 v,
∂t
ρ
(A.2)
∂T
∂2T
∂2T
∂2T
+ v · ∇T = α[ 2 +
+
].
∂t
∂x1
∂x22
∂x23
Se definen los par´
ametros adimensionales.
v1∗ = νd v1 , x∗1 =
v2∗ = νd v2 , x∗2 =
v3∗ = νd v3 , x∗3 =
x1
d ,
x2
d ,
x3
d .
(A.3)
ν
t,
d2
T −T0
∗
T = ∆T ,
t∗ =
Despejando las variables dimensionales como funci´on de los par´
ametros adimensionales:
v1 =
v2 =
v3 =
v12 ν
d ,
v22 ν
d ,
v32 ν
d ,
d2 t∗
ν ,
= T ∗ ∆T
x1 = dx∗1 , t =
x2 = dx∗2 , T
x3 =
+ T0 ,
dx∗3 .
En el caso estudiado, el sistema gira en torno al eje z, por lo tanto, los t´erminos de las fuerzas
de Coriolis y centr´ıfuga se simplifican, como se muestra en el ap´endice B.
56
Adimensionalizaci´
on de las ecuaciones
A.1.
Ecuaci´
on de continuidad
∂v1
∂v2
∂v3
+
+
=0
∂x1 ∂x2 ∂x3
sustituyendo los par´
ametros adimensionales:
(A.4)
ν ∂v1∗
ν ∂v2∗
ν ∂v3∗
+
+
=0
∗
∗
dd ∂x1 dd ∂x2 dd ∂x∗3
(A.5)
∂v2∗
∂v3∗
ν ∂v1∗
+
+
=0
∗
∗
d2 ∂x1 ∂x2 ∂x∗3
(A.6)
∂v2∗
∂v3∗
∂v1∗
+
+
=0
∂x∗1 ∂x∗2 ∂x∗3
(A.7)
∇∗ · v ∗ = 0
(A.8)
o bien:
A.2.
Ecuaci´
on de cantidad de movimiento
∂v1
∂ 2 v1 ∂ 2 v1 ∂ 2 v1
∂v1
∂v1
∂v1
−1 ∂P
+v1
+v2
+v3
−2v2 ω3 −ω32 r1 =
+g1 +ν
i1 (A.9)
+
+
∂t
∂x1
∂x2
∂x3
ρ ∂x1
∂x21
∂x22
∂x23
∂v2
∂ 2 v2 ∂ 2 v2 ∂ 2 v2
∂v2
∂v2
∂v2
−1 ∂P
i2 (A.10)
+
+
+v1
+v2
+v3
+2v1 ω3 −ω32 r2 =
+g2 +ν
∂t
∂x1
∂x2
∂x3
ρ ∂x2
∂x21
∂x22
∂x23
∂v3
∂ 2 v3 ∂ 2 v3 ∂ 2 v3
∂v3
∂v3
∂v3
−1 ∂P
i3 .
+ v1
+ v2
+ v3
=
+ g3 + ν
+
+
∂t
∂x1
∂x2
∂x3
ρ ∂x3
∂x21
∂x22
∂x23
(A.11)
Desarrollando el t´ermino de la primera direcci´
on:
∂ 2 v1 ∂ 2 v1 ∂ 2 v1
∂v1
∂v1
∂v1
−1 ∂P
∂v1
+v1
+v2
+v3
−2v2 ω3 −ω32 r1 =
+g1 +ν
+
+
. (A.12)
∂t
∂x1
∂x2
∂x3
ρ ∂x1
∂x21
∂x22
∂x23
Sustituyendo los par´
ametros adimensionales:
∗
ν ν ∂v1
d2 d ∂t∗
∂v∗
∂v∗
∂v∗
1
2
3
1 1 ∂P
+ β(T − To )g1
ρ d ∂x∗1
ν ∂ 2 v1∗
ν ∂ 2 v1∗
ν ∂ 2 v1∗
+
+
,
+ν
dd2 ∂x∗2
dd2 ∂x∗2
dd2 ∂x∗2
1
2
3
ν
ν ∗ ν
ν ∗ ν
ν ∗
2
1
1
1
+ νd v1∗ dd
∂x∗ + d v2 dd ∂x∗ + d v3 dd ∂x∗ − 2 d v2 ω3 − ω3 r1 = −
A.2 Ecuaci´
on de cantidad de movimiento
57
∗
∗
∗
ν 2 ∂ 2 v1∗
1 ∂P
ν ∗
ν 2 ∂v1∗
∂ 2 v1∗
∂ 2 v1∗
∗ ∂v1
∗ ∂v1
∗
∗ ∂v1
2
+
v
+
v
+β(T
∆T
)g
+
+
v
v
ω
−ω
r
=
−
−2
+
+
,
1
3
1
2
3
1
2
3
d3 ∂t∗
∂x∗1
∂x∗2
∂x∗3
d
ρd ∂x∗1
d3 ∂x∗2
∂x∗2
∂x∗2
1
2
3
(A.13)
∗
∗
∗
d3
d3 2
d2 ∂P
∂v1∗
d2 ∗
∗ ∂v1
∗ ∂v1
∗ ∂v1
+
v
+
v
+
β
+
v
v
ω
−
ω
r
=
−
∆T g1 T ∗ + ∇∗2 v1∗ .
−
2
3
1
2
3
1
2
3
∗
∗
∗
∗
∂t∗
∂x1
∂x2
∂x3
ν
ν2
ρν 2 ∂x1
ν2
(A.14)
Si se define una presi´on adimensional y un radio adimensional:
P∗ =
P
ρν 2
d2
r1∗ =
r1
d
se puede reescribir la ecuaci´
on (A.14) como:
Dv1∗ 2d2 ω3 ∗ d4 ω32 ∗
∂P ∗ βd3 ∆T g1 ∗
+
−
v
−
r
=
−
T + ∇∗2 v1∗ ,
2
Dt∗
ν
ν2 1
∂x∗1
ν2
(A.15)
y recordando la definici´on de los n´
umeros adimensionales:
Ta =
4Ω2 d4
ν2
Ra =
βd3 ∆T g
ν2
la ecuaci´
on (A.15) finalmente queda escrita de la siguiente manera:
1
Dv1∗
∂P ∗
Ta ∗
∗
2v −
+ RaT ∗ + ∇∗2 v1∗ .
−
T
a
r
=
−
2
Dt∗
4 1
∂x∗1
(A.16)
En la segunda direcci´
on:
∂ 2 v2 ∂ 2 v2 ∂ 2 v2
∂v2
∂v2
∂v2
−1 ∂P
∂v2
+v1
+v2
+v3
+2v1 ω3 −ω32 r2 =
+g2 +ν
+
+
. (A.17)
∂t
∂x1
∂x2
∂x3
ρ ∂x2
∂x21
∂x22
∂x23
Se reescribe la ecuaci´
on (A.17) como:
∂P ∗ βd3 ∆T g2 ∗
Dv2∗ 2d2 ω3 ∗ d4 ω32 ∗
+
v
−
+
r
=
−
T + ∇∗2 v2∗ ,
1
2
Dt∗
ν
ν2
∂x∗2
ν2
(A.18)
1
Dv2∗
∂P ∗
Ta ∗
∗
2v −
+ RaT ∗ + ∇∗2 v2∗ .
+
T
a
r
=
−
2
1
Dt∗
4
∂x∗2
(A.19)
o bien:
En la tercera direcci´
on:
∂ 2 v3 ∂ 2 v3 ∂ 2 v3
∂v3
∂v3
∂v3
−1 ∂P
∂v3
+
+
,
+ v1
+ v2
+ v3
=
+ g3 + ν
∂t
∂x1
∂x2
∂x3
ρ ∂x3
∂x21
∂x22
∂x23
(A.20)
58
Adimensionalizaci´
on de las ecuaciones
´o:
Dv3∗
∂P ∗
+ RaT ∗ + ∇∗2 v3∗ .
=
−
Dt∗
∂x∗3
A.3.
(A.21)
Ecuaci´
on de la energ´ıa
∂2T
∂2T
∂T
∂2T
∂T
∂T
∂T
+
+
.
+ v1
+ v2
+ v3
= αl
∂t
∂x1
∂x2
∂x3
∂x21
∂x22
∂x23
(A.22)
Sustituyendo los par´
ametros adimensionales:
∆T ν ∂T ∗ ν ∗ ∆T ∂T ∗ ν ∗ ∆T ∂T ∗ ν ∗ ∆T ∂T ∗
∆T ∂ 2 T ∗ ∆T ∂ 2 T ∗ ∆T ∂ 2 T ∗
+ 2
+ 2
,
+
+
=
α
+
v
v
v
d2 ∂t∗ d 1 d ∂x∗2 d 2 d ∂x∗2 d 3 d ∂x∗3
d2 ∂x∗2
d ∂x∗2
d ∂x∗2
1
2
3
(A.23)
ν
∗
∗
∗
∂T ∗
∂2T ∗ ∂2T ∗ ∂2T ∗
∗ ∂T
∗ ∂T
∗ ∂T
+
+
,
+
v
+
v
+
v
=
α
2
3
1
∂t∗
∂x∗1
∂x∗2
∂x∗3
∂x∗2
∂x∗2
∂x∗2
1
2
3
(A.24)
o bien:
∗
∗
∗
∂T ∗
α ∂2T ∗ ∂2T ∗ ∂2T ∗
∗ ∂T
∗ ∂T
∗ ∂T
+
v
+
v
=
+
v
+
+
.
2
3
1
∂t∗
∂x∗1
∂x∗2
∂x∗3
ν ∂x∗2
∂x∗2
∂x∗2
1
2
3
(A.25)
Retomando la definici´on del n´
umero de Prandtl:
∗
∗
∗
1 ∂2T ∗ ∂2T ∗ ∂2T ∗
∂T ∗
∗ ∂T
∗ ∂T
∗ ∂T
+
v
+
v
=
+
v
+
+
,
2
3
1
∂t∗
∂x∗1
∂x∗2
∂x∗3
P r ∂x∗2
∂x∗2
∂x∗2
1
2
3
(A.26)
o de forma vectorial:
1
DT ∗
=
∇∗2 T ∗ .
Dt∗
Pr
(A.27)
Ap´
endice B
Simplificaci´
on de los t´
erminos de
rotaci´
on
Debido a que el sistema gira alrededor del eje vertical, los t´erminos Ω × (Ω × r) y 2Ω × v de
la ecuaci´
on (2.15) se simplifican de la siguiente manera:
Desarrollando en notaci´
on ´ındice el t´ermino de aceleraci´
on centr´ıfuga se tiene:
Ω × (Ω × r) = ωi ωj ri ij − ωi ωi rj ij
= ω1 ωj r1 ij + ω2 ωj r2 ij + ω3 ωj r3 ij − ω1 ω1 rj ij − ω2 ω2 rj ij − ω3 ω3 rj ij
= (ω1 r1 ω1 + ω2 r2 ω1 + ω3 r3 ω1 − ω1 ω1 r1 − ω2 ω2 r1 − ω3 ω3 r1 )i1
+(ω1 r1 ω2 + ω2 r2 ω2 + ω3 r3 ω2 − ω1 ω1 r2 − ω2 ω2 r2 − ω3 ω3 r2 )i2
+(ω1 r1 ω3 + ω2 r2 ω3 + ω3 r3 ω3 − ω1 ω1 r3 − ω2 ω2 r3 − ω3 ω3 r3 )i3 .
De igual manera para el t´ermino de aceleraci´
on de Coriolis:
(B.1)
60
Simplificaci´
on de los t´
erminos de rotaci´
on
2Ω × v = 2ǫijk ωi vj ik
= 2ǫ1jk ω1 vj ik + 2ǫ2jk ω2 vj ik + 2ǫ3jk ω3 vj ik
= 2ǫ11k ω1 v1 ik + 2ǫ12k ω1 v2 ik + 2ǫ13k ω1 v3 ik
+2ǫ21k ω2 v1 ik + 2ǫ22k ω2 v2 ik + 2ǫ23k ω2 v3 ik
+2ǫ31k ω3 v1 ik + 2ǫ32k ω3 v2 ik + 2ǫ33k ω3 v3 ik
= (2ω2 v3 − 2ω3 v2 )i1
+(2ω3 v1 − 2ω1 v3 )i2
+(2ω1 v2 − 2ω2 v1 )i3 .
(B.2)
Considerando que el sistema gira con una velocidad constante en torno al eje x3 (z), es decir,
que las componentes del vector velocidad angular Ω en la direcci´
on x1 (x) y en la direcci´
on x2
(y) son iguales a cero y la componente en la direcci´
on z tiene un valor constante; entonces se
puede reescribir el t´ermino de la fuerza centr´ıfuga de la siguiente forma:
(ω3 ω3 r1 )i1 + (ω3 ω3 r2 )i2 .
(B.3)
De manera semejante el t´ermino de Coriolis se escribe como:
(−2ω3 v2 )i1 + (2ω3 v1 )i2 .
(B.4)
Para agregar estas fuerzas al programa de c´omputo, los t´erminos de las fuerzas centr´ıfuga
(ec. (B.3)) y de Coriolis (ec. (B.4)) se escribieron de la siguiente forma:
Fx = ρ0 ω 2 x
Fy = ρ0 ω 2 y
centrif uga,
(B.5)
61
Fx = −2ρ0 ωvy
Fy = 2ρ0 ωvx
Coriolis
(B.6)
donde la ω el la velocidad angular con la que gira el sistema en torno al eje vertical (z) , vx
y vy son las componentes en x y y de la velocidad en cada punto, respectivamente.
Ap´
endice C
C´
alculo del n´
umero de Nusselt para
el caso difusivo.
Para calcular el n´
umero de Nusselt es necesario conocer la distribucii´on de temepraturas a
lo largo del radio para el caso difusivo, para ello se parte de la ecuaci´
on diferencial de la energ´ıa
en coordenadas esf´ericas:
∂
∂r
r2
∂T
∂r
=0
(C.1)
cuya soluci´
on general es:
A
+B
r
(C.2)
T1 en r = r1
(C.3)
T (r) =
sujeta a las condiciones de frontera,
T2 en r = r2
y se obtiene el sistema de ecuaciones,
A
+B
r1
A
+B
T2 =
r2
T1 =
(C.4)
de donde se puede despejar a la constante B de la primera ecuaci´
on,
B = T1 −
A
r1
(C.5)
64
C´
alculo del n´
umero de Nusselt para el caso difusivo.
y substituirla en la segunda ecuaci´
on,
T2 =
A
A
+ T1 −
r2
r1
A
A
−
r2 r1
=
+ T1
(C.6)
simplificando,
1
1
−
r2 r1
T2 = A
+ T1
(C.7)
y reacomodando los t´erminos se puede escribir como:
A
1
1
−
r2 r1
= T2 − T1
(C.8)
con esto se obtiene una expresi´
on para la constante A,
A=
T2 − T1
1
r2
−
1
r1
.
(C.9)
Posteriormente, sustituyendo la ecuaci´
on (C.9) en la ecuacion (C.5), se obtiene:
B = T1 −


1  T2 − T1 
.
1
1
r1
r2 − r1
(C.10)
Con esto se tiene una expresi´
on de la constante B en t´erminos de los radios y las temperaturas,

B = T1 − 

T2 − T1 
.
r1
r2 − 1
(C.11)
Sustituyendo las ecuaciones (C.9) y (C.11) en la ecuaci´
on (C.2), se obtiene,
T (r) =




1  T2 − T1 
T2 − T1 
+ T1 − 
r1
1
1
r
r2 − r1
r2 − 1
(C.12)
que es la distribuci´
on de la temperatura como funci´on del radio para todo el ´anulo.
Para hacer el c´
alculo anal´ıtico del n´
umero de Nusselt en las fronteras interior y exterior del
´anulo, es necesario definir algunas variables basadas en las dimensiones del modelo f´ısico:
r1 = rinterior , r2 = rexterior
D = r2 − r1 ,
∆T = Tinterior − Texterior .
Una vez definidas las nuevas variables, el an´
alisis parte de la definici´on del n´
umero de Nusselt:
65
Nu =
hD
k
(C.13)
y se sustituye la expresi´
on del coeficiente convectivo,
Nu =
✁k
− ∆T
∂T
∂r r=R
D
=−
k
✓
1 ∂T
∆T ∂r
D.
(C.14)
r=R
Si se eval´
ua la derivada,
∂T
∂r
se simplifica en:
r=R



 



∂ 1
T2 − T1
T2 − T1 
+ T1 −
=
 r1 − 1 
∂r r  1 − 1 
r2
r1
r2
∂T
∂r
=−
r=R
1
r2
T2 − T1
1
1
r2 − r1
(C.15)
r=R
.
(C.16)
El n´
umero de Nusselt queda entonces como:
Nu = −
1
∆T
−
1
r2
T2 − T1
1
1
r2 − r1
(r2 − r1 ) .
(C.17)
Retomando la expresi´
on (C.17) y tomando en cuenta el radio en la superficie interna, (r = r1 ),
el n´
umero de Nusselt queda como:
N uint = −
1
∆T
−
T2 − T1
1
1
r2 − r1
1
r12
(r2 − r1 )
(C.18)
considerando que si ∆T = T1 − T2 , entonces −∆T = T2 − T1 , por lo que:
✟)
(−✟
∆T
1
1
r2 − r1
1
1
N uint = − ✟ − 2
∆T
r1
✟
(r2 − r1 )
(C.19)
simplificando los signos:
N uint = −
1
r12
1
r2
1
−
1
r1
(r2 − r1 )
(C.20)
simplificando los parentesis:
N uint =
−1
r12
1
r2
−
1
r1
realizando la multiplicaci´
on en el denominador:
(r2 − r1 )
(C.21)
66
C´
alculo del n´
umero de Nusselt para el caso difusivo.
N uint =
−1
r12
r2
− r1
(r2 − r1 )
(C.22)
(r2 − r1 )
(C.23)
obteniedo el com´
un denominador:
−1
N uint =
r12 −r1 r2
r2
el n´
umero de Nusselt en la frontera interna queda como:
N uint =
r12
−r2
(r2 − r1 )
− r1 r2
(C.24)
si se multiplica esta expresi´
on del n´
umero de Nusselt por la relaci´
on adimensional de los
radios
r1
r2 ,
el n´
umero de Nusselt interno queda como:
N uint =
−r2
(r2 − r1 )
r12 − r1 r2
r1
r2
✘
=
✘r✘
r2✘
(r✘
1−
2)
✘
✘
✘
r1✘
(r✘
−
r
1
2)
r1
r2
(C.25)
que se reduce a
N uint = 1.
(C.26)
De manera similar, se hace el an´
alisis para la frontera exterior, donde (r = r2 );
N uext = −
1
∆T
−
1
r22
T2 − T1
1
1
r2 − r1
(r2 − r1 )
(C.27)
✟)
(−✟
∆T
1
1
r2 − r1
(r2 − r1 )
(C.28)
simplificando la diferencia de temperaturas
1
1
N uext = − ✟ − 2
∆T
r2
✟
se escribe,
N uext =
−1
r2 −
r22
r1
(r2 − r1 )
(C.29)
o bien, haciendo la expansi´
on del denominador queda como:
N uext =
−1
r1 r2 −r22
r1
(r2 − r1 )
por lo tanto el n´
umero de Nusselt en la frontera externa se expresa de la forma:
(C.30)
67
N uext =
−r1
(r2 − r1 )
r1 r2 − r22
(C.31)
si se multiplica esta expresi´
on por la relacion de radios inversa
r2
r1 ,
tambi´en adimensional, el
Nusselt externo queda como:
N uext
−r1
=
(r2 − r1 )
r1 r2 − r22
r2
r1
✘
=
✘✘
r1✘
(r✘
1 − r2 )
✘
✘r✘
r2✘
(r✘
1−
2)
r2
r1
,
(C.32)
que se simplifica en:
N uext = 1.
(C.33)
Bibliograf´ıa
[Alfe et al., 2002] Alfe, D., Gillan, M., and Price, G. (2002). Composition and temperature of
the earth’s core contrained by combining ab initio calculations and seismic data. Earth and
Planetary Science Letters, 195:91–98.
[Ardes et al., 1997] Ardes, M., Busse, F., and Wicht, J. (1997). Thermal convection in rotating
sphericall shells. Physics of the Earth and Planetary Interiors, 99:55–57.
[Arpaci and Larsen, 1984] Arpaci, V. S. and Larsen, P. S. (1984). Convection heat transfer.
Prentice Hall.
[Aubert et al., 2008] Aubert, J., Amit, H., Hulot, G., and Olson, P. (2008). Thermochemical
flows couple the earth’s inner core growth to mantle heterogeneity. Nature, 454:758–761.
[Avila, 2008] Avila, R. (2008). The spectral element method and the meshless local petrov
galerkin method for the solution of the fluid equations.
[Busse et al., 2003] Busse, F., Zaks, M., and Brausch, O. (2003). Centrifugally driven thermal
convection at high prandtl numbers. Physica D, 184:3–20.
[Busse, 1970] Busse, F. H. (1970). Thermal instabilities in rapidly rotating systems. Journal of
Fluid Mechanics, 44:441–460.
[Cabello, 2009] Cabello, A. (2009).
Simulaci´
on num´erica del flujo en el interior de esferas
conc´entricas con rotaci´
on y cambo de fase.
[Chandrasekhar, 1961] Chandrasekhar, S. (1961). Hydrodynamic and hydromagnetic stability.
Oxford University Press.
[Dumas and Leonard, 1994] Dumas, G. and Leonard, A. (1994). A divergence-free spectral expansions method for three-dimensional flows in spherical-gap geometries. Journal of Computational Physics, 111:205–219.
BIBLIOGRAF´
IA
70
[Evonuk and Glatzmaier, 2006] Evonuk, M. and Glatzmaier, G. (2006). A 2d study of the effects
of the size of a solid core on the equatorial flow in giant planets. Icarus, 181:458–464.
[Futterer et al., 2010] Futterer, B., Egbers, C., Dahley, N., Koch, S., and Jehring, L. (2010).
First identification of sub-and supercritical convection patterns from ’geoflow’, the geophysical
flow simulation experiment integrated in fluid science laboratory. Acta Astronautica, 66:193–
200.
[Futterer et al., 2008] Futterer, B., Gellert, M., von Larcher, T., and Egbers, C. (2008). Thermal convection in rotating spherical shells: An experimental and numerical approach within
geoflow. Acta astronautica, 62:300–307.
[Greenspan, 1968] Greenspan, H. P. (1968). The theory of rotating fluids. Cambridge Press.
[Hernlund and Tackley, 2008] Hernlund, J. W. and Tackley, P. J. (2008). Modeling mantle
convection i n the spherical annulus. Physics of the Earth and Planetary Interiors, 171:48–54.
[Hollerbach et al., 2006] Hollerbach, R., Junk, M., and Egbers, C. (2006). Non-axisymmetric
instabilities in basic state spherical couette flow. Fluid Dynamics Research, 38:257–273.
[Jordan, 1979] Jordan, T. H. (1979). Structural geology of the earth’s interior. Proc. Natl. Acad.
Sci., 76(9):4192–4200.
[Karniadakis and Sherwin, 1999] Karniadakis, G. E. and Sherwin, S. J. (1999). Spectral /hp
Element Methods for CFD. Oxford University Press.
[Kelley et al., 2007] Kelley, D. H., Triana, S. A., Zimmerman, D., Tilgner, A., and Lathrop, D.
(2007). Inertial waves driven by differential rotation in a planetary geometry. Geophysical
and Astrophysical Fluid Dynamics, 101(5-6):469–487.
[King et al., 2009] King, E., Stellmach, S., Noir, J., Hansen, U., and Aurnou, J. (2009). Boundary layer control of rotating convection systems. Nature, 457:301–304.
[Lorenzani and Tilgner, 2001] Lorenzani, S. and Tilgner, A. (2001). Fluid instabilities in precessing spheroidal cavities. Journal of Fluid Mechanics, 447:111–128.
[Marcus and Tuckerman, 1987] Marcus, P. and Tuckerman, L. (1987). Simulation of flow between concentric rotating spheres. part 1. steady states. Journal of Fluid Mechanics, 185:1–30.
[Monnereau et al., 2010] Monnereau, M., Calvet, M., Margerin, L., and Souriau, A. (2010). Lopsided growth of earth’s inner core. Science, 328:1014–1017.
BIBLIOGRAF´
IA
[Nair et al., 2005] Nair, R. D., Thomas, S. J., and Loft, R. D. (2005). A discontinous galerkin
transport scheme on the cubed sphere. Monthly Weather Review, 133:814–828.
[Noir et al., 2001] Noir, J., Jault, D., and Cardin, P. (2001). Numerical study of the motions
within a slowly precessing sphere at low ekman number. Journal of Fluid Mechanics, 437:283–
299.
[Olson, 2011] Olson, P. (2011). Laboratory experiments on the dynamics of the core. Physics
of the Earth and Planetary Interiors, 187:1–18.
[Patera, 1983] Patera, A. T. (1983). A spectral element method for fluid dynamics: laminar flow
in a channel expansion. Journal of Computational Physics, 54:468–488.
[Roberts, 1968] Roberts, P. H. (1968). On the thermal instability of a rotating fluid sphere
containing heat sources. Phil. Trans. R. Soc. Lond. Series A, 263:93–117.
[Ronchi et al., 1996] Ronchi, C., Iacono, R., and Paolucci, P. S. (1996). The ”cubed sphere”: A
new method for the solution of partial differential equations in spherical geometry. Journal
of Computational Physics, 124:93–114.
[Simitev and Busse, 2003] Simitev, R. and Busse, F. (2003). Patterns of convection in rotating
spherical shells. New Journal of Physics, 5:97.1–97.20.
[Stixrude et al., 1997] Stixrude, L., Wasserman, E., and Cohen, R. (1997). Composition and
temperature of earth’s inner core. Journal of Geophysical Research, 102:24,729–24,739.
[Tilgner, 1996] Tilgner, A. (1996). High rayleigh-number convection in spherical shells. Physical
Review, 53(5):4847–4851.
[Tilgner, 2007] Tilgner, A. (2007). Kinematic dynamos with precession driven flow in a sphere.
Geophysical and Astrophysical Fluid Dynamics, 101(1):1–9.
[Travnikov et al., 2001] Travnikov, V. V., Rath, H. J., and Egbers, C. (2001). Stability of natural
convection between spherical shells: energy theory. Heat and Mass Transfer, 45:4227–4232.
[White, 2003] White, F. M. (2003). Fluid Mechanics. Mc Graw Hill, 5th Edition.
[Wimmer, 1978] Wimmer, M. (1978). Experiments on a viscous fluid flow beetween concentric
rotating spheres. Journal of Fluid Mechanics, 78:317–335.
[Wu and Roberts, 2009] Wu, C.-C. and Roberts, P. (2009). On a dynamo driven by topographic
precession. Geophysical and Astrophysical Fluid Dynamics, 103(6):467–501.
71
BIBLIOGRAF´
IA
72
[Zhang, 1994] Zhang, K. (1994). On coupling between the poincar´e equation and the heat
equation. Journal of Fluid Mechanics, 268:211–229.
[Zhang et al., 2010] Zhang, K., Chang, K., and Liao, X. (2010). On fluid flows in precessing
spheres in the mantle frame of reference. Physics of Fluids, 22:116604.
[Zhang and Liao, 2004] Zhang, K. and Liao, X. (2004). A new theory for convection in rapidly
rotating spherical systems. XXI ICTAM.