Condiciones de Contorno en Magnetohidrodinámica
Condiciones de Contorno en Magnetohidrodinámica
Tesis Doctoral
Gasdinámicas y Magnetohidrodinámicas.
Corona Solar
10 de Diciembre de 2015
Condiciones de Contorno basadas en Características
Magnetohidrodinámicas.
Solar
por
Director
Comisión Asesora
Dra. Andrea Costa
Asociación
FCEFYN-UNC
FCEFYN-UNC
Córdoba, Argentina
10 de Diciembre de 2015
A mi abuelo Miguel, por inspirarme a ser lo que soy hoy.
Agradecimientos
incondicionales hicieron mucho más fácil mi vida como doctorando. Sobre todo a mis
También me siento muy agradecido por tener muchos muy buenos amigos, que
frustración. Siempre tuvieron buenos consejos para mí, compartimos muy lindos
eterna. A los que están cerca, a los que están lejos , y a los que ya no están.
codirector el Dr. Sergio Elaskar por su conanza y apoyo durante esta tesis, así
esta investigación.
ix
Resumen
Palabras claves: Magnetohidrodinámica - Volúmenes Finitos- Condiciones de Con-
torno - Corona Solar .
Continuando con la línea de estudio de trabajos anteriores, empleamos un esque-
aproximado del tipo de Roe. Para el caso del modelo MHD 1D usamos un solver
las variables y vectores propios del sistema, para evitar la superposición no física de
especíca.
del problema, para luego estudiar la evolución del sistema. Se estudió, por otro lado,
la interacción entre los ujos convectivos, los difusivos y los efectos de los términos
las observaciones.
xi
Abstract
Keywords: Magnetohydrodynamics - Finite Volume Method- Boundary Conditions-
Solar Corona.
Continuing with the previous work of our research group, we used a high reso-
lution TVD scheme to numerically model the dynamics of solar coronal loops. We
used the one dimensional Godunov type ux developed by Harten and Yee, for both
Riemann solver to solve interphase uxes. For the MHD system we used an exact Roe
solver developed by Cargo and Gallice. We also normalised all variables, eigenvalues
We also included source and diusive terms that account for nonlinear heat trans-
fer, radiation emission of energy phenomena specical to coronal loops and energy
dary condition model, which allowed to deal with incoming and outgoig waves at
the system at the boundaries. We also used a TVD backward Euler implicit scheme
robust enough to solve for the complex and demanding source and diusive terms
variables to study the time evolution of the system. We then studied the interaction
between convective and diusive processes, and the inuence of source terms on
them, obtaining dierent wave patterns. Those patterns had velocities and periods
xiii
Zussamenfassung
Schlusselwörter: Magnetohydrodynamik - Finite-Volumen-Verfahren - Grenzbedin-
gungen -Sonnen Corona.
Um an die früheren Arbeiten unserer Forschungsgruppe anzuschlieÿen, haben wir
Dynamik der koronalen Bögen zu modellieren. Wir haben dafür die 1D Godunov-
Typ Flussfunktion von Harten und Yee verwendet, sowohl für die Gasdynamik als
sen, haben wir einen Roe Löser benutzt. Für das MHD-System verwendeten wir einen
exakten Roe Löser, der von Cargo und Gallice entwickelt wurde. Des Weiteren haben
wir alle Variablen, Eigenwerte und Eigenvektoren des Systems normalisiert, um die
Gleichzeitig haben wir auch Quellen und diusive Eekte, nichtlineare Wärmeü-
bertragung und spezische Strahlungsemission für die koronalen Bögen und Energie-
Darüber hinaus haben wir für das beschriebene Schema ein Verfahren zur Her-
leitung der Grenzbedingungen entwickelt, das sich auf die charakteristischen Wellen
Zusätzlich wird eine korrekte Auösung der ein- und auslaufenden Wellen im Ar-
Zeitintegrationsschema genutzt.
Nachdem wir unser Modell kalibriert und getestet hatten, haben wir hydrostatis-
che Lösungen aus der Literatur verwendet, um unter dem Einuss von Störungen die
Dynamik des Systems zu studieren. Wir haben eine starke Interaktion zwischen den
xv
xvi
Símbolos
Escalares
Griegos
γ : Exponente isoentrópico
η : Resistividad eléctrica
κ : Conductividad térmica
λD : Longitud de Debye
xvii
xviii
(l)
λj+1/2 : valor propio asociado al autovector l−ésimo en la interfaz de las celdas j y
j+1
ρ : Densidad de masa
φ : Potencial eléctrico
Latinos
±
Cj+1/2 : Coecientes de expansión del ujo numérico para una ley de conservación
D : Difusividad térmica
EH : Función de calentamiento
FO : Número de Fourier
celda j
kB : Constante de Boltzmann
L : Longitud característica
Lrad : Balance de calor por unidad de volumen cedido al exterior por radiación
n : Densidad de partículas
p : Presión
T : Temperatura
vación escalar
Vectores y Matrices
Griegos
Latinos
A : Matriz jacobiana del vector de ujo hiperbólico respecto a las variables conser-
vativas
et al. (1985)
g : Aceleración de la gravedad
j : Densidad de corriente
conservadas
xxi
n̂ : Versor normal
v : Vector velocidad
Abreviaturas y acrónimos
BC : Boundary Conditions, condiciones de contornos
Agradecimientos viii
Resumen x
Abstract xii
Zussamenfassung xiv
1. Introducción 1
1.1. Plasma en la corona solar . . . . . . . . . . . . . . . . . . . . . . . . 1
3. Modelos Físicos 17
3.1. Términos fuente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
xxiii
xxiv Índice general
4. Esquema Numérico 67
4.1. Discretización en volúmenes nitos . . . . . . . . . . . . . . . . . . . 67
5. Condiciones de Contorno 91
5.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
operadores L . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6. Resultados 125
6.1. Casos gasdinámicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Bibliografía 209
Capítulo 1
Introducción
vista. Esta región se conoce como Fotosfera. Existe otra región por encima de la
supercie solar, donde la densidad del medio es mucho menor, la temperatura mucho
más elevada y donde ocurren una gran variedad de fenómenos dinámicos. Dicha
región se conoce como Corona solar, y está compuesta principalmente por hidrógeno
completamente ionizado, más algunos oligoelementos responsables de la mayoría de la
onda extremo ultravioletas (EUV), rayos X blandos (SXR) y rayos X duros (HXR).
de los 90 varias misiones espaciales (entre ellas las Yohkoh de Japón (Ogawara et al.,
1991), la SoHo y la TRACE (Handy et al., 1999)) permitieron obtener gran cantidad
1
2 Capítulo 1. Introducción
Figura 1.1: Esquema de las distintas regiones presentes en la atmósfera solar (tomado
del sitio ocial de NASA: [Link]
de datos en las logitudes EUV y SXR. En la corona solar también existen fenómenos
corona solar. Generalmente suele distinguirse entre la llamada Región activa, donde
la actividad magnética es mucho más intensa y existen puntos de concentración del
1.1. Plasma en la corona solar 3
campo magnético de diferente polaridad. Por esta razón suelen formarse estructuras
suelen existir regiones donde las líneas de campo son abiertas, produciéndose ujo
como soft X-Ray jets. Los fenómenos anteriormente mencionados ocurren de manera
lenta y no producen grandes desequilibrios energéticos locales. Sin embargo, existen
energía magnética, tales como las fulguraciones (o ares ) o las eyecciones de masa
p
β= B2
(1.1)
2µ0
el ujo del plasma en la dirección de sus líneas, haciendo que éstas actúen como
condiciones, por ejemplo dentro de un arco coronal de gran longitud. Esto suele traer
del arco se tienen regiones con densidades mucho menores, lo que puede inducir que
adecuados.
ticos. Ellos se originan en regiones donde el campo magnético posee gran intensidad,
más recientes conrman que son más bien haces de múltiples estructuras de pequeños
4 Capítulo 1. Introducción
do el campo magnético como barrera. De esta manera, cada arco tiene sus propias
ña atmósfera aislada del resto. La dinámica del plasma en los arcos coronales está
La temperatura en la cromosfera suele ser del orden de los 7000 K, pero una vez
que el plasma uye en el interior de un arco ésta crece hasta valores de 106 K en el
ápice del mismo (Aschwanden, 2004). Este fenómeno es aún discutido en la comuni-
dad cientíca. Si bien los modelos que pretenden explicar el calentamiento son muy
néticas que emergen a la corona. Por otro lado, la radiación electromagnética juega
importancia en las zonas cercanas a las bases (Aschwanden, 2004). Por otro lado, la
arcos debido a que la conductividad térmica es no lineal y muy elevada, haciendo que
en algunos casos este proceso sea más efectivo que la convección para el transporte
de energía
Para modelar de forma correcta a los arcos coronales es necesario tener en cuenta
ralización hecha por Serio et al. (1981). Estos trabajos seminales permitieron obtener
1.1. Plasma en la corona solar 5
soluciones numéricas.
imposibles (como presiones negativas) (Einfeldt et al., 1991). Además, en general los
grandes gradientes espaciales y temporales que ocurren en los arcos coronales suelen
Por otro lado, los fenómenos energéticos no lineales que ocurren agregan com-
plejidad. Uno de ellos es la emisión de energía por radiación. Ésta se produce por
mentos más pesados. Para modelar estos fenómenos es necesario emplear un modelo
multiespecies con una base de datos con las tasas de reacción. Esto vuelve al sistema
les robustos y pasos de tiempo más pequeños. Müller et al. (2003) emplean las rutinas
HAO-DIAPER para obtener las tasas de emisión de radiación de cada una de las
dominios de grandes dimensiones y son necesarias discretizaciones muy nas del mis-
ecuaciones por celda, lo cual insume un gran costo computacional. En los trabajos
de alta performance.
fenómeno astrofísico. Para obtener mayor precisión en el cálculo de los términos con-
vectivos usamos un esquema TVD de alta resolución, que permite resolver disconti-
donde la solución es contínua. Para tratar adecuadamente los términos fuente fue
condiciones iniciales más sosticadas que las empleadas en trabajos anteriores del
más sosticado para tener en cuenta fenómenos transitorios y modelar la física del
los operadores desarrollados por Thompson (1987), varios autores han adaptado y
1.2. Organización de esta tesis 7
como por ejemplo los trabajos de Poinsot y Lele (1992), Sutherland y Kennedy
(2003), [Link] (2004), etc. Estos modelos han sido utilizados en problemas de
ejar o transmitir correctamente distintos tipos de ondas que llegan a los contornos.
condiciones de no reexión de ondas para MHD con términos fuente basado en ex-
en características para resolver un problema similar al del viento solar, pero para
correctamente la dinámica de los términos parabólicos y fuente. Por eso, buena parte
ral y espacial del sistema se debe a familas de ondas no lineales que se propagan a
partir de perturbaciones. Por esta razón una buena comprensión de las propiedades
matemáticas de los sistemas hiperbólicos permite diseñar y entender los métodos nu-
En el Capítulo 3 describimos los modelos físicos con los que usualmente se es-
8 Capítulo 1. Introducción
tudian los arcos coronales: los términos fuente existentes, el modelo gasdinámico de
de los valores y vectores propios del sistema MHD para salvar casos en los que el
modelos físicos: Describimos el esquema TVD de Harten Yee utilizado para las simu-
En el Capítulo 6 discutimos primero los casos que usamos para validar las herra-
micos con distintos tipos de términos fuente, luego soluciones hidrostáticas para arcos
coronales.
la propiedad de que pueden expresarse como leyes de conservación. Por esta razón, se
han desarrollado gran cantidad de esquemas numéricos para resolver problemas rela-
Por esta razón, resulta pertinente hacer un breve resumen de las propiedades de las
leyes de conservación
La función de ujos es una función vectorial, que para el caso cartesiano tridimen-
9
10 Capítulo 2. Sistemas Hiperbólicos y sus Propiedades
sional puede expresarse mediante tres funciones escalaresf (u), g(u), h(u) asociadas
diferencial:
∂u
+ ∇ · F(u) = 0 (2.3)
∂t
Para una variable vectorial, el concepto de ley de conservación puede generali-
∂U
+ ∇ · F (U) = 0 (2.4)
∂t
donde
u1 f1 (u1 , u2 , . . . , um ) g1 (u1 , u2 , . . . , um ) h1 (u1 , u2 , . . . , um )
u2 f2 (u1 , u2 , . . . , um ) g2 (u1 , u2 , . . . , um )
h2 (u1 , u2 , . . . , um )
U = u3 ; F (U) = f3 (u1 , u2 , . . . , um ) g3 (u1 , u2 , . . . , um ) h3 (u1 , u2 , . . . , um )
... ... ... ...
um fm (u1 , u2 , . . . , um ) gm (u1 , u2 , . . . , um ) hm (u1 , u2 , . . . , um )
(2.5)
la ecuación.
∂U
+ ∇ · F (U) = S(U) (2.6)
∂t
Si el vector de términos fuente a su vez es función del gradiente del vector de
mayoría de los problemas de mecánica de uidos y plasmas el ujo está dominado por
∂F (U) ∂F ∂U ∂U
= = A (U) (2.7)
∂x ∂U ∂x ∂x
matriz jacobiana del sistema para la dirección x. De la misma manera, para las
∂U ∂U ∂U ∂U
+ A (U) + Ay (U) + Az (U) =0 (2.8)
∂t ∂x ∂y ∂z
∂U ∂U
En el caso de que el sistema sea unidimensional, las derivadas y son nulas,
∂y ∂z
por lo tanto sólo es necesario calcular la jacobiana A (U)
A continuación se presenta el problema de Riemann para analizar las distintas
x=0
∂U ∂F
+ =0 {−∞ < x < ∞, t > 0} (2.9a)
∂t ∂x
U(i) si x < 0
L
U(t = 0, x) = (2.9b)
U(i) si x ≥ 0
R
rísticas, cuyas velocidades de propagación son los valores propios λi del sistema.
12 Capítulo 2. Sistemas Hiperbólicos y sus Propiedades
Si todos los valores propios del sistema son reales y cada uno posee un vector
de ondas, y los vectores propios el cambio en las variables de estado a través de una
onda característica. Todo sistema hiperbólico posee dos tipos de vectores propios:
matriz L, cuyas las son los vectores propios izquierdos, dichas matrices cumplen la
relación:
R = L−1 (2.12)
Genuinamente no lineal
Un campo característico es genuinamente no lineal si, para alguna normaliza-
ción de R(i) , se cumple que el gradiente del valor propio asociado es monótono
∇U λi · R(i) 6= 0 (2.13)
Linealmente degenerado
2.3. Problema de Riemann 13
∇U λi · R(i) = 0 (2.14)
No convexo
Si el campo característico considerado para algunos valores de U es genuina-
mente no lineal, y en otros es linealmente degenerado, se dice que el campo
característico es no convexo
ondas de la misma familia viajen en la misma dirección, sin alcanzarse nunca. Este
Cuando los vectores propios del sistema son linealmente independientes, pueden
expresarse los estados izquierdo y derecho como combinación lineal de los mismos:
m
X m
X
(i)
UL = αi R ; UR = βi R(i) (2.15)
i=1 i=1
Deniendo la matriz R como la matriz cuyas columnas son los vectores propios del
sistema, puede obtenerse un nuevo conjunto de variables mediante la transformación
W = R−1 U (2.16)
∂wi ∂wi
+ λi =0 (2.18a)
∂t ∂x
α si x < 0
i
wi (t = 0, x) = (2.18b)
β si x ≥ 0
i
m
X I
X
U(x, t) = αi R(i) + βi R(i) (2.20)
i=I+1 i=1
∂U ∂V ∂Ui
=P donde Pij = (2.21)
∂t ∂t ∂Vj
2.4. Formulación en variables primitivas 15
∂U ∂F
P−1 + P−1 =0
∂t ∂x
∂U ∂F ∂V
P−1 + P−1 =0
∂t ∂V ∂x
∂V ∂V
+ Ap (V) =0 (2.22)
∂t ∂x
concordancia con las leyes de conservación físicas del problema. Esto puede verse
Sin embargo, se demuestra que los valores propios del sistema primitivo coinciden
con los del sistema conservativo, y los vectores propios derechos pueden transformarse
donde R(i) son los vectores propios derechos del sistema en variables conservativo, y
(i)
r son los del sistema en variables primitivas. De la misma manera, la transformación
Donde L(i) son los vectores propios izquierdos del sistema en variables conservativo,
(i)
y l son los del sistema en variables primitivas.
Capítulo 3
Modelos Físicos
El modelo Magnetohidrodinámico
Surge de una combinación de las ecuaciones de Euler con las ecuaciones de
que 1 (Aschwanden, 2004). De esta manera, las líneas de campo magnético son lo
uye a lo largo de las líneas de campo. Por otro lado, en la presente tesis se tuvieron
17
18 Capítulo 3. Modelos Físicos
se realiza a través de difusión térmica por conducción de calor a través de sus bases,
y por radiación hacia el exterior a lo largo del arco. Además, existe un incremento
de energía dentro de la corona cuyas causas aún no están del todo claras (Coen,
2008), pero que también cumple un rol crucial en la dinámica energética. Dichos
soluciones analíticas para casos simplicados, que proporcionan leyes de escala y es-
tabilidad para distintos parámetros del problema (Serio et al., 1981). En el presente
trabajo se tuvieron en cuenta dichos efectos como términos fuentes. Se aplican sobre
propuestas en la literatura.
literatura (Aschwanden, 2004). Éste expresa a las pérdidas radiativas como función
emisiones observada . Dicha función fue obtenida por Rosner et al. (1978) como una
a trozos.
ajustamos una relación potencial para 104,0 < T < 104,2 . Las funciones de correlación
para Λ(T ) se presentan a continuación:
−58,70 8,81
10
T ; 104,0 < T < 104,2
10−21,68 ; 104,2 < T < 104,3
10−21,85 ; 104,3 < T < 104,6
10−31 T 2 ;
104,6 < T < 104,9
Λ(T ) = (3.2)
10−21,2 ; 104,9 < T < 105,4
10−10,4 T −2 ; 105,4 < T < 105,75
10−21,94 ; 105,75 < T < 106,3
10−17,73 T −2/3 ;
106,3 < T < 107
20 Capítulo 3. Modelos Físicos
La conducción de calor sigue la Ley de Fourier, que dene al vector ujo de calor
como:
q = −κ∇T (3.3)
variables termodinámicas.
el modelo de Spitzer, basado en el modelo desarrollado por el mismo autor para ob-
tener la conductividad eléctrica (Spitzer, 1962). Dicho modelo propone una relación
no lineal para calcular el ujo de calor en la dirección paralela a las líneas del campo
calor debido a la conducción se debe a las colisiones entre electrones y iones. Luego,
5/2
1,84 · 10−5 Te
κ= (3.4)
lnΛ
−1/2 Te
lnΛ = 29,7 + ln n (3.5)
106
1,84 · 10−5
κsp = ≈ constante (3.6)
lnΛ
2 d2 T 7/2
d 5/2 dT
Q̇cond = −∇ · q = κsp T = κsp (3.8)
dx dx 7 dx2
con la hipótesis necesaria. Sin embargo, existe evidencia de que la ley de Spitzer
Poedts, 2004).
calor por conducción está saturado. Los autores Cowie y McKee (1977) describen el
modelo para la saturación del ujo de calor. El valor máximo del ujo de calor por
3
qmax = ne kB Te vc (3.9)
2
corriente neta. Los autores demostraron que si la fuente de calor posee una función
1/2
2kB Te
qsat = 0,4 ne kB Te (3.10)
πme
del sonido isotérmica como ciso = (ne kB Te )1/2 e introduciendo un factor φ para tener
22 Capítulo 3. Modelos Físicos
en cuenta las incertidumbres, puede expresarse el ujo de calor por conducción como:
conducción de la forma
qsat
qcond lim = qcond (3.12)
qsat + kqcond k
Normalizando los ujo de calor por conducción respecto al ujo saturado según
obtenemos la expresión
1
qcond lim = q (3.14)
1 + kqcond k cond
En la Fig. 3.2 comparamos ambos ujos de calor normalizados
0.8
0.6
0.4
0.2
q cond lim
−0.2
q cond lim
−0.4 q cond
−0.6
−0.8
−1
−5 −4 −3 −2 −1 0 1 2 3 4 5
q cond
Figura 3.2: Comparación entre los ujos de calor por conducción limitado y conven-
cional
3.2. Ecuaciones de Euler 23
T < 104 K) hacia la Corona (donde T ≈ 106 K) a través de una pequeña zona
por el cual esta transformación se lleva a cabo no se conoce con certeza todavía.
varios autores modelan dicho efecto de calentamiento mediante funciones simples con
un solo máximo en toda la longitud del arco. En el presente trabajo se empleó una
(t − t0 )2
x − x0
EH (x, t) = EH 0 exp exp 2
para0 < x < L/2 (3.15)
Hcal τcal
permitir concentrar el efecto de calentamiento hacia las bases o el ápice. Este tipo
en varios trabajos, como (Aschwanden y Tsiklauri, 2009; Rosner et al., 1978; Serio
v(t, x) := [u, v, w]
masa expresa que la masa de una partícula material arrastrada por el ujo no cambia
∂ρ
+ ∇ · ρv = 0 (3.16)
∂t
ujo es igual a la sumatoria de fuerzas externas que actúan sobre dicha partícula
material.
∂ρvT T
+ ∇ · ρvT ⊗ v + Ip = Scm (U) (3.17)
∂t
donde p(t, x) es la presión del uido. Y donde Scm (U ) son términos fuente asociados
propio del uido como términos fuente, que se expresan en función de la aceleración
diferencial como:
∂ρ ei + 12 v · v
1
+∇· ρei + ρv · v + p v = ∇ · q + Se (U) (3.19)
∂t 2
3.2. Ecuaciones de Euler 25
entalpia:
p
H = ei + (3.20)
ρ
y la entalpia total
1 2
Ht = H + |v| (3.21)
2
Se demuestra que para un ujo adiabático la entalpia total se conserva a lo largo de
q = −κ∇(T ) (3.22)
hacia el exterior por emisión de radiación, Ec.(3.1). Además debe incluirse el efecto
Las tres Leyes de conservación (3.16), (3.17) y (3.19) forman el sistema de ecua-
1
et = ei + v2 (3.24)
2
26 Capítulo 3. Modelos Físicos
∂U
∂t
+ (∇ · F) = S(U) en x ∈ Ω, t > 0
(3.25a)
U(t = 0, x) = φ(x)
donde
ρ ρv
U = ρvT , F = ρvT ⊗ v + Ip (3.25b)
ρet (ρet + p) v
Una ecuación de estado (EOS) dada por relaciones termodinámicas provee la ecua-
ción adicional necesaria para la clausura del sistema. Si se considera en todos los
casos que el uido es una mezcla de gases térmicamente perfectos, puede emplearse
p
ei = (3.26)
ρ(γ − 1)
Cp
donde γ= es el exponente isoentrópico para el gas o mezcla de gases considerados,
Cv
y Cp y Cv son los calores especícos a presión y volumen constante, respectivamente.
ei
T = (3.27)
Cv
R
Cv = (3.28)
γ−1
3.2. Ecuaciones de Euler 27
0 1 0 0 0
−u + γ−1
2
2 v
2
(3 − γ) u (1 − γ) v (1 − γ) w (γ − 1)
Ac =
−uv v u 0 0
(3.29)
−uw w 0 u 0
γ−1
−γ etρu + (γ − 1) uv2 γ eρt
2 2
− 2 v + 2u (γ − 1)uv (γ − 1)uw γu
1 1 0 0 1
u
u−a
0
0
u+a
R1c = v ; R2c = v ; R3c = 1 ; R4c = 0 ; R5c = v
w w 0 1 w
1 2
Ht − ua v v w Ht + ua
2
(3.30)
donde a es la velocidad del sonido y Ht la entalpía total del ujo, Ec. (3.21).
1 0 0 0 0
u ρ 0 0 0
P= v 0 ρ 0 0 (3.31)
w 0 0 ρ 0
1 2 1
2 (u + v2 + w2 ) ρu ρv ρw γ−1
u ρ 0 0 0
u2 2ρu 0 0 1
Q = uv ρv ρu 0 0 (3.32)
uw ρw 0 ρu 0
1 2 1 2 γ γ
2 uv 2 ρv + ρu2 + γ−1 p ρuv ρuw γ−1 u
28 Capítulo 3. Modelos Físicos
sección varía, pueden extenderse las ecuaciones de Euler para tener en cuenta dichos efectos.
expresarse como:
∂ρA ∂ρAu
+ =0 (3.33)
∂t ∂x
Se observa que en este caso la magnitud transportada por el uido es un producto entre
∂ρAu ∂ ∂A
+ (ρAu2 + pA) − p =0 (3.34)
∂t ∂x ∂x
Y la conservación de la energía
∂ρet A ∂
+ [u(ρet A + pA)] = 0 (3.35)
∂t ∂x
ρA ρAu 0
∂ ∂
ρAu + ρAu2 + pA = p ∂A (3.36)
∂t ∂x ∂x
ρet A (ρet A + pA)u 0
Se observa que las variables conservativas y los ujos son los de las ecuaciones de Euler
la del sistema de Euler con área constante, por ende es independiente del área. Lo mismo
siguiente manera:
ρA
u ρA 0
ρA
0
∂ 1 ∂
p ∂A
u + 0 u u = (3.37)
∂x
∂t ρA ∂x
ρA
pA 0 γpA u pA ∂A
(γ − 1)pu ∂x
3.2. Ecuaciones de Euler 29
1 0 0 u ρA 0
P= u ρA 0 ; Q = u2 2ρAu 1 (3.38)
1 2 1 1 3 3 2 γ γ
2u ρAu γ−1 2u 2 ρAu + γ−1 pA γ−1 u
Como en la presente tesis se trataron efectos de ujos difusivos y términos fuente, resul-
a que las ecuaciones que componen el sistema de leyes de conservación poseen diferentes
unidades físicas, ocurre que las diferentes componentes de los vectores propios pueden te-
ner órdenes de magntud muy distintos. Esto resulta problemático desde el punto de vista
fuente y difusivos quedan afectados por números adimensionales que dan información extra
sobrelineado en esta sección para identicar a las variables adimensionales, y el subjo ref
para identicar los valores de referencia.
A
=A
Aref
A continuación, se obtienen los valores por los que tienen que reemplazarse las variables
Longitudes
neral suele ser igual a la longitud o semilongitud del dominio considerado. Para cada una
30 Capítulo 3. Modelos Físicos
xi
xi = → xi = Lxi (3.39)
L
expresarse como:
∂ 1
∇= = L∇ → ∇ = ∇ (3.40)
∂xi L
Densidad
estar relacionada a alguna condición física, como el valor en el contorno o un valor mínimo
esperable.
ρ
ρ= → ρ = ρρref (3.41)
ρref
Debe tenerse en cuenta que, como los datos experimentales y valores de referencia están
dados en función de la densidad de iones ni , éstos deben multiplicarse por la masa promedio
de todas las especies presentes para transformar a la variable ρ.
Velocidades
Optamos por adimensionalizar las componentes del vector velocidad vi mediante una
velocidad de sonido a,
γp
a2 = (3.42)
ρ
ésta es la velocidad a la que se propaga la información física en el presente modelo. Como
en el contorno
vi
vi = → vi = aref vi (3.43)
aref
donde
γpref
a2ref = (3.44)
ρref
3.2. Ecuaciones de Euler 31
Tiempo
aref L
t=t →t=t (3.45)
L aref
∂ ∂ L ∂ ∂ aref ∂
= aref = → = (3.46)
∂t L ∂t
aref ∂t ∂t L ∂t
L
tref =
aref
Presión y energía
entalpia total. Analizando las dimensiones de esta ecuación observamos que la presión y
y al de la energía como:
entonces:
p
p= → p = ρref a2ref p (3.49)
ρref a2ref
ei
ei = → ei = a2ref ei (3.50)
a2ref
32 Capítulo 3. Modelos Físicos
Temperatura
L2
[T ] = 2
t
1
T =T → T = T a2ref (3.51)
a2ref
κref a2ref
q=− k · ∇T (3.52)
L
Aceleración de la gravedad
para la física del problema, que de acuerdo a la literatura (Aschwanden y Schrijver, 2002;
p0 z
p = p0 exp −
ρ0 g
p0
Λ0 =
ρ0 g
donde ρ0 y p0 corresponden a los valores de dichas variables para una altura de referencia,
ρ0 L
g=g
p0
3.2. Ecuaciones de Euler 33
poder adimensionalizar con la velocidad del sonido. Esto trae como ventaja una ecuación
más simple en variables adimensionales, que a su vez sigue respetando las magnitudes físicas
L
g=g (3.54)
a2ref
sionales, Ecs. (3.39), (3.40), (3.41), (3.51), (3.49), (3.45), y sacando factor común en cada
caso, se obtienen las expresiones adimensionales de cada una de las ecuaciones de conser-
vación.
Para simplicar las expresiones se expresan dichas ecuaciones en notación indicial (por
Continuidad
aref
Sacando factor común
L ρref en la Ecuación (3.16) de conservación de masa, permite
expresarla de manera adimensional como:
aref ∂ρ
ρref + ∇ · (ρvi ) = 0 (3.55)
L ∂t
dimensional
∂ρ
+ ∇ · (ρvi ) = 0 (3.56)
∂t
ρref a2ref
Sacando factor común en el miembro izquierdo de la ecuación de cantidad de
L
moviemento y resolviendo se obtiene la expresión adimensional de esta ley de conservación.
Para el caso ideal, al ser un sistema homogéneo la ecuación adimensional tiene la misma
forma que la dimensional. Sin embargo, en el sistema no ideal (no homogéneo) los términos
34 Capítulo 3. Modelos Físicos
∂ρvi
+ ∇ · [ρvi vj − δij (p)] = ρgi (3.58)
∂t
Conservación de energía
ρref a3ref
∂ρ et
+ ∇ · [(ρ et + p) vi ]
L ∂t
aref 2
κref 2
= 2
aref ∇ · qi + EH + Lrad + ρref aref ρgi · vi
L L
uL
Pe = (3.59)
D
κ
donde D= Cp ρ es la difusividad térmica.
∂ρ et 1 L L
+ ∇ · [(ρ et + p) vi ] = ∇ · qi + 3 EH + Lrad + ρgi · vi (3.60)
∂t Pe ρref aref ρref a3ref
Reemplazando los resultados anteriores (Ecs. (3.56), (3.58) y (3.60)) en las Ecuaciones
de conservación (3.25):
ρ ρv
∂ 1
ρv + ∇ · ρv ⊗ v + I (p) =
∂t L
ρ et (ρ et + p) v
0 0
ρgi
∇· 0 + (3.61)
L L
1 E + L + ρgi · vi
Pe∇ · (κ∇T ) ρref a3ref H ρref a3ref rad
3.3. Modelo Magnetohidrodinámico 35
plasma. Podemos identicar como plasmas a la corona solar, el medio interestelar, el medio
interplanetario, el interior de muchas estrellas, el medio con el que se conna las reacciones
nucleares en un reactor.
mientos colectivos, gracias a la gran interacción que ocurre entre ellas debida a las fuerzas
debido a campos magnéticos externos. Las partículas son libres en tanto la energía cinética
es mayor que la energía de interacción entre ellas, aunque las fuerzas que las
vinculan son de largo alcance. Como están cargadas, apartamientos de carga apreciables
que haya grandes apartamientos de cargas se traduce, en general, en que los plasmas son
es aproximadamente el mismo para longitudes mayores que una longitud característica del
Si se entrega energía a un gas mediante una descarga eléctrica o una elevada temperatura
éste se ioniza produciendo a su vez cargas que son capaces de ionizar. La cantidad relativa de
puede ser obtenida mediante la ecuación de Saha (Goedbloed y Poedts, 2004), asumiendo
36 Capítulo 3. Modelos Físicos
3/2
T 3/2
ni 2πme kB Ui
= exp − (3.62)
ne h2 ni kB T
del gas considerado. Teniendo en cuenta que la energía de ionización del hidrógeno es de
13,6 eV, el exponente de la expresión es −15,8 · 10−4 /T . Esto signica que para que un
gas se encuentre en estado de plasma se requieren temperaturas elevadas, que son poco
comunes en las situaciones de la vida cotidiana. Por esta razón, el estudio de los plasmas
En general en los plasmas suelen ocurrir procesos de intercambio de energía que invo-
lucran radiación, recombinación e intercambio de calor entre especies. Los electrones libres
emiten radiación electromagnética al ser acelerados y los iones, átomos y moléculas emiten
vez produce ionización. Eventualmente se llega a un equilibrio entre los diferentes proce-
sos dadas ciertas condiciones de interacción con las condiciones externas. Los procesos de
ionización más importantes (que involucran átomos neutros A, iones i y electrones e) son,
mico estos procesos se equilibran, cada proceso y su inverso producen la misma cantidad
de reacciones por unidad de tiempo, lo que produce una población relativa de especies es-
tacionaria que es función de la temperatura de cada especie. Sin embargo, puede asociarse
una única temperatura a todos los procesos. Muy rara vez el plasma es lo sucientemente
extenso y denso como para retener la radiación y establecer un equilibrio con ella (es el
caso de los interiores estelares). Sin embargo, en la mayoría de los plasmas es posible tener
3.3. Modelo Magnetohidrodinámico 37
Saha (Ec. (3.62)) permite obtener la población relativa de especies. Sin embargo, en las
condiciones de la corona solar (y también en las condiciones habituales del plasma de labo-
como un continuo debe cumplirse la condición de que el camino libre medio entre colisio-
nes resulte muy pequeño frente a las dimensiones lineales características. Muchos de los
plasmas de interés son tales que las colisiones binarias (de a dos partículas) ocurren poco
frecuentemente; se dice que dichos plasmas son no colisionales (Inan y Golkowski, 2011). Sin
embargo, aun así puede utilizarse la hipótesis del continuo, ya que existen otros mecanismos
físicos distintos a las colisiones binarias que hacen que los plasmas posean comportamiento
Los campos eléctrico E y megnético B son los que se encargan de uniformizar las
propiedades del plasma para longitudes mayores que una cierta longitud característica,
permitiendo tratar al plasma como medio continuo. En ese caso se introduce una magnitud
más apropiada que el camino libre medio que permite distinguir entre el comportamiento
colectivo que genera dicha cuasi neutralidad y el comportamiento en el que dicha neutralidad
podemos asumir que el mismo está compuesto por iones de Hidrógeno y electrones, y que la
el potecial eléctrico φ
qi φ(x)
ni (x) = no i exp − (3.63)
kB T
es posible determinar la distribución espacial de la densidad de carga ρq i = qi n(x) y del
∂ 2 φ(x) qi ni (x)
− 2
= (3.64)
∂x 0
38 Capítulo 3. Modelos Físicos
bución de cargas (por ejemplo que en una pequeña región existen sólo electrones) es posible
encontrar una solución linealizada para el campo eléctrico, empleando las Ecs. (3.63) y
(3.64).
kxk
φ(x) = K exp − (3.65)
λD
Donde λD es la longitud de Debye y K una constante. Dicha longitud cuantica el alcance
longitud de Debye λD pueden ignorarse los efectos de dichas perturbaciones. En este caso,
el plasma puede considerarse como cuasi neutro, y por ende, como un medio contínuo y
homogéneo.
Ello se debe a que las interacciones dominantes son de largo alcance y que pequeños
campos eléctricos locales, generando una dinámica que afecta colectivamente a gran número
de ellas. Es decir, para problemas en los que interesa describir fenómenos con longitudes
características mayores a la longitud de Debye de un plasma, éste puede ser descrito como
un medio continuo en el que prevalece la cuasi neutralidad y para fenómenos cuya longitudes
son menores a la longitud de Debye del plasma, este debe ser descrito como un conjunto
oscilaciones de cargas. Existe una frecuencia asociada a cada especie que lo componen . En
caso de un plasma de hidrógeno existirá una asociada a los iones y otra a los electrones:
s
2
ne,i qe,i
ωp e,i = , (3.67)
0 me,i
los plasmas.
3.3. Modelo Magnetohidrodinámico 39
La fuerza restitutiva lleva a oscilaciones espontáneas del plasma que ocurren a la frecuencia
de plasma. Dado que los electrones se mueven mucho más rápidamente que los iones, éstos
son los responsables fundamentales de la oscilación. En tal caso se puede considerar el mo-
delo sencillo en que los iones están quietos y los electrones no tienen movimiento térmico.
partícula con carga q poseerá un movimiento circular alrededor de las líneas de campo con
cular v⊥
v⊥
rc = (3.69)
Ωc
Las Ecs. (3.68) y (3.69) imponen restricciones sobre las escalas temporales y espaciales
del problema a estudiar: si la escala espacial del problema es mucho mayor al radio de
Larmor o si la escala temporal es mucho mayor que la inversa de frecuencia del ciclotrón,
Para entender al plasma como sistema dinámico hace falta notar que las partículas
cargadas son aceleradas por fuerzas de Lorenz y gradientes de presión, pero a su vez el
partículas cargadas. Existen distintos niveles de aproximación para modelizar las ecuaciones
probabilidad fi (x, w, t), que representa la densidad de probabilidad para que existan par-
tículas con una velocidad w determinada en una posición x y un instante de tiempo dados.
Planteando la conservación de la derivada total de dicha función para cada especie, se
∂fi wi ∂f
i
+ wi · ∇fi + qi E + ×B · =0 (3.70)
∂t c ∂wi
conocidas como ecuaciones de Vlasov, que combinadas con las ecuaciones de Maxwell per-
sus productos con las velocidades estadísticas w, asumiendo algún tipo de distribución para
así obtener ecuaciones de conservación de masa, de cantidad de movimiento y de energía.
Z ∞ Z ∞
ni = fi dw; n i vi = wfi dw; (3.71)
0 0
donde vi es la velocidad local promedio para la especie considerada. Por otro lado, expre-
sando la velocidad probabilística w como una suma de una componente media v más una
0
componente aleatoria w con media igual a 0.
Z ∞ Z ∞
vi · vi + w0 · w0 + 2v · w0 fi dw
mi w · wfi dw = mi
0
Z0 ∞
vi · vi + w0 · w0 fi dw
= mi
0
Z ∞
w0 · w0 fi dw
= mi ni vi · vi + mi (3.73)
0
El término mi ni vi · vi representa la energía cinética del centro de masa del elemento con-
R∞ 0 0
y Pi = m i
siderado,
0 (w · w ) fi dw da origen a un tensor de presiones. Se obtienen así
expresiones para la conservación de la masa y cantidad de movimiento para cada una de
3.3. Modelo Magnetohidrodinámico 41
las i especies.
∂ni
+ ∇ · n i vi = 0 (3.74)
∂t
∂vi
mi ni + (vi · ∇)vi = −∇ · Pi + ni qi (E + vi × B) + Pij (3.75)
∂t
colisiones entre las especies i y j. Nuevamente, combinadas con las ecuaciones de Maxwell,
y considerando los iones y los electrones como dos especies independientes, se obtiene un
set de ecuaciones que se conoce como teoría de dos uidos. Dichas ecuaciones representan
forma relativa a la velocidad media de la especie considerada. Para darle clausura al sistema
es necesario denir una ecuación de energía y una ecuación de estado, o simplemente alguna
al plasma como un solo uido, haciendo un promedio ponderado de todas las especies según
X
ρ= m i ni , (3.76)
X
ρv = m i n i vi , (3.77)
X
j= n i qi v i , (3.78)
donde ρ es la densidad en masa, v es la velocidad del centro de masa del elemento conside-
permite reducir el número de ecuaciones a resolver e ignorar los efectos de las oscilaciones de
alta frecuencia. En el caso de un plasma, en general pueden existir partículas con diferente
carga eléctrica neta y la interacción entre partículas de distinta carga. Dicho modelo surge
trabajo algebraico, permite obtener un nuevo sistema de ocho ecuaciones escalares. Como
las relaciones entre el campo eléctrico y el campo magnético para un medio o sistema de
partícula que compone el sistema, permiten determinar el campo eléctrico E(x, t), el campo
magnético B(x, t) y las fuerzas electromagnéticas que actúan en cada partícula ([Link],
2006).
Ley de Gauss
ρe
∇·E= (3.79)
0
∇·B=0 (3.80)
∂B
∇×E= (3.81)
∂t
Ley de Ampere
∂E
∇ × B = µ0 j + µ0 0 (3.82)
∂t
Teniendo en cuenta la relación
1
µ0 0 =
c2
donde c es la velocidad de la luz, es posible armar que si las velocidades características
cargas( por ejemplo: iones de diversos tipos, electrones libres, etc). Por esta razón, debe
expresarse una ecuación de conservación de masa para cada uno de los conjuntos de partí-
culas existentes. Para el caso más simple de un plasma compuesto sólo por iones con carga
3.3. Modelo Magnetohidrodinámico 43
∆N
nα = lı́m (3.83)
∆V →0 ∆V
∂ni
+ ∇ · (ni vi ) = 0 (3.84)
∂t
∂ne
+ ∇ · (ne ve ) = 0 (3.85)
∂t
la cantidad total de iones, puede denirse una densidad de partículas total de la siguente
manera
n = ne + ni (3.86)
Cuando se considera un plasma globalmente neutro, en virtud de que las escalas del pro-
blema considerado son mucho mayores que la longitud de Debye (Ec. (3.66)), se cumple
ne ≈ ni → n ≈ 2ne
por su masa molecular promedio. Por otro lado, se dene la densidad del plasma como
ρ = mn (3.87)
Donde m es una masa molecular promedio de todas las especies. Para el caso de un plasma
compuesto sólo por iones y electrones:
ρ = ne me + ni mi ≈ ni mi (3.88)
Ya que la masa de los electrones resulta despreciable respecto de la de los iones. Luego,
multiplicando la ecuación de conservación de los iones por una masa promedio de las partí-
44 Capítulo 3. Modelos Físicos
culas, se obtiene una ecuación de continuidad análoga a la del caso hidrodinámico (Euler),
Ec. (3.16).
Euler, con la diferencia que en el miembro derecho aparece una fuerza electromagnética de
∂ρvT T
+ ∇ · ρvT ⊗ v + Ip = j × B + ρg (3.89)
∂t
donde j es la densidad de corriente, dada por la Ec. (3.78). Para un plasma de Hidrógeno
j = qe (ni vi − ne ve ) (3.90)
µ0 j × B = (∇ × B) × B
Usando la identidad
1
(∇ × B) × B = ∇ · (B ⊗ B) − ∇(B · B)
2
se obtiene
1 1
j×B= ∇ · (B ⊗ B) − ∇(B · B)
µ0 2
∂ 1 1
(ρv) + ∇ · (ρv ⊗ v − B ⊗ B) + ∇(p + B · B) = ρg (3.91)
∂t µ0 2µ0
3.3. Modelo Magnetohidrodinámico 45
Para expresar la ley de Inducción Magnética (Ec. (3.81)) en función sólo del campo
por qe /me y qi /mi respectivamente. Teniendo en cuenta las Ecs. (3.76), y teniendo en cuenta
que me mi puede obtenerse:
∂j qe ne qe2 ne qe2
=− ∇ · Pe + (E + v × B) + qe me j × B − η·j (3.92)
∂t me me me
E+v×B=η·j (3.93)
donde η representa un tensor de resistividad general, que representa los efectos anisotrópicos
en la conductividad eléctrica.
permite obtener:
∂B 1
+ ∇ · (v ⊗ B − B ⊗ v) = −∇ × η (∇ × B) (3.95)
∂t µ0
De manera análoga al caso del sistema de Euler, esta relación representa el primer prin-
cinética.
∂ 1 2 1 2
ρv + ρei +∇· v ρv + ρei + ∇ · (pv) = ∇ · (k · ∇T ) + E · j (3.96)
∂t 2 2
46 Capítulo 3. Modelos Físicos
v2 B2
et = ei + + (3.97)
2 2ρµ0
Luego, para introducir los términos magnéticos se suma a ambos miembros la ecuación
B2
ρ∂et B
+ ∇ · v ρet + p + − (v · B) =
∂t 2µ0 µ0
1
∇ · k∇T − [η (∇ × B)] × B + ρvg (3.98)
µ0
Combinando las ecuaciones (3.16), (3.91) y (3.98), puede escribirse el sistema MHD en
ρv
ρ
ρv ⊗ v − 1 B ⊗ B + I p + B 2
∂ ρv
+∇· µ0 2µ0
=
∂t
B
v ⊗ B − B ⊗ v
2
ρet B 1
ρet + p + 2µ 0
v − µ0 (B · v) B
0 0
0 ρg
∇· +
1
(3.99)
0 −∇ × µ0 η (∇ × B)
k · ∇T LRad + EH + ρvg
A su vez deben agregarse dos ecuaciones más para la clausura del sistema:
∂
(∇ · B) = 0 (3.100)
∂t
el tiempo.
no los elimina explícitamente. Para salvar esta condición existen diferentes enfoques
(staggered)
se simplica, ya que al sólo existir derivadas no nulas en una dirección , en este caso
la dirección x:
sionaliza el sistema de Euler, y además porque las expresiones de las matrices jacobianas
y los vectores propios son más sencillas en variables adimensionales, ya que desaparece el
proceso.
Longitudes
manera que para el sistema de Euler, de acuerdo a las Ecs. (3.39), (3.40).
Densidad
los datos y valores de referencia están dados en función de la densidad de partículas n, éstos
deben multiplicarse por la masa promedio de las especies consideradas para transformar a
la variable ρ.
Velocidades
importante (en virtud de que el parámetro β < 1); y a que las velocidades características
del sistema dependen fuertemente del mismo, las velocidades se adimensionalizan mediante
vi
vi = → vi = bxref vi (3.102)
bxref
como:
Bref
bxref = √ (3.103)
µ0 ρref
3.4. Adimensionalización de las ecuaciones 49
µ0 posee otro valor. Existen tres tipos de sistemas CGS: electrostático, electromagnético, y
gaussiano. La física deja claro que existen sólo dos constantes independientes de entre 0 ,
µ0 y c, dependiendo de cómo se dena el sistema de unidades empleado. Las ecuaciones de
∇ · E = 4πkc ρq (3.104)
∇·B=0 (3.105)
∂B
∇ × E = −αl (3.106)
∂t
αb ∂E
∇ × B = 4παb j + (3.107)
kc ∂t
En la tabla 3.1 se expresan los distintos factores de conversión para pasar de un sistema
kc kc
sistema kc αb 0 µ0 ka = c2
αl = αB c2
cgs electrostático(ESU) 1 c−2 1 c−2 c−2 1
2 −2
cgs electromagnético(EMU) c 1 c 1 1 1
CGS Gaussiano 1 c−1 1 1 c−2 c−1
CGS de Lorenz Heavyside
1
4π
1
4πc
1 1
1
4πc2
c−1
c2 1 b 4π 1
SI 1
b b 4πc2 b b
Tabla 3.1: Conversión de sistemas de unidades para las ecuaciones del electromag-
netismo
Campos magnéticos
Para ujos unidimensionales, se cumple que Bx = cte, por lo que el valor de Bx resulta una
elección natural.
Bi
Bi = → Bi = Bi Bref (3.108)
Bref
50 Capítulo 3. Modelos Físicos
Tiempo
bxref L
t=t →t=t
L bxref
∂ ∂ L ∂ ∂ bx ∂
= bxref
= → = ref
∂t ∂t bxref ∂t ∂t L ∂t
L
L
tref =
bxref
Presión
2
Bref
p
p= 2 →p= p
Bref µo
µo
Temperatura
1
T =T → T = T b2xref (3.109)
b2xref
Energía
Por último, la energía total (compuesta por la suma de las energías interna, cinética y
además una densidad de energía (energía por unidad de masa) asociada al campo magnético.
et
et = 2
Bref
→ et = b2xref et
ρµ0
3.4. Adimensionalización de las ecuaciones 51
Continuidad
bxref
Sacando factor común
L ρref en la ecuación de conservación de masa expresada en
bxref
∂ρ
ρref + ∇ · (ρvi ) = 0
L ∂t
dimensional
De la misma manera que para el sistema de Euler, para el caso ideal la ecuación adi-
mensional tiene la misma forma que la dimensional. Pero el miembro derecho del sistema
ρref b2x
ref
Dividiendo a ambos miembros por , y multiplicando y dividiendo el segundo miem-
L
bro por u, se obtiene que los términos no homogéneos quedan multiplicados por:
∂ρvi Bj Bj
+ ∇ · ρvi vj − Bi Bj + δij (p + ) = ρgi
∂t 2
Inducción electromagnética
Reynolds Magnético :
µ0 uL
ReM = (3.111)
η
que representa una medida de la intensidad del acoplamiento entre el ujo y el campo
bxref
Al = (3.112)
u
52 Capítulo 3. Modelos Físicos
obtenemos
∂Bi 1
+ ∇ · (vi Bj − Bi vj ) = − ∇ × η ∇ × Bi (3.114)
∂t RM Al
Conservación de energía
2
bxref Bref
Para la ecuación de la energía se saca factor común
a µ0 del miembro izquierdo,
obteniéndose:
2
bxref Bref
∂et Bj Bj
+∇· et + p + vi − Bj · vj Bi =
L µ0 ∂t 2
3
bxref
kref 2
bx ∇ · q i + EH + L rad + ρref ρgi · vi (3.115)
L2 ref L
Los términos fuente resultan idénticos a los de la ecuación de energía del modelo de Eu-
ler, con la diferencia de que están afectados por un factor de adimensionalización distinto
Teniendo en cuenta las Deniciones del número de Péclet (3.59) y el de Alfvén (3.112)
obtenemos:
∂et Bj Bj
+∇· et + p + vi − Bj · vj Bi
∂t 2
1 L L
= ∇ · qi + 3
EH + Lrad + ρgi · vi (3.116)
P e Al ρref bxref ρref b3xref
3.4. Adimensionalización de las ecuaciones 53
ρv
ρ
2
B
ρv ⊗ v − B ⊗ B + I (p + )
∂ ρv
2
+ ∇ · =
∂t
B
v ⊗ B − B ⊗ v
2
et et + p + B2 v − B · v B
0 0
ρg
0
∇· + (3.117)
1
−∇ × µ0 η ∇ × B
0
L
1
∇·q ρ b3
EH + ρ Lb3 Lrad + ρg · v
P e Al ref xref ref xref
los sistemas MHD adimensionales de siete y ocho ondas respectivamente. Para simplicar la
sional
ρu ρv ρw
ρ ρu2 − B 2 + p + B
2
x ρuv − Bx By ρuw − Bx Bz
ρu
2
B2
ρv
ρuv − Bx By ρvv − By By + p + ρvw − By Bz
2
∂ ρw
+∇· B2
ρuw − Bx Bz ρvw − By Bz ρww − Bz Bz + p + =0
∂t
Bx
2
0 vBx − By u wBx − Bz u
By
uB − Bx v 0 wBy − Bz v
y
Bz
uBz − Bx w vBz − By w 0
et
C1 u − C2 Bx C1 v − C2 By C 1 w − C 2 Bz
(3.118)
54 Capítulo 3. Modelos Físicos
donde:
B2
C1 = e t + p + ; C2 = Bx u + By v + Bz w; B 2 = Bx2 + By2 + Bz2
2
La naturaleza de la ecuación de la inducción (Ec. (3.95)) implica que el vector de ujo tenga
porque aparece un valor propio nulo en el sistema de vectores y valores propios. Los autores
Peyrard y Villedieu (1999) y Sokolov et al. (2008) proponen posibles soluciones para este
inconveniente.
x se expresa como:
0 1 0 0 0 0 0 0
−2u2 + 1−γ 2
2 v (3 − γ) u (1 − γ) v (1 − γ) w 0 (2 − γ) By (2 − γ) Bz (γ − 1)
−uv v u 0 0 −Bx 0 0
−uw w 0 u 0 0 −Bx 0
Ac =
0 0 0 0 0 0 0 0
uBy By
− ρ + vBρ x − Bρx u
ρ 0 0 0 0
− uBz + wBx Bz
0 − Bρx 0 0 u 0
ρ ρ ρ
Ac8,1 Ac8,2 Ac8,3 Ac8,4 0 Ac8,5 Ac8,6 γu
(3.119)
donde
et u 2 − γ B 2 Bx
Ac8,1 = −γ + (γ − 1) uv2 − u + (B · v)
ρ 2 ρ ρ
2 − γ B2 B2
et γ − 1 2
Ac8,2 = γ − 2
v + 2u + − x
ρ 2 2 ρ ρ
Bx By
Ac8,3 = (1 − γ) uv −
ρ
Bz Bx
Ac8,4 = (1 − γ) wu −
ρ
Ac8,6 = (2 − γ) By u − Bx v
Ac8,7 = (2 − γ) Bz u − wBx
La matriz jacobiana para el caso de dos dimensiones o más tendrá ocho autovalores en
donde
2
s 2
2 2
1 γp B γp B γpB
c2f,s = + ± + − 4 2 x (3.121)
2 ρ ρ ρ ρ ρ
λa = u ± bx
donde
Bx
bx = √
ρ
λs = u
λB = 0
El sistema conservativo, así como está planteado, no posee solución para un problema de
Riemann con cualquier tipo de condición inicial, ya que para la matriz jacobiana uno de
los autovalores es 0 y uno de los vectores propios queda indeterminado. Por esta razón es
ρ u ρ 0 0 0 0 0 0 ρ
By
u 0 u 0 0 − Bρx Bz 1
u
ρ ρ ρ
By
v
0 0 u 0 − ρ − Bρx 0 0
v
Bx Bx
∂ w
+ 0 0 0 u − ρ 0 ρ 0 ∂
w =0
∂t
Bx
0 0 0 0 0 0 0 0
∂x
Bx
By
0 By −Bx 0 −v u 0 0
By
Bz 0 Bz 0 −Bx −w 0 u 0 Bz
p 0 γp 0 0 (γ − 1)v · B 0 0 u p
(3.122)
ρ u ρ 0 0 0 0 0 0
By
u 0 u 0 0 − Bρx Bz 1
ρ ρ ρ
By Bx
v
0 0 u 0 − ρ ρ 0 0
− Bρz Bx
w 0 0 0 u 0 ρ 0
V= ; Ap = (3.123)
Bx
0 0 0 0 0 0 0 0
By
0 By −Bx 0 −v u 0 0
Bz 0 Bz 0 −Bx −w 0 u 0
p 0 γp 0 0 (γ − 1)v · B 0 0 u
Los valores propios del sistema primitivo son idénticos a los del sistema conservativo, y las
matrices de transformación entre ambos sistemas, dada por la relación (2.21) son:
1 0 0 0 0 0 0
u ρ 0 0 0 0 0
v 0 ρ 0 0 0 0
P=
w 0 0 ρ 0 0 0
(3.124)
0 0 0 0 1 0 0
0 0 0 0 0 1 0
1 2 1
2 (u + v2 + w2 ) ρu ρv ρw By Bz γ−1
3.5. Sistema unidimensional 57
u ρ 0 0 0 0 0
2
u 2ρu 0 0 By Bz 1
uv ρv ρu 0 −Bx 0 0
uw ρw
Q= 0 ρu 0 −Bx 0
0 By −Bx 0 u 0 0
0
0 −Bz By w −v 0
γ
Q7,1 Q7,2 ρuv − Bx By ρuw − Bx Bz 2uBy − vBx 2uBz − wBx γ−1 u
(3.125)
donde
1 1 γ
Q7,1 = u(u2 + v 2 + w2 ) Q7,2 = ρ(u2 + v 2 + w2 ) + ρu2 + p + By2 + Bz2
2 2 γ−1
0 1 0 0 0 0 0
−2u2 + 1−γ
2 v
2
(3 − γ) u (1 − γ) v (1 − γ) w (2 − γ) By (2 − γ) Bz (γ − 1)
−uv v u 0 −Bx 0 0
Ac = −uw w 0 u −Bx 0
uBy By
− ρ + vBρ x − Bρx
ρ 0 u 0 0
− uB
ρ + ρ
z wBx Bz
0 − Bρx 0 u 0
ρ
Ac8,1 Ac8,2 Ac8,3 Ac8,4 Ac8,5 Ac8,6 γu
(3.127)
u ρ 0 0 0 0 0
By Bz 1
0 u 0 0
ρ ρ ρ
0 0 Bx
u 0 ρ
0 0
Ap = 0 0 Bx
0 u 0 0 (3.128)
ρ
0 By −Bx 0 u 0 0
0 Bz 0 −Bx 0 u 0
0 γp 0 0 0 0 u
Bx By Bz
bx = √ ; by = √ ; bz = √ ; (3.129)
ρ ρ ρ
primitiva:
3.5. Sistema unidimensional 59
" #
−cf −cf c2f by c2f bz 1
l±f = 0, ±cf , ± 2 2
bx by , ± 2 2
bx bz , 2 2
√ , 2 2
√ ,
(cf − bx ) (cf − bx ) (cf − bx ) ρ (cf − bx ) ρ ρ
(3.132)
ρ
±cf
−c
± 2 f 2 bx by
(cf −bx )
−cf
r±f = ± (c2f −b2x ) bx bz (3.133)
2 c f 2 by √ ρ
2
(cf −bx )
c2f √
(c2f −b2x ) bz ρ
ρa2
c2s c2s
−cs −cs by bz 1
l±s = 0, ±cs , ± 2 2
b x b y , ± 2 2
b x b z , 2 2
√ , 2 2
√ ,
(cs − bx ) (cs − bx ) (cs − bx ) ρ (cs − bx ) ρ ρ
(3.134)
ρ
±cs
± −cs b b
(c2s −b2x ) x y
−cs
r±s ± (c22s −b2x ) bx bz
= (3.135)
√
(c2c−b
s
2 ) by ρ
s x
c2s √
(c2s −b2x ) bz ρ
ρa2
Bz By
l±a = 0, 0, −Bz , By , ± √ , ∓ √ , 0 (3.136)
ρ ρ
60 Capítulo 3. Modelos Físicos
0
0
−Bz
r±a = By (3.137)
√
± ρBz
√
∓ ρBy
0
λe = u (3.138)
−1
le = 1, 0, 0, 0, 0, 0, 2 (3.139)
a
1
0
0
re = 0 (3.140)
0
0
res propios asociados a ondas diferentes se aproximan entre sí, y sus vectores propios
asociados (en el caso de que dichos valores propios coincidan) se vuelven singulares,
ρ u ρ 0 0 0 0 0 ρ
By Bz 1
u
0 u 0 0 ρ ρ ρ
u
Bx
v 0 0 u 0 0 0 v
ρ
∂
Bx
∂
w + 0 0 0 u 0 0 w =0 (3.141)
∂t ρ ∂x
By 0 By −Bx u By
0 0 0
Bz
0 Bz 0 −Bx 0 u 0
Bz
p 0 γp 0 0 0 0 u p
Teniendo en cuenta las Ecs. (3.129), (3.130), y (3.131) puede obtenerse una nueva
matriz jacobiana cuyos vectores y valores propios ponen en evidencia casos de inde-
terminación. Los valores propios del sistema son, usando la nomenclatura propuesta
λ1 = u − cf ; λ2 = u − bx λ3 = u − cs λ4 = u
(3.142)
λ5 = u + cs λ6 = u + bx λ7 = u + cf
Donde a es la velocidad del sonido gasdinámica, dada por la Ec. (3.42). Luego,
resolviendo esta ecuación de cuarto grado pueden expresarse las velocidades magne-
tosónicas como:
s 2
1 γp 1 γp γp
c2f = 2
+b + +b 2 − 4 b2x
2 ρ 2 ρ ρ
s 2
1 γp 1 γp γp 2
c2s = + b2 − + b2 −4 b
2 ρ 2 ρ ρ x
Además, dichas soluciones cumplen las siguientes identidades:
c4f − a2 b2x = c2f c2f − c2s ; c4s − a2 b2x = c2s c2s − c2f
(3.144a)
62 Capítulo 3. Modelos Físicos
Se observa que las curvas de cf y cs constantes son elipses e hipérbolas cofocales res-
son bx /a y b⊥ /a .
s 2
γp γp γp 2
2c2f − + b2 = + b2 −4 b
ρ ρ ρ x
2 2
a2 a2
bx b⊥
+ 2 =1
c2f
a cf − a2 a
De la misma manera para la onda magnetosónica lenta, se obtiene la expresión para
2 2
a2 a2
bx b⊥
− 2 =1
c2f
a c f − a2 a
b y = bz = 0
ocurre:
1 1p 2
c2f = (a2 + b2x ) + (a − b2x )2 = a2
2 2
1 1p 2
c2s = (a2 + b2x ) − (a − b2x )2 = b2x
2 2
Luego, cuando bx = a ambas velocidades magnetosónicas serán iguales, lo que produ-
ce que el sistema de vectores propios se vuelva indeterminado. Los vectores propios
asociados a las ondas magnetosónicas son los que resultan entonces indeterminados,
hiperbólico.
3.5. Sistema unidimensional 63
2.5
cs=0.05 cs=0.15
cs=0.25
cs=0.35
2 cs=0.45
cs=0.55
cs=0.65
1.5 cs=0.75
bn /a
cs=0.85
cs=0.95
0.5 cf=2.405
cf=2.205
cf=2.005
cf=1.805
cf=1.605
cf=1.205 cf=1.405
cf=1.005
0
0 0.5 1 1.5 2 2.5
bx/a
a una misma onda característica, dados por las Ecs. (3.132), (3.133),(3.134) y (3.135)
, se obtiene:
2c2f (c2f − c2s )
Lf · Rf =
(c2f − b2x )
(3.145)
2c2s (c2f − c2s )
Ls · Rs =
(c2s − b2x )
Estas expresiones pueden ser singulares o indeterminadas en algunos casos. Por
esta razón, y para garantizar que los vectores propios sean ortonormales, es necesario
multiplicar a dichos vectores por sendos factores kl y kr , cuyo producto sea igual a
(c2 − b2x )
kl kr =
2c2 (c2f − c2s )
a2 − c2s c2f − a2
αf2 = αs2 = (3.146)
c2f − c2s c2f − c2s
64 Capítulo 3. Modelos Físicos
αf2 + αs2 = 1
Los autores normalizan los vectores izquierdos por 1/a2 para que tengan las
1 αf βy a −αf βz a αs
ls = 2 0, ±αs cs , ±αf cf βy SBx , ±αf cf βz SBx , − √ , √ , (3.147a)
2a ρ ρ ρ
1 αs βy a αf βz a αf
lf = 2 0, ±αf cf , ∓αs cs βy SBx , ∓αs cs βz SBx , √ , √ , (3.147b)
2a ρ ρ ρ
αs ρ αf ρ
±αs cs ±αf cf
±αf cf βy SBx ∓αs cs βy SBx
rs = ±αf cf βz SBx rf = ∓αs cs βz SBx (3.148)
−αf √ρβy a αs √ρβy a
√ √
−αf ρβz a αs ρβz a
αs ρa2 αf ρa2
donde
by bz
βy = , βz = , SBx = signo(Bx ) (3.149)
b⊥ b⊥
Sin embargo, estas expresiones son indeterminadas cuando b⊥ = 0 , es decir cuando
permanecen sin cambios. A pesar de ello, a los vectores de las ondas de Alfvén
anteriormente denidos:
1 βz SBx βy SBx
la = 0, 0, ±βz , ∓βy , − √ , √ , 0 (3.150)
2 ρ ρ
0
0
±βz
ra = ∓βy (3.151)
√
−βz SBx ρ
√
βy SBx ρ
0
Casos de indeterminación
Existen seis casos posibles, en función de los valores que puedan tomar las varia-
bles a, b x , y b⊥ .
1. Si bx ≈ 0 pero b⊥ 6= 0, los coecientes αf y αs quedan:
a2 b2x
c2f ≈ a2 + b2⊥ ; c2s ≈ (3.152)
a2 + b2⊥
a2 b2⊥
αf2 ≈ ; αs2 ≈ ;
a2 + b2⊥ a2 + b2⊥
αf2 ≈ 1; αs2 ≈ 0;
αf2 ≈ 1; αs2 ≈ 0;
66 Capítulo 3. Modelos Físicos
De la misma forma que para el caso anterior, las ondas magnetosónicas rápidas
tienden a ondas sónicas, pero las ondas lentas tienden a ondas de Alfvén
4. Si |bx | > a + y b⊥ ≈ 0
c2f ≈ b2x ; c2s ≈ a2 (3.155)
αf2 ≈ 0; αs2 ≈ 1
En este caso ocurre lo inverso para el caso anterior, las ondas lentas tienden a
magnetosónicas se superponen:
minados
αf2 ≈ 0; αs2 ≈ 0
Sin embargo, como cumplen con la condición αs2 + αf2 = 1, los vectores pro-
bx ≈ (a + 1)signo(bx ); b⊥ = 2
b⊥ b⊥
! !
atan( ) atan( )
|bx |−a |bx |−a
αf ≈ sin + δf ; αs ≈ cos + δs ;
2 2
b⊥
donde además demostraron que |δf,s | 6 4a
.
Esquema Numérico
de estado por funciones constantes a trozos. Es decir, se asume que las variables
sistema no conservativo, dada por la Ec. (2.6) puede escribirse para cada celda Ωj
en la forma discreta como:
∆Uj 1 h i
=− F 1 − Fj− 1 + S(U)j (4.1)
∆t ∆x j+ 2 2
función de ujo numérico, evaluada en la interfaz de dos celdas contiguas, que debe
Dicha función, por lo tanto debe ser consistente con el ujo del sistema continuo,
y además debe cumplir con otras propiedades matemáticas para garantizar la es-
67
68 Capítulo 4. Esquema Numérico
presente tesis.
∂u ∂f (u)
+ =0 (4.2)
∂t ∂x
u(x, 0) = φ(x)
como: Z
∂u
T V = dx (4.3)
∂x
Si la Variación Total no aumenta en el tiempo
crearán nuevos extremos locales, los máximos existentes no aumentarán, y los míni-
cercanas a discontinuidades.
Asumiendo además que existe una función de entropía (ya sea física o numérica)
η(u) con un ujo asociado Ψ(u) tal que también cumple la ley de conservación:
∂η(u) ∂Ψ(u)
+ =0 (4.5)
∂t ∂x
de la forma:
∆t
un+1 = unj −
j f j+1/2 − f j−1/2 (4.6)
∆x
Un esquema numérico general de 2k + 1 celdas puede expresarse , mediante un
±
operador diferencial discreto L(u) en función de coecientes Cj+1/2 :
−
un+1
j
+
= Lun = unj + Cj+1/2 ∆j+1/2 u − Cj−1/2 ∆j−1/2 u (4.7)
+ −
Los coecientes Cj+1/2 y Cj−1/2 son en general función de los valores de la variable
u en otras celdas, esto permite que el esquema numérico pueda ser no lineal.
−
+
Cj+1/2 = C + (uj−k+1 , ..., uj+k ); Cj−1/2 = C − (uj−k , ..., uj+k−1 ) (4.8)
+ − + −
Cj+1/2 ≥ 0; Cj−1/2 ≥ 0; Cj+1/2 + Cj−1/2 ≤1 (4.9)
1 ∆x ∆t
f j+1/2 = f (uj+1 ) + f (uj ) − ψ aj ∆j+1/2 u (4.10)
2 ∆t ∆x
f (uj+1 )−f (uj ) si ∆j+1/2 u 6= 0
∆j+1/2 u
aj+1/2 = (4.11)
a(x
j+1/2 ) si ∆j+1/2 u = 0
el esquema resultante es TVD, siempre que la función ψ(x) sea una función escalar
Sin embargo, este esquema es de primer orden de precisión. Para lograr una apro-
ximación de segundo orden, Harten se vale de que se demuestra que una aproximación
∂u ∂f (u) ∆x ∂g(u)
+ =
∂t ∂x ∆t ∂x
Pasando al primer miembro el término que contiene a g(u) se obtiene una ley de
∂u ∂ ∆x
+ f (u) − g(u) =0 (4.13)
∂t ∂x ∆t
g(u) = O(∆x)
1h
f j+1/2 = f (uj+1 ) + f (uj )+
2
∆x ∆x ∆t
(gj + gj+1 ) − ψ aj+1/2 + γj+1/2 ∆j+1/2 u (4.14)
∆t ∆t ∆x
ciones:
" 2 #
∆t ∆t
gj+1 + gj = ψ aj+1/2 − aj+1/2 ∆j+1/2 u + O(∆x2 )
∆x ∆x
Entonces el esquema conservativo (Ec. (4.6)) con el ujo dado por la Ec. (4.14)
donde
1
σ(x) = ψ(x) − x2 > 0 (4.17)
2
O en su forma discreta
∆t
Sj+1/2 = signo σ aj+1/2 ∆j+1/2 u (4.18)
∆x
Por otro lado, el ujo de Harten cumple con la propiedad de preservar la po-
1991).
(Yee, 1987) introdujo una modicación en el ujo de Harten para volverlo más
exible en la elección de la función limitadora a emplear, y para lograr que sea menos
difusivo. Dicho ujo modicado se conoce como ujo de Harten-Yee. Sacando como
72 Capítulo 4. Esquema Numérico
1
f j+1/2 = f (uj+1 ) + f (uj ) + σ(aj+1/2 )(gj + gj+1 )−
2
∆t
ψ aj+1/2 + γj+1/2 ∆j+1/2 u (4.19)
∆x
∆t
gj = Sj+1/2 max 0, min
∆x a j+1/2 ∆j+1/2 u
, Sj+1/2 λaj−1/2 ∆j−1/2 u (4.20)
positividad en regiones de baja densidad, lo que lo hace una opción adecuada para
R−1 AR = Λ
W = R−1 U (4.21)
Wt + ΛWx = 0 (4.22)
dichas ecuaciones desacopladas. La forma nal del ujo numérico TVD de Yee que
4.2. Método TVD de Harten para una ley de conservación unidimensional 73
variables originales.
1
Fi+1/2 = Fi+1 + Fi + Ri+1/2 Φi+1/2 (4.23)
2
expresan como:
(l) (l) (l) (l) (l) (l) (l) (l)
φi+1/2 =σ (λi+1/2 )(gi+1 + gi ) −ψ λi+1/2 + γi+1/2 αi+1/2 (4.24)
la siguiente manera:
nh io
(l) (l) (l) (l) (l)
gi = Si+1/2 max 0, min (λi+1/2 )|αi+1/2 |, (λi−1/2 )αi−1/2 Si+1/2 ) (4.25)
(l)
Si+1/2 = signo(αi+1/2 ) (4.26)
(l) (l)
gj+1(l)−gj
si
(l)
αj+1/2 6= 0
(l) αj+1/2
γj+1/2 = (4.27)
0 (l)
si αj+1/2 =0
(l) (l) 1 (l) ∆t (l) 2
σ (λi+1/2 ) = ψ(λi+1/2 ) − (λ ) (4.28)
2 ∆x i+1/2
(l)
Y la función ψ(λi+1/2 ) es la función encargada de introducir la disipación numérica
propiamente dicha. Para respetar la condición de entropía y lograr la convergencia
a una solución físicamente correcta (Harten y Hyman, 1983), esta función se dene
como:
(l) (l)
(l)
|λi+1/2 | si |λi+1/2 | ≥ δ1
ψ(λi+1/2 ) =
(l)
(l) (4.29)
(λi+1/2 )2 + δ12 /2δ1 si |λi+1/2 | < δ1
la gravedad
La función de calentamiento EH
unidimensional como:
0
S=
−ρg
(4.30)
2 ∂ 2 T 7/2
−ρgu + κsp + EH + Lrad
7 ∂x2
Z
1
S(U)dΩ ≈ Sj (Uj ) (4.31)
Ωj Ωj
Z Z
∇ · (κ∇T ) dΩ = n · (κ∇T ) dA (4.33)
Ωj Aj
Z
∂ ∂T ∂T ∂T
κ dΩ = κA − κA (4.34)
Ωj ∂x ∂x ∂x j−1/2 ∂x j+1/2
ción térmica:
5/2 5/2
κj−1 + κj Tj−1 + Tj
κj−1/2 = = κsp
2 2
5/2 5/2
κj+1 + κj Tj + Tj+1
κj+1/2 = = κsp (4.35)
2 2
∂T Tj−1 − Tj
Fc j−1/2 = κA = κj−1/2 Aj−1/2
∂x j−1/2 ∆x
∂T Tj − Tj+1
Fc j+1/2 = κA = κj+1/2 Aj+1/2 (4.36)
∂x j+1/2 ∆x
Fc j−1/2 − Fc j+1/2
Q̇cond = (4.37)
∆x
D∆t
FO = (4.38)
∆x2
76 Capítulo 4. Esquema Numérico
κ
donde D= Cp ρ
es la difusividad térmica, que se estima con los valores de la celda
que el esquema sea estable en el tiempo y para tratar de manera consistente los
∂
[F(U, x)] = S(U, x) (4.39)
∂x
diendo del instante de tiempo en el que se evalúen las funciones de ujos numéricos y
Un+1 − Un ∂F(Un )
+ − S(Un ) = 0 (4.40)
∆t ∂x
Un+1 − Un ∂F(Un+1 )
+ − S(Un+1 ) = 0 (4.41)
∆t ∂x
la mayoría de las veces esto no puede calcularse de forma exacta, es necesario hacer
(n) (n)
(n+1) (n+1) ∂Fj+1/2 (Uj , Uj+1 ) ∂Uj
Fn+1
j+1/2 (Uj , Uj+1 ) ≈ Fj+1/2 (Unj , Unj+1 ) + ∆t
∂Uj ∂t
(n) (n)
∂Fj+1/2 (Uj , Uj+1 ) ∂Uj+1
+ ∆t + O(∆t2 ) (4.43)
∂Uj+1 ∂t
∂Sn ∂Uj
Sn+1
j (Un+1 n+1
j−1 , Uj , Un+1 n
j+1 ) ≈ S + ∆t
∂Uj ∂t
∂Sn ∂Uj+1 ∂Sn ∂Uj−1
+ ∆t + ∆t + O(∆t2 ) (4.44)
∂Uj+1 ∂t ∂Uj−1 ∂t
∂Fnj+1/2 ∂Fnj+1/2
= Acj (Uj , Uj+1 ); = A+
j (Uj , Uj+1 ) (4.45)
∂Uj ∂Uj+1
∂Snj
= AScj (Uj−1 , Uj , Uj+1 ) (4.46)
∂Uj
∂Snj
= AS+
j (Uj−1 , Uj , Uj+1 ) (4.47)
∂Uj+1
∂Snj
= AS−
j (Uj−1 , Uj , Uj+1 ) (4.48)
∂Uj−1
Pero, por otro lado, los métodos implícitos tienen condiciones de estabilidad para la
integración temporal menos severas y son menos sensibles a errores numéricos cuando
notación
∆Uj = Un+1
j − Unj (4.49)
∆Uj Acj A+
j Acj A+j
+ ∆Uj + ∆Uj+1 − ∆Uj−1 − ∆Uj (4.50)
∆t ∆x ∆x ∆x ∆x
−
F (Unj+1/2 − F(Unj−1/2 )
− AScj ∆Uj + AS+ + S(Unj )
j ∆U j+1 + AS j ∆U j−1 ∆t = −
∆x
∆t c ∆t + c
I+ (A − (A − ASj ∆t ∆Uj
∆x j ∆x j−1
∆t + + ∆t c −
+ A − ASj ∆t ∆Uj+1 − A − ASj ∆t ∆Uj−1 =
∆x j ∆x j−1
∆t n
Fj−1/2 − Fnj+1/2 + Sn ∆t
(4.51)
∆x
0 0 0
−g
c
ASj =
0 0
n−1 2,5 n−1 2,5
c n ρj Tj uj K 2 κsp Tj uj n ρj Tj K 2 κsp Tj
ASj 3 3 − 2 2
+ 2 −g 2 2
− 2
Cv mh µ Cv ∆x ρj Cv mh µ Cv ∆x ρj
(4.52a)
n−1 n
2
c n ρ T
j j uj − 2 C v (n − 2) ρ T
j j K κsp Tj uj − 2 Cv κsp Tj 3,5
2,5 2
ASj 3 3 = −
2 Cv mh 2 µ2 Cv ∆x2 ρj
(4.52b)
4.4. Integración temporal 79
0 0 0
AS+ =
0 0 0
j 2,5 3,5
2
κsp Tj+1 2,5 uj+1 κsp Tj+1 2,5
κsp Tj+1 uj+1 − 2 Cv κsp Tj+1
−
2 Cv ∆x2 ρj+1 Cv ∆x2 ρj+1 Cv ∆x2 ρj+1
(4.53)
0 0 0
AS− =
0 0 0
j 2,5 3,5
2
κsp Tj−1 2,5 uj−1 κsp Tj−1 2,5
κsp Tj−1 uj−1 − 2 Cv κsp Tj−1
−
2 Cv ∆x2 ρj−1 Cv ∆x2 ρj−1 Cv ∆x2 ρj−1
(4.54)
donde se empleó una función Λ(T ) genérica, que describe en general a todas las
ecuaciones asociadas a la función de pérdidas por radiación (Ec. (3.2)), para obtener
Λ(T ) = KT n
Por otro lado, si empleamos una discretización del término de conducción de calor
por volúmenes nitos, como la dada por la Ec. (4.37), las expresiones analíticas de las
Otra opción posible es emplear una evaluación numérica aproximada para dichas
en la integración espacial sin presentar oscilaciones, debe vericarse que los términos
convectivos cumplan la propiedad TVD. (Yee et al., 1985) demostraron que para una
∂u ∂f (u)
+ =0
∂t ∂x
80 Capítulo 4. Esquema Numérico
∆t n+1 n+1
∆t n n
un+1
j + η f j+1/2 − f j−1/2 = unj − (1 − η) f j+1/2 − f j−1/2
∆x ∆x
con el ujo numérico f j+1/2 de Harten (Ec. (4.14)) es TVD de segundo orden en el
∆t
≤ ∆t ψ λj+1/2 ≤ 1
λj+1/2
∆x
∆x 1−η
demuestra para este último caso que dicho esquema es incondicionalmente estable
∆t n+1 n+1
un+1
j + η f j+1/2 − f j−1/2 = unj
∆x
∆t − ∆t +
un+1
j − C (λ + γ)n+1
j+1/2 ∆j+1/2 u
n+1
+ C (λ + γ)n+1
j−1/2 ∆j−1/2 u
n+1
= unj (4.55)
∆x ∆x
± 1
Cj+1/2 = ψ(λj+1/2 + γj+1/2 ) ± (λj+1/2 + γj+1/2 ) (4.56)
2
Donde las expresiones de λ, γ y gj están dadas por las Ecs. (4.11), (4.15) y (4.16)
respectivamente.
1 ∆t 1
σ(λ) = ψ(λ) + η− λ2
2 ∆x 2
plearse:
1
σ(λ) = [ψ(λ)]
2
debido a que el esquema es más estable, permitiendo pasos de integración mayores.
Debe notarse que el esquema dado por la Ec. (4.55) es un sistema no lineal de
C − (λ + γ)n+1 − n
j+1/2 ≈ C (λ + γ)j+1/2 (4.57)
De esta manera se obtiene el sistema linealizado más simple posible, pero se pierde
dj = un+1
j − unj (4.58)
∆t n n
E1 dj−1 + E2 dj + E3 dj+1 = − fj+1/2 − fj−1/2 (4.59)
∆x
82 Capítulo 4. Esquema Numérico
donde
∆t −
E1 = − C (λ + γ)nj+1/2
∆x
∆t −
C (λ + γ)nj+1/2 + C + (λ + γ)nj+1/2
E2 = 1 + (4.60)
∆x
∆t +
C (λ + γ)nj+1/2
E3 = −
∆x
El esquema es independiente del valor del paso de tiempo ∆t para los términos
convectivos
En la presente tesis se empleó el ujo de Harten Yee (Ec. (4.19)) en lugar del
ujo original de Harten (Ec. (4.14)), ya que el primero, al ser una generalización del
generalizan mediante:
∆t −
E1 = − J (λ + γ)nj+1/2
∆x
∆t −
J (λ + γ)nj+1/2 + J+ (λ + γ)nj+1/2
E2 = 1 + (4.61)
∆x
∆t +
J (λ + γ)nj+1/2
E3 = −
∆x
Donde
J± n ± l l n −1 n
j+1/2 = Rj+1/2 diag(C (λ + γ )j+1/2 )(R )j+1/2 (4.62)
como los coecientes αi de expansión espectral. Para el caso del modelo de Euler, se
utilizó el solver de Roe convencional (Roe, 1981), y para el caso MHD se empleó el
el estado izquierdo y derecho. Dicha matriz debe poseer las siguientes propiedades
(Roe, 1981):
mente independientes.
Para el caso de un sistema MHD unidimensional (de siete ondas), el estado in-
U(UR , UL ) = U(ρ, u, v, w, H ∗ , By , Bz )
84 Capítulo 4. Esquema Numérico
de Roe: √ √ √ √
ρL ξR + ρR ξL ρR ξR + ρL ξL
ξ= √ √ ; ξ= √ √
ρR + ρL ρR + ρL
A su vez, se dene un incremento linealizado del producto entre dos variables como:
donde ∆ξ = ξL − ξR .
Para el caso general del problema de Riemann de ocho ondas no existe un solver
de Roe general, por lo que varios autores utilizan distintos tipos de promedios para
una matriz de Roe para el sistema MHD unidimensional isoentrópico, que a su vez
∂u ∂ρu
+ =0 (4.63a)
∂t ∂x
∂ ∂
ρu2 + p + B 2 /2 = 0
(ρu) + (4.63b)
∂t ∂x
∂B ∂Bu
+ =0 (4.63c)
∂t ∂x
Cuya matriz jacobiana es:
0 1 0
2 2
−u + a + X 2u + Y Z (4.64)
−u Bρ B
ρ
u
derecha:
2 BR + BL
∆B = 2 ∆B
2
El sistema de ecuaciones puede dar valores propios imaginarios. Por esta razón los
B2
∆ = X∆ρ + Y ∆(ρu) + Z∆B (4.65)
2
donde los coecientes X, Y y Z deben evaluarse. Como se espera que el sistema sea
invariante en una transformación galileana y que los valores propios sean simétricos
respecto de la velocidad u, el término Y debe ser igual a 0. Por otro lado, para obtener
se obtiene que
(∆B)2
X= √ √ 2 (4.66)
2 ρL + ρR
lo que permite expresar el promedio de Roe para el salto de la presión magnética
como:
B2
∆ = X∆ρ + B∆B (4.67)
2
A partir de este resultado, puede obtenerse la matriz jacobiana de Roe para el sistema
Donde
γ−1 2
A2,1 = −u2 + (2 − γ)X + v
2
Bx
A7,1 = uH∗ + u(A2,1 + u2 ) + (B · v)
ρ
Bx2
A7,2 = H∗ + u(A2,2 − 2u) −
ρ
Bx 2
A7,3 = uA2,3 − By
ρ
Bx 2
A7,4 = uA2,4 − Bz
ρ
A7,5 = uA2,5 − Bx v
A7,6 = uA2,6 − Bx w
A7,7 = u + uA2,7
λ 1 = u − cf λ 2 = u − ca λ 3 = u − cs λ 4 = u
(4.69)
λ 5 = u + cs λ 6 = u + ca λ 7 = u + cf
donde:
2
c a 2 = bx (4.70)
q
1 2
cf 2 = a∗2 + a∗4 − 4a2 bx (4.71)
2
q
1 2 2
cs 2 = 2 4
a∗ − a∗ − 4a bx (4.72)
2
y " #
Bx By Bz
b= √ , √ , √
ρ ρ ρ
(4.73)
2
2 2
a = (2 − γ)X + (γ − 1) H∗ − v /2 − b
2
a∗2 = a2 + b
4.5. Solver de Riemann aproximado de Roe para MHD 87
ραf
ραf (u − cf )
ρ(αf v + αs cs βy S)
1
R(1) = 2 ρ(αf w + αs cs βz S)
ρa √
ραs aβy
√
ραs aβz
√
h i
B2
ραf (H∗ − ρ
− ucf + ραs cs S(βy v + wβz ) + ραs a B⊥ )
(4.74)
0
0
−ρβz
R(2) = ρβy (4.75)
√
−S ρβz
√
S ρβy
−ρ(vβz − wβy )
ραs
ραs (u − cs )
ρ(αs v − αf cf βy S)
1
R(3) = 2 ρ(αs w − αf cf βz S) (4.76)
ρa √
− ραf aβy
√
i − ραf aβz
√
h
B2
ραs (H∗ − ρ
− ucs − ραf cf S(βy v + wβz − ραf a B⊥ )
88 Capítulo 4. Esquema Numérico
1
u
v
1 w
R(4) = 2 (4.77)
a
0
0
γ−2
v2 /2 + X
γ−1
ραs
ραs (u + cs )
ρ(αs v + αf cf βy S)
1
R(5) = 2 ρ(αs w + αf cf βz S) (4.78)
ρa √
− ραf aβy
√
i − ραf aβz
√
h
B2
ραs (H∗ − ρ
+ ucs + ραf cf S(βy v + wβz − ραf a B⊥ )
0
0
ρβz
R(6) = −ρβy (4.79)
√
−S ρβz
√
S ρβy
ρ(vβz − wβy )
ραf
ραf (u + cf )
ρ(αf v − αs cs βy S)
1
R(7) = 2 ρ(αf w − αs cs βz S) (4.80)
ρa √
ραs aβy
√
ραs aβz
√
h i
B2
ραf H∗ − ρ
+ ucf − ραs cs S(βy v + wβz + ραs a B⊥ )
4.5. Solver de Riemann aproximado de Roe para MHD 89
1
α1 = [αf (X∆ρ + ∆p) + ραs cs S(βy ∆v + βz ∆w) (4.81a)
2
−ραf cf ∆u + ραs a(βy ∆By + βz ∆Bz )]
p
1 S
α2 = [βy ∆w − βz ∆v + √ (βy ∆Bz − βz ∆By ) (4.81b)
2 ρ
1
α3 = [αs (X∆ρ + ∆p) − ραf cf S(βy ∆v + βz ∆w) (4.81c)
2
−ραs cs ∆u − ραf a(βy ∆By + βz ∆Bz )]
p
1 S
α6 = [−βy ∆w + βz ∆v + √ (βy ∆Bz − βz ∆By ) (4.81f )
2 ρ
1
α7 = [αf (X∆ρ + ∆p) − ραs cs S(βy ∆v + βz ∆w) (4.81g)
2 p
+ραf cf ∆u + ραs a(βy ∆By + βz ∆Bz )]
Capítulo 5
Condiciones de Contorno
5.1. Introducción
Un modelado consistente y físicamente correcto de las condiciones de contorno
es fundamental para lograr la estabilidad numérica del modelo, así como para lograr
∂U ∂F
+ + C0 = 0 para xL ≤ x ≤ xR
∂t ∂x
Con una discretización del dominio físico en M celdas, y en los contornos con
una formulación general con k celdas, numeradas para el extremo izquierdo (x < xL )
(−k . . . − 1, 0), y para el extremo derecho (x > xR ) (M + 1, M + 2 . . . M + k). En la
Figura 5.1 se esquematiza el problema considerado.
que se propagan hacia afuera, y otras que se transmiten desde el contorno hacia el
interior del dominio. Las características que se propagan hacia afuera dependen sólo
de las variables en el interior del dominio (y por ello no pueden ser inuenciadas
por el valor de las variables en el contorno). Por esta razón, las ondas características
91
92 Capítulo 5. Condiciones de Contorno
λn-j
λn-j
λn-1 λ2 λn-1
λn
λ2 λn λ1
λ1
dominio
de cálculo V
V
xL xR
celda m celda m+1
celda 0 celda 1
las relaciones impuestas por las ondas características salientes en cada instante.
De la misma forma, las características entrantes sólo dependen del valor de las
dominio. Por ello deben ser establecidas a priori de acuerdo a la física del problema
que se desee estudiar, como se expresa en los libros de Homann y Chiang (2000b)
y Blazek (2005).
las restantes n−b variables para poder evaluar ujos, u otras operaciones necesa-
rias para cálculos dentro del dominio. Por estas razones, algunos autores hablan de
La forma más simple de estimar las variables asociadas a las características salien-
tes es mediante una extrapolación de orden cero o de orden uno. Esta extrapolación
diferencias nitas explícitos o implícitos. Una discusión bastante completa sobre dis-
(Yee, 1981). Para imponer BCs con extrapolación de orden cero en el espacio para
5.2. Condiciones de contorno para el modelo gasdinámico 93
el de las funciones limitadoras gi (dadas por la Ec. (4.25)) en las celdas que colindan
U = U ; g0 = g1 para x < xL
0 1
(5.1)
U
m+1 = Um ; gm+1 = gm para x > xR
que
U = U
m m+1
(5.2)
U =U
m+2 m−1
Sin embargo, esta forma en general no suele ser consistente con la física del
variables primitivas:
∂V ∂V
+A +C=0 (5.3)
∂t ∂x
y teniendo en cuenta la transformación entre el sistema conservativo y primitivo,
∂U ∂V
=P
∂t ∂t
∂F ∂V
=Q (5.4)
∂x ∂x
luego
A = P−1 Q
C = P−1 C0
Premultiplicando la Ec. (5.3) por los vectores propios izquierdos l(i) puede expre-
sarse la evolución del sistema en función de los vectores propios izquierdos y de las
velocidades características:
∂V ∂V
l(i) · + λ(i) l(i) · + l(i) · C = 0 (5.5)
∂t ∂x
ρ 0 −a2 0
T T T
V = u ; l(1) = −ρa ; l(2) = 0 l(3) = ρa
p 1 1 1
En el contorno izquierdo las ondas con λ(i) ≤ 0 serán ondas salientes y las ondas
(i)
con λ ≥0 serán entrantes. En el contorno derecho es a la inversa. El criterio de
las amplitudes de las ondas entrantes debe ser constante en el tiempo. Deniendo
5.2. Condiciones de contorno para el modelo gasdinámico 95
ticas
∂Z(i) (i) ∂Z
(i)
+λ =0 (5.7)
∂t ∂x
Se observa que para este sistema la amplitud de cada onda, dada por Z(i) es constante
dx
a lo largo de una curva C en el plano x−t denida por
dt
= λ(i) . Debe tenerse en
cuenta, sin embargo, que esta parametrización con la variable Z es exacta sólo si
ser una relación exacta, aún así este método provee de una buena estimación para
Luego, si para imponer las condiciones de contorno se desea que las amplitudes de
las ondas entrantes al dominio sean iguales a cero en el contorno, basta con establecer
Para determinar cuáles ondas son entrantes y cuáles salientes, basta con evaluar
saliente al dominio; y si fuera mayor a cero, se trata de una onda entrante. Se observa
que las ondas λ(1) , λ(2) . . . λ(j) son salientes y poseen velocidades menores a cero (ya
∂V ∂V
l(i) ·
∂t
+ λi l(i) · ∂x
+ l(i) · C x=x = 0 si λ(i) ≤ 0
L
(5.9)
l(i) · ∂V
∂t
+ l(i) · C x=x = 0 si λ(i) > 0
L
∂p ∂u
− ρa + L(1) = 0 (5.11a)
∂t ∂t
∂p ∂ρ
− a2 + L(2) = 0 (5.11b)
∂t ∂t
∂p ∂u
+ ρa + L(3) = 0 (5.11c)
∂t ∂t
(1) (1) ∂p ∂u (2) (2) ∂p ∂ρ ∂p ∂u
L =λ − ρa ; L =λ − a2 ; L (3)
=λ (3)
+ ρa
∂x ∂x ∂x ∂x ∂x ∂x
(5.12)
Pueden estimarse los operadores L(i) en las celdas j−ésimas de los contornos mediante
5.2. Condiciones de contorno para el modelo gasdinámico 97
derivadas descentradas :
L(1) = (uj − aj ) 1 [pj+1 − pj − ρj aj (uj+1 − uj )] , j < 1, uj − aj < 0
j ∆x
(5.13a)
L(1) = (u − a ) 1 [p − p
j j j ∆x j j−1 − ρj aj (uj − uj−1 )] , j > M, uj − aj > 0
h i
L(2) = uj 1 pj+1 − pj − a2 (ρj+1 − ρj ) , j < 1, uj < 0
j ∆x j
h i (5.13b)
L(2) = u 1 p − p − a 2 (ρ − ρ
j j ∆x j j−1 j j j−1 , j > M,
) uj > 0
L(3) = (uj + aj ) 1 [pj+1 − pj + ρj aj (uj+1 − uj )] , j < 1, uj + aj < 0
j ∆x
(5.13c)
L(3) = (u + a ) 1 [p − p + ρ a (u − u )] , j > M, uj + aj > 0
j j j ∆x j j−1 j j j j−1
(i)
Lj =0 en los demás casos (5.13d)
Una vez determinados los valores del operador L en cada caso, se puede evaluar la
Estas relaciones permiten, mediante el mismo integrador temporal utilizado para las
celdas del dominio, obtener la evolución temporal de las variables primitivas en el contorno.
∂ρu ∂ρ ∂u
=u +ρ
∂t ∂t ∂t
∂et 1 ∂ρ ∂u 1 ∂p
= u2 + ρu +
∂t 2 ∂t ∂t γ − 1 ∂t
98 Capítulo 5. Condiciones de Contorno
∆ρj 1 (2) 1 (3) (1)
= Lj − Lj + Lj (5.15a)
∆t a2i 2
∆ (ρu)j
uj (2) 1 (3) (1) 1 (3) (1)
= Lj − L j + L j − Lj − L j (5.15b)
∆t a2j 2 2aj
1 u2j
∆et j (2) 1 (3) (1) uj (3) (1)
= Lj − Lj + L j − Lj − L j (5.15c)
∆t 2 a2j 2 2aj
1
(3) (1)
− Lj + Lj
2(γ − 1)
los términos del miembro derecho con las variables de estado en el paso de tiempo n, y la
derivada temporal se calcula con una diferencia nita. Para utilizar el esquema de integra-
la matriz cuyas las son los vectores propios izquierdos del sistema en varaibles primitivas
y L al vector cuyos elementos son los L(i) , puede escribirse el sistema en forma matricial
∂U n+1
= − Pl−1 L + C0
(5.16)
∂t
Teniendo en cuenta las leyes de transformación para los vectores propios izquierdos y
derechos entre el sistema privitivo y el sistema conservativo, dadas por las Ecs. (2.12),
conservativo.
Los términos del miembro derecho se evalúan mediante una expansión en serie de Taylor,
de la misma forma que para la Ec. (4.43). Debe tomarse especial cuidado en evaluar las
derivadas en los nodos del stencil que usa esta metodología. Por un lado, los vectores propios
se evalúan sólo en la celda considerada. Pero por otro, los operadores L(i) se evalúan de la
misma manera que para el esquema explícito, en función de que las ondas sean entrantes
o salientes. Teniendo en cuenta que las derivadas espaciales se evalúan mediante esquemas
5.2. Condiciones de contorno para el modelo gasdinámico 99
L(i) = L(i) (Uj , Uj+1 ), para j<1
(5.18)
L(i) = L(i) (Uj−1 , Uj ), para j>M
Por estas razón se emplean dos tipos de extrapolaciones, según se trate del extremo izquierdo
o derecho:
∂ ∂Uj ∂ ∂Uj+1
[RL]n+1 ≈ [RL]n + (RL)n ∆t + (RL)n ∆t, para j<1 (5.19a)
∂Uj ∂t ∂Uj+1 ∂t
∂ ∂Uj ∂ ∂Uj−1
[RL]n+1 ≈ [RL]n + (RL)n ∆t+ (RL)n ∆t, para j>M (5.19b)
∂Uj ∂t ∂Uj−1 ∂t
misma manera que para el interior del dominio, o sea mediante la Ec. (4.44).
Luego, para integrar las condiciones de borde de forma implícita es necesario evaluar la
matriz jacobiana
∂
donde el operador (k) implica la derivada respecto a la k−ésima componente del vector
∂Uj
de variables conservativas evaluada en la j−ésima celda del dominio. Luego, la colum-
∂
(k)
RL =
∂Uj
∂L(1) ∂L(2) ∂L(3)
∂R11 (1) ∂R12 (2) ∂R13 (3)
(k) L + (k) R11 + (k) L + (k) R12 + (k) L + (k) R13
∂Uj ∂Uj ∂Uj ∂Uj ∂Uj ∂Uj
∂R21 (1) ∂L(1) ∂R22 (2) ∂L(2) ∂R23 (3) ∂L(3)
(k) L (k) R21 (k) L (k) R22 (k) L (k) R23
+ + + + + (5.21)
∂Uj ∂Uj ∂Uj ∂Uj ∂Uj ∂Uj
∂R31 (1) ∂L(1) ∂R32 (2) ∂L(2) ∂R33 (3) ∂L(3)
(k) L + (k) R31 + (k) L + (k) R32 + (k) L + (k) R33
∂Uj ∂Uj ∂Uj ∂Uj ∂Uj ∂Uj
El sistema de ecuaciones para integrar implícitamente las variables del contorno puede
100 Capítulo 5. Condiciones de Contorno
expresarse como:
h i
∂ ∂
I+ ∂Uj R(Uj )L(Uj , Uj+1 )∆t ∆Uj + ∂Uj+1 R(Uj )L(Uj , Uj+1 )∆t∆Uj+1 =
−R(Uj )L(Uj , Uj+1 )∆t para j < 1 (5.22)
h i
∂
I + ∂U j
R(Uj )L(Uj−1 , U j )∆t ∆Uj + ∂U∂j−1 R(Uj )L(Uj , Uj−1 )∆t∆Uj−1 =
−R(Uj )L(Uj , Uj−1 )∆t para j>M (5.23)
Como los operadores L(i) tienen dentro un operador diferencial discreto, no puede acep-
tarse que conmuten con el operador diferencial respecto de las variables de estado, empleado
∂ ∂
Jc = Rj L(Uj , Uj+1 ); J+ = Rj L(Uj , Uj+1 ); (5.24)
∂Uj ∂Uj+1
∂ ∂
Jc = Rj L(Uj , Uj−1 ); J− = Rj L(Uj , Uj−1 ); (5.25)
∂Uj ∂Uj−1
La obtención de esta matriz jacobiana exacta resulta relativamente sencilla para el mo-
delo anteriormente descrito. Sin embargo, para un modelo más complejo (como el modelo
MHD) las expresiones analíticas resultan demasiado extensas y complejas de derivar, in-
cluso mediante un procesador simbólico. Además, deben obtenerse 49 elementos para cada
una de las tres matrices jacobianas (en j − 1 , j,y j + 1 ) que se necesitan. Esta metodología
resulta entonces propensa a los errores, tanto en la obtención de las expresiones como en
su implementación en el algoritmo. Por esta razón, se optó por emplear para los casos más
Keyes, 2004), basada en las aproximaciones a productos de matrices y vectores usadas para
Dada una función vectorial F(U) y un vector v, este método permite obtener de manera
aproximada el producto de la matriz jacobiana asociada a F(U) y el vector v mediante:
donde ε es un parámetro que representa una pequeña perturbación. Debe elegirse cuida-
valor es muy pequeño, los errores de redondeo pueden llevar a resultados inexactos.
Para obtener el jacobiano aproximado se emplea una serie de vectores ek cuyos elementos
son iguales a 0, excepto el k−ésimo, que es igual a 1. Dichos vectores representan una
F(U + εk ek ) − F(U)
J · ek ≈ J(k) = (5.28)
εk
donde el parámetro εk se estima, de acuerdo a Knoll y Keyes (2004), para cada variable
(k)
Uj como:
(k)
εk = bUj +b (5.29)
√
donde b = O εtrunc . Para una computadora de 64 bits y doble precisión basta con hacer
b= 10−6 .
(k)
Sin embargo, en el caso de que Uj = −1 el valor de εk se anula, ocurriendo una
singularidad en el valor de los jacobianos aproximados. Esto ocurre por ejemplo con el
campo magnético transversal en el tubo de choque de Brio y Wu, caso que se discutirá más
(k)
adelante. Por esta razón modicamos esta fórmula tomando el valor absoluto de Uj = −1
(k)
εk = bkUj k + b (5.30)
Como para el modelo gasdinámico unidimensional de Euler pudimos obtener las matri-
ces jacobianas exactas, comparamos los resultados obtenidos con ambos métodos. Como se
Por otro lado, también se probó calcular las matrices jacobianas asociadas a los términos
fuente en el interior del dominio de cálculo con esta metodología. Como se obtuvieron muy
buenos resultados, con errores relativos máximos del orden del 2 % o 3 % comparados a la
expresión analítica de los casos más simples, optamos por emplear esta metodología para
del cambio de presión en el contorno debe viajar aguas arriba hasta que se encuentre
un equilibrio. Sin embargo, especicar más condiciones de contorno que las estrictamente
denidas por la física del problema produce oscilaciones e inestabilidades en el método. Con
que sólo puede conocerse para pocos problemas sencillos. En los trabajos de (Poinsot y
Lele, 1992) y de (Yee, 1981) se presenta una discusión muy completa sobre la física y el
Por estas razones, en un segundo trabajo de (Thompson, 1990) se generalizó esta me-
todología para distintos tipos de condiciones y para cualquier sistema de ecuaciones hiper-
bólico. En dicho trabajo se estudió el sistema de Euler 3D con términos fuente debidos a la
gravedad. En este caso, por la naturaleza de los estudios de esta tesis, se emplea un modelo
cuasi-unidimensional con área variable A(x) y con un vector de términos fuente genérico
cantidad de movimiento.
∂ρA ∂ρuA
+ + C 0(1) = 0 (5.31a)
∂t ∂x
∂ρuA ∂ρu2 A ∂pA
+ + + C 0(2) = 0 (5.31b)
∂t ∂x ∂x
∂ρet A ∂
+ [u (ρet A + pA)] + C 0(3) = 0 (5.31c)
∂t ∂x
∂ρA ∂ρA ∂u
+u + ρA + C (1) = 0 (5.32a)
∂t ∂x ∂x
∂u ∂u 1 ∂pA
+u + + C (2) = 0 (5.32b)
∂t ∂x ρA ∂x
∂pA ∂pA ∂u
+u + γpA + C (3) = 0 (5.32c)
∂t ∂x ∂x
En este caso los vectores propios izquierdos y los operadores L son análogos a los
del sistema unidimensional, con la diferencia que están valuados en función de las nuevas
(1) (1) ∂pA ∂u
L =λ − ρAa (5.33a)
∂x ∂x
∂pA ∂ρA
L(2) = λ(2) − a2 (5.33b)
∂x ∂x
(3) ∂pA ∂u
L(3) =λ + ρAa (5.33c)
∂x ∂x
dores Li como:
∂ρA 1 1 (3)
+ 2 −L(2) + L + L(1) + C (1) = 0 (5.34a)
∂t a 2
∂u 1
+ L(3) − L(1) + C (2) = 0 (5.34b)
∂t 2ρAa
∂pA 1 (3)
+ L + L(1) + C (3) = 0 (5.34c)
∂t 2
Además resultan útiles las siguientes relaciones para la variación temporal de las va-
riables conservadas dadas por la Ec. (5.15); y otras que relacionan la variación temporal y
104 Capítulo 5. Condiciones de Contorno
C (3)
∂T T γ − 1 (3)
+ 2 L(2) + L + L(1) = − (5.35a)
∂t ρa 2 ρR
" !#
∂T T L(2) γ − 1 L(3) L(1)
= 2 + + (5.35b)
∂x ρa u 2 u+a u−a
Trabajando algebraicamente con estas ecuaciones pueden obtenerse diferentes tipos de con-
Esta nueva manera de emplear los operadores L(i) se integra implícitamente empleando
en la sección anterior.
En este caso la condición está dada por u = 0 ∀t. Para el caso de x = xL , se tiene que
λ(1) = −a y L
(1) será una onda saliente, por lo tanto:
(1) (1) ∂pA ∂u
L =λ − ρAa
∂x ∂x
1 (3)
L − L(1) + C (2) = 0 → L(3) = L(1) − 2ρAaC (2) (5.36)
2ρAa
(3) (3) ∂pA ∂u
L =λ + ρAa
∂x ∂x
5.3. Generalización del Modelo 105
1 (3)
L − L(1) + C (2) = 0 → L(1) = L(3) + 2ρAaC (2)
2ρAa
5.3.2. Salida
u<0 para x = xL
u>0 para x = xR
Por lo tanto, los operadores L(2) y L(3) deberá obtenerse siempre mediante la relación para
ondas salientes dadas por la Ec. (5.13), siempre y cuando los términos fuente no incluyan
términos con derivadas de segundo orden en las variables conservadas (como por ejemplo
caso de que esto no ocurra, se debe determinar sólo el valor del operador L(1) . Dicho valor
Flujo supersónico
u<a<0 para x = xL
0<a<u para x = xR
Para x = xL todas las velocidades características son menores a cero. Por lo tanto, no
puede imponerse ninguna condición de contorno y los valores de los L(i) deben calcularse
Flujo subsónico
onda 3
(3) ∂W (3) (3) ∂pA ∂u
l · +L +l ·C = + ρAa + L(3) + ρAaC (2) − C (3) = 0
∂t x=xL ∂t ∂t
vimiento en variables primitivas (Ec. (5.32b)), si exigimos que la sumatoria de fuerzas sea
cero en el contorno, ello implica que un elemento de uido es transportado por advección
simple:
∂u ∂u
+u =0 (5.38)
∂t ∂x
La condición necesaria para que se cumpla la condición anterior es
1 ∂pA
+ C (2) = 0 (5.39)
ρA ∂x
Esta condición resulta útil cuando no se conoce con certeza el comportamiento de alguna
1 (3) ∂u
L − L(1) + C (2) = u
2ρAa ∂x
(3) (1) ∂u (2)
L =L + 2ρAa u −C (5.40)
∂x
5.3. Generalización del Modelo 107
(1) (3) ∂u (2)
L =L − 2ρAa u −C
∂x
∂u
El término
∂x se evalúa con una derivada descentrada, de la misma forma que para los
operadores L asociados a ondas salientes.
1 (3)
L + L(1) + C (3) = 0
2
Y para el derecho:
y un cierto valor de presión p∞ , asociado a un ujo en el campo lejano o far eld. Este
igual a cero si ambas presiones son iguales. Este esquema resulta útil para modelizar ujos
inestacionarios que pasan a través de salidas o entradas donde las condiciones pueden variar
tiempo. (Poinsot y Lele, 1992) emplean la siguiente expresión L(1) en el contorno derecho
donde K = σ(1−M 2 )aj /L es una constante que cuantica la intensidad de la onda reejada.
Ella depende de M: el número de Mach máximo del ujo, σ: un parámetro denido por el
usuario aj : la velocidad del sonido en el contorno, y L: una longitud característica del ujo.
∂u
Velocidad constante Esta condición implica
∂t = 0 ∀t, al igual que la condición de
pared sólida. Por esta razón se emplea la Relación (5.36) para determinar el operador L
108 Capítulo 5. Condiciones de Contorno
saliente, con la diferencia de que se emplea un valor de u distinto de cero como condición
a imponer.
!
L(2) γ − 1 L(3) L(1) ρa2 ∂T
+ + = ; (5.44)
u 2 u + a (u − a) T ∂x
( )
ρa2 ∂T L(2) 2 L(1)
L(3) = (u + a) − − (u + a) (5.45)
T ∂x u γ−1 u−a
Y para el derecho:
( )
ρa2 ∂T L(2) 2 L(3)
L(1) = (u − a) − − (u − a) (5.46)
T ∂x u γ−1 u+a
γ − 1 (1)
L(2) + (L + L(3) ) = −γC (3) ; (5.47)
2
pudiendo despejarse cualquiera de los L(i) en función de otros dos. Luego, para el contorno
( )
(3) (1) (3) L(2) 2
L = −L − γC + (5.48)
u γ−1
y para el derecho: ( )
L(2) 2
L(1) = −L(3) − γC (3) + (5.49)
u γ−1
5.3. Generalización del Modelo 109
5.3.3. Entrada
u>0 para x = xL
u<0 para x = xR
En este caso es necesario denir al menos dos operadores L en función de dos condiciones
de contorno físicas:
Flujo Supersónico:
En este caso todas las ondas características son entrantes al dominio. Por esta razón
será necesario además especicar el valor de L(1) para xL y el de L(3) para xR . Si se quiere
∂ρ ∂u ∂p
= = =0
∂t ∂t ∂t
L(1) = ρAaC (2) − C (3) ; L(2) = a2 C (1) − C (3) ; L(3) = −ρAaC (2) − C (3) ; (5.50)
∂ρ ∂u ∂p
Sino, pueden establecerse valores para
∂t , ∂t y ∂t y resolver el sistema de Ecuaciones (5.34)
para obtener los valores de los L(i) .
Flujo Subsónico:
Existen dos opciones posibles para establecer los valores de los dos operadores asociados
Establecer que la onda asociada a L(2) no se reeja, lo que equivale a establecer que
onda entrante en función de alguna condición física que se desee imponer (presión
(L
(3) para xL y L(1) para xR ), y con el valor obtenido establecer un valor para L(2)
en función de otra condición física que desee imponerse (caudal másico constante,
L(3) y L(1)
de donde
Las otras opciones posibles para denir L(2) necesitan de alguna condición sobre la onda
1 h i
L(2) = 2ρAa2 C (2) + 2a2 uC (1) − (u − a) L(1) − (u + a) L(3) (5.53)
2u
5.4. Extensión al modelo MHD 111
valor de L
(2) de la Ec. (5.35a), obteniéndose:
γ − 1 (3)
L(2) = − (L + L(1) ) − C (3) γ(γ − 1) (5.54)
2
valor de L
(2) de la Ec. (5.35b):
" #
(2) ρa2 γ − 1 L(3) L(1)
L =u − ( + ) (5.55)
T 2 u+a u−a
∂U ∂F
+ + C0 = 0 para xL ≤ x ≤ xR
∂t ∂x
∂V ∂V
+A +C=0
∂t ∂x
donde la matriz jacobiana A para el sistema primitivo normalizado está dada por la Ec.
(3.141), y sus vectores propios derechos e izquierdos, dados por las Ecs. (3.139), (3.147) y
∂U ∂V ∂Ui
=P →P=
∂t ∂t ∂Vj
están dadas por la Relación (3.124). Luego, para el tratamiento de los términos fuente en
C = P−1 C0
112 Capítulo 5. Condiciones de Contorno
0
gx
gy
C=
gz
(5.56)
0
0
2
(γ − 1) κ ∂∂xT2 + Lrad + EH
De la misma manera que para el caso de las ecuaciones de Euler, aplicando la formula-
ción de Thompson para condiciones de contorno de no reexión, dadas por las Ecs. (5.8) y
h α β c S ∂w α β c S ∂v
s z s s y s αf cf ∂u αf ∂p
L(1) = (u − cf ) 2
+ 2
− 2
+ 2
2a ∂x 2a ∂x 2a ∂x 2a ρ ∂x
αs βz ∂Bz αs βy ∂By i
+ √ + √ (5.57a)
2a ρ ∂x 2a ρ ∂x
h β ∂w β ∂v βy S ∂Bz βz S ∂By i
y z
L(2) = (u − bx ) − + √ − √ (5.57b)
2 ∂x 2 ∂x 2 ρ ∂x 2 ρ ∂x
h α β c S ∂w α β c S ∂v αs cs ∂u αs ∂p
f z f f y f
L(3) = (u − cs ) − 2
− 2
− 2
+ 2
2a ∂x 2a ∂x 2a ∂x 2a ρ ∂x
αf βz ∂Bz αf βy ∂By i
− √ − √ (5.57c)
2a ρ ∂x 2 a ρ ∂x
(4) ∂ρ 1 ∂p
L =u − 2 (5.57d)
∂x a ∂x
h α β c S ∂w α β c S ∂v αs cs ∂u αs ∂p
f z f f y f
L(5) = (u + cs ) + + +
2a2 ∂x 2a2 ∂x 2a2 ∂x 2a2 ρ ∂x
αf βz ∂Bz αf βy i
− √ − √ ∂By (5.57e)
2a ρ ∂x 2a ρ ∂x
(6) βy ∂w βz ∂v βy S ∂Bz βz S ∂By
L = (u + bx ) − + + √ − √ (5.57f )
2 ∂x 2 ∂x 2 ρ ∂x 2 ρ ∂x
5.4. Extensión al modelo MHD 113
h α β c S ∂w α β c S ∂v
s z s s y s αf cf ∂u αf ∂p
L(7) = (u + cf ) − − + +
2a2 ∂x 2a2 ∂x 2a2 ∂x 2a2 ρ ∂x
αs βz ∂Bz αs βy ∂By i
+ √ + √ (5.57g)
2a ρ ∂x 2a ρ ∂x
En este caso la formulación es más extensa que para el caso de Euler, ya que se tienen
siete ondas y la expresión de los vectores propios resulta más compleja debido a la norma-
cf , bx y cs . Sin embargo, como esta metodología posee la gran ventaja de resolver cada
Haciendo el producto escalar entre el vector propio l(1) del sistema MHD, dado por la
Ec. (3.147), con la derivada temporal del vector de variables primitivas y con el vector de
de la onda 1. De la misma manera se obtienen expresiones análogas para las ondas restantes.
A continuación se muestran las relaciones asociadas a las siete ondas características del
sistema:
Onda 1
(1) ∂V (1) (1) αf cf ∂u αs βy cs S ∂v αs βz cs S ∂w
l · +L +l ·C = − + gx + + gy + + gz
∂t 2 a2 ∂t 2a2 ∂t 2a2 ∂t
∂2T
αf ∂p
+ 2 + (γ − 1)κ 2 + (γ − 1)Lrad + EH (γ − 1)
2a ρ ∂t ∂x
αs βz ∂Bz αs βy ∂By
+ √ + √ + L(1) = 0 (5.58)
2a ρ ∂t 2a ρ ∂t
Onda 2
(2) ∂V βy ∂w
l · + L(2) + l(2) · C = + gz −
∂t 2 ∂t
βz ∂v βy S ∂Bz βz S ∂By
+ gy + √ − √ + L(2) (5.59)
2 ∂t 2 ρ ∂t 2 ρ ∂t
114 Capítulo 5. Condiciones de Contorno
Onda 3
(3) ∂V (3) (3) αs cs ∂u αf βz cf S ∂w αf βy cf S ∂v
l · +L +l ·C = − 2 ( +gx )− + gz − + gy
∂t 2a ∂t 2a2 ∂t 2a2 ∂t
αf βy ∂By αf βz ∂Bz
− √ − √ +
2a ρ ∂t 2a ρ ∂t
∂2T
αs ∂p
+ (γ − 1)κ + (γ − 1)L rad + EH (γ − 1) + L(3) = 0 (5.60)
2a2 ρ ∂t ∂x2
Onda 4
∂V ∂ρ
l(4) · + L(4) + l(4) · C =
∂t ∂t
∂2T
1 ∂p
− 2 + (γ − 1)κ 2 + (γ − 1)Lrad + EH (γ − 1) + L(4) = 0 (5.61)
a ∂t ∂x
Onda 5
∂V αf βz cf S ∂w αf βy cf S ∂v αs cs ∂u
l(5) · +L(5) +l(5) ·C = + gz + + gy + + gx
∂t 2a2 ∂t 2a2 ∂t 2a2 ∂t
∂2T
αs ∂p
+ + (γ − 1)κ (γ − 1)Lrad + EH γ − 1)
2 a2 ρ ∂t ∂x2
αf βz ∂Bz αf βy ∂By
− √ − √ + L(5) = 0 (5.62)
2a ρ ∂t 2a ρ ∂t
Onda 6
(6) ∂V βy ∂w βz ∂v
l · + L(6) + l(6) · C = − + gz + + gy
∂t 2 ∂t 2 ∂t
βy S ∂Bz βz S ∂By
+ √ − √ + L(6) = 0 (5.63)
2 ρ ∂t 2 ρ ∂t
Onda 7
(7) ∂V (7) (7) αs βz cs S ∂w αs βy cs S ∂v αf cf ∂u
l · +L +l ·C = − + gz − + gy + 2 + gx
∂t 2a2 ∂t 2a2 ∂t 2a ∂t
∂2T
αf ∂p
+ 2 + (γ − 1)κ 2 (γ − 1)Lrad + EH (γ − 1)
2a ρ ∂t ∂x
αs βz ∂Bz αs βy ∂By
+ √ + √ + L(7) = 0 (5.64)
2a ρ ∂t 2a ρ ∂t
5.4. Extensión al modelo MHD 115
de los operadores L
Teniendo en cuenta la expresión de los operadores L (dada por la Ec. (5.57)) es posible
∂V
K =L (5.65)
∂x
donde K es una matriz cuyas las son los vectores propios izquierdos l(i) multiplicados
(i)
por su velocidad característica asociada λ . De esta manera, premultiplicando a ambos
∂V
miembros por la inversa de K puede obtenerse las derivadas parciales
∂x en función de los
operadores L:
Estas relaciones resultan útiles si se desea imponer condiciones de contorno que impliquen
Empleando las Relaciones (5.66) y las Identidades (3.144), es posible expresar la varia-
ción temporal de las variables primitivas V en función del vector de términos fuente y de
los operadores L:
∂ρ
+ αf ρ L(7) + L(1) + αs ρ L(5) + L(3) + L(4) + C (1) = 0 (5.67a)
∂t
∂u
+ αf cf L(7) − L(1) + αs cs L(5) − L(3) + C (2) = 0 (5.67b)
∂t
∂v
− αs βy cs L(7) − L(1) + βz L(6) − L(2) + αf βy cf L(5) − L(3) + C (3) = 0
∂t
(5.67c)
∂w
− αs βz cs L(7) − L(1) − βy L(6) − L(2) + αf βz cf L(5) − L(3) + C (4) = 0
∂t
(5.67d)
√
∂By ρ
+ (a αs βy u − αs βy bx cs − αf by cf ) L(1)
∂t u − cf
√ √
ρ (3) ρ
− (aαf βy u + αs by cs − αf βy bx cf ) L − (aαf βy u − αs by cs + αf βy bx cf ) L(5)
u − cs u + cs
√
ρ √
+ (aαs βy u + αs βy bx cs + αf by cf ) L(7) − βz ρ L(2) + L(6) = 0 (5.67e)
u + cf
√
∂Bz ρ
+ (aαs βz u + αs βz bx cs + αf bz cf ) L(7)
∂t u + cf
√ √
ρ ρ
− (aαf βz u − αs bz cs + αf βz bx cf ) L(5) − (aαf βz u + αs bz cs − αf βz bx cf ) L(3)
u + cs u − cs
√
ρ √
+ (aαs βz u − αs βz bx cs − αf bz cf ) L(1) + βy ρ L(2) + L(6) = 0 (5.67f )
u − cf
∂p
+ ρa2 αf L(7) + L(1) + ρa2 αs L(5) + L(3) + C (7) = 0 (5.67g)
∂t
De esta manera, imponiendo condiciones sobre estas relaciones pueden modelarse los mis-
mos tipos de condiciones de contorno que para el modelo de Euler. Como primera medida,
5.4.4. Salida
Cuando se trata de una salida en uno de los extremos del dominio existirán al menos 4
Las ondas asociadas a L(4) , L(5) , L(6) , L(7) para el extremo derecho x = xR
Las ondas asociadas a L(1) , L(2) , L(3) , L(4) para el extremo izquierdo x = xL
Las ondas restantes pueden ser salientes o entrantes dependiendo de su velocidad caracte-
La onda 3 será saliente al dominio si u > cs : en este caso se tratará de un ujo súper
magnetosónico lento.
La onda 2 será saliente si u > bx : en este caso se tratará de un ujo súper alfvénico.
La onda 1 será saliente si u > cf : en este caso se tratará de un ujo súper magneto-
sónico rápido.
La onda 5 será saliente al dominio si |u| > cs : en este caso se tratará de un ujo super
magnetosónico lento.
La onda 6 será saliente si |u| > bx : en este caso se tratará de un ujo super alfvénico.
La onda 7 será saliente si |u| > cf : en este caso se tratará de un ujo super magne-
tosónico rápido.
De esta forma, existen 3 operadores posibles sobre los que se puede imponer condiciones
de contorno particulares y tres potenciales casos súper o sub. Para establecer condiciones
sobre las variables gasdinámicas (la presión, densidad, temperatura y velocidad) decidimos
emplear el operador L(1) por dos razones. Primero porque, como demostraron los auto-
res (Roe y Balsara, 1996), cuando el sistema degenera al caso gasdinámico puro la onda
aparece, de acuerdo a la Ec. (3.153). Además, si se establecen condiciones para las ondas
ujos super alfvénicos. Las ondas de Alfvén en este caso sólo transportan información de las
En este caso todas las ondas son salientes al dominio, y no debe imponerse ninguna
condición. Todos los operadores L(i) deben establecerse de acuerdo a las Relaciones de
compatibilidad (5.57).
Cada una de las condiciones sobre variables gasdinámicas en el modelo de Euler para
de L
(7) de la Ec. (5.67g) se obtiene una expresión análoga a la del contorno de presión
αs L(5) αs L(3) 1 ∂p
L(1) = − − − L(7) − 2 (5.69)
αf αf a αf ρ ∂t
para gasdinámica (Ec. (5.40)), con la diferencia que en la suma de fuerzas aparecen fuerzas
de origen magnético. La relación que debe satisfacer L(7) en el extremo izquierdo para
αs cs L(5) αs cs L(3)
(7) (1) 1 ∂u (2)
L =L − + + u −C (5.70)
αf cf αf cf αf cf ∂x
αs cs L(5) αs cs L(3)
(1) (7) 1 ∂u (2)
L =L + − − u −C (5.71)
αf cf αf cf αf cf ∂x
que para el sistema gasdinámico, esta condición se emplea para imponer paredes sólidas
∂u ∂u
(estableciendo
∂t =0 y u = 0) o salidas de velocidad constante (estableciendo
∂t =0 y
5.4. Extensión al modelo MHD 119
rápidas análogas a las del modelo gasdinámico (dadas por la Ec. (5.36)). Para el contorno
αs cs (3) αs cs (5) 1 ∂u
L(7) = L(1) + L − L − C (2) + (5.72)
αf cf αf cf αf cf ∂t
y para el derecho:
(1) (1) αs cs (5) αs cs (3) 1 (2) ∂u
L =L + L − L + C + (5.73)
αf cf αf cf αf cf ∂t
Flujo subalfvénico
Si el ujo es subalfvénico, entonces al menos dos condiciones deben ser impuestas so-
bre las variables en el contorno. Una de ellas está asociada a cualquiera de las variables
nicas rápidas, como las Ecs. (5.68), (5.70) o (5.72). Por otro lado, existen dos componentes
co (By y Bz ) que pueden restringirse mediante condiciones de contorno. Sin embargo, las
análogas en la dirección z . De esta manera, si ocurre una variación temporal en sólo una de
dichas variables transversales asociadas a la dirección y (por ejemplo, v ), necesariamente
ocurrirá una variación temporal para la otra (en este caso By ), pero las variables asociadas a
la otra dirección no serán perturbadas. Esto se observa en el tubo de choque de (Brio y Wu,
120 Capítulo 5. Condiciones de Contorno
1988), donde si como condición inicial se establece una perturbación sobre By solamente,
deberse a una particularidad del sistema MHD unidimensional, debido a que sólo aparecen
a la onda magnetosónica lenta (Ec. (5.62)) permite imponer vínculos entre variables tanto
gasdinámicas como magnéticas. Sin embargo, sabiendo cómo es la degeneración del sistema
MHD al sistema gasdinámico, creemos que no es una buena idea imponer restricciones sobre
variables gasdinámicas a través de L(5) , ya que podrían imponerse condiciones que, en caso
de que el ujo sea localmente gasdinámico, sean no físicas, redundantes o contradictorias.
A pesar de que el sistema permite una cierta exibilidad, siempre debe vericarse que
todas las condiciones impuestas a través de operadores L(i) satisfagan las ecuaciones de
5.4.5. Entrada
En el caso de tener una entrada al menos cuatro ondas son entrantes al dominio, y
las ondas salientes deben ser determinadas de acuerdo a las Ecs. (5.57). Además, la onda
asociada al cuarto vector propio (y al transporte de la entropía del ujo) será entrante. De
esta manera, si la ecuación de la energía no tiene términos fuente, el ujo será isoentrópico y
el valor físicamente correcto para L(4) será cero. Sin embargo, de manera análoga al modelo
∂T
=
∂t" #
L(4)
− αf (γ − 1) L(7) + αs (γ − 1) L(5) − + αs (γ − 1) L(3) + αf (γ − 1) L(1) T (5.75)
ρ
5.4. Extensión al modelo MHD 121
∂ρu
= −αf ρ (u + cf ) L(7) − αs ρ (u + cs ) L(5) −
∂t
uL(4) − αs ρ (u − cs ) L(3) − αf ρ (u − cf ) L(1) (5.76)
Para las tres ondas restantes pueden existir tres patrones de onda posibles, los cuales
se describen a continuación:
rápida. Luego, todas las ondas características serán entrantes al dominio. Una vez que el
ujo alcanza esta condición, ninguna información física puede viajar desde el dominio aguas
arriba de la entrada. Esto requiere que todos los operadores L sean iguales a cero, lo que
∂U
implica que en el contorno existirá una solución estacionaria tal que
∂t = 0.
sale del dominio. Si la entrada es isoentrópica, entonces puede establecerse el caudal másico
de entrada a través de L(1) o L(7) usando la Ec. (5.76). Empleamos dichos operadores
En esta tesis optamos por establecer condiciones de no reexión para todas las ondas de
Alfvén y magnetosónicas lentas entrantes. Sin embargo, este tema sigue abierto a discusión,
nes de contorno similar para modelar problemas astrofísicos de viento solar utilizando un
poladas para todas las variables a 300 radios solares de la supercie del sol (lo que equivale
a modelar esta región como una salida súper magnetosónica rápida). Para los contornos
entrada submagnetosónica lenta, deben imponerse cinco restricciones. El autor empleó tres
de las cinco restricciones para imponer que la componente normal del campo magnético B
122 Capítulo 5. Condiciones de Contorno
∂Bx − ∇ × (V × B) = 0
∂t
(5.77)
∂B
x =0
∂t
condición igualando a cero la derivada temporal de los términos del rotor de la Ec. (5.77).
∂ (V × B)y = 0 → u ∂Bz + Bz ∂u − Bx ∂w = 0
∂t ∂t ∂t ∂t
(5.78)
∂ (V × B) = 0 → u ∂By + B ∂u − B ∂v = 0
∂t z ∂t y ∂t x ∂t
incógnitas son los operadores L. Por otro lado, estas relaciones de vínculo no permiten
exibilidad, ya que si el ujo se acelerara y fuera super magneosónico lento sería necesario
imponer más condiciones sobre los L entrantes, modicando así el sistema de ecuaciones
para el caso anterior. Esto puede traer inconvenientes desde el punto de vista numérico y de
la implementación. El autor propone tres opciones diferentes para las otras dos restricciones
1. Limitar el ujo de masa cuando éste exceda un valor máximo crítico, a través de
relaciones de la forma:
∂
ρu = (ρu)crit ; (ρu) = 0 (5.79)
∂t
∂ p
(5.80)
∂t ρk
Sin embargo, no queda claro cómo funciona este modelo cuando se establece la velocidad
y su derivada temporal iguales a cero. En ese caso debería cesar todo ujo de masa hacia
5.4. Extensión al modelo MHD 123
da temporal de las variables primitivas en función de los operadores L(i) y de los términos
También es necesario hacer una extrapolación por serie de Taylor para RL, mediante
las Ecs. (5.19). Como se explicó anteriormente, dada la complejidad de las expresiones de
los operadores L(i) y de los vectores propios, las matrices jacobianas asociadas a las Ecs.
(5.21) se calculan de manera aproximada con la metodología propuesta por (Knoll y Keyes,
2004), ya que las expresiones analíticas son muy difíciles de obtener y propensas a errores
de implementación.
Capítulo 6
Resultados
problema de acuerdo a:
ρ p et u x tai
ρ= ; p= ; et = ; u= ; x= ; t= ;
ρi ρi a2i ρi a2i ai L L
Utilizamos en todos los casos como condición inicial una interpolación lineal entre el
celdas más 2 en cada uno de los contornos, el número de iteraciones varió de acuerdo al tipo
125
126 Capítulo 6. Resultados
Para este caso sólo deben denirse las condiciones a la entrada para las tres variables
de estado. Las condiciones a la salida dependerán sólo de los valores en el dominio, por lo
tanto deberán obtenerse mediante una relación de compatibilidad (Poinsot y Lele, 1992).
Para este caso sólo puede emplearse la condición de entrada supersónica para la entrada,
y la de salida supersónica. Comparamos los resultados con los obtenidos para condiciones de
contorno extrapoladas convencionales. Este es el caso más sencillo, y como era de esperarse,
el modelo de BCs no tuvo inuencia en los resultados ya que las BCs extrapoladas son
En ambos casos con menos de 500 iteraciones la solución convergió de manera práctica
al estado estacionario, ya que el residuo de la densidad tomaba valores del orden de 1 · 10−8 .
Resulta llamativo que ambos esquemas de BCs exhiben casi exactamente la misma curva
solución analítica predice gradiente nulo en la salida para todas las variables, condición que
es satisfecha de forma exacta por ambos modelos. Además, como las condiciones iniciales
son de ujo supersónico en todo el dominio los cambios durante el transitorio no pueden
Por esta razón, para esta condición el modelo de condiciones de contorno no tiene
Para este caso las condiciones a la entrada son idénticas al caso anterior, con la diferencia
que se impone una presión a la salida mucho mayor, de manera tal que aparezca una onda
6.1. Casos gasdinámicos 127
1
0.7
0.9
0.6
0.8
0.7 0.5
p
ρ
0.6
0.4
0.5
0.3
0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
0
1.9 10
1.85
1.8 −5
10
1.75
1.7 −10
10
eρ
u
1.65
1.6
−15
10
1.55
1.5
−20
10
0 0.2 0.4 0.6 0.8 1 200 400 600 800 1000 1200 1400
x n
Exacta Caracteristica Extrap Caracteristica Extrap
de choque en x = 5 ft
Mi = 1,5; pi = 2000 lb/ft2 ; ρi = 2,241 · 10− 3 slug/ft3 ; ai = 1118 ft/s;
en la entrada
pe = 4930 lb/ft2 ;
en la salida
128 Capítulo 6. Resultados
En todos los casos corrimos 2500 iteraciones con CF L = 1,8, los resultados obtenidos
2 1.8
1.8 1.6
1.4
1.6
1.2
1.4
1
ρ
p
1.2
0.8
1 0.6
0.8 0.4
0.2
0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
−2
2 10
1.8
1.6
−4
10
1.4
1.2
eρ
u
1 −6
10
0.8
0.6
−8
0.4 10
0 0.2 0.4 0.6 0.8 1 500 1000 1500 2000 2500
x n
p=cte force free et = cte p∞ no ref extrap
En este caso la elección de las condiciones de contorno juega un rol más importante en
directamente de cuánto sea el valor de la presión a la salida, el valor inicial de las variables
(5.40)) produjo, para CF L ≤ 1,8, menores oscilaciones durante el transitorio que los demás
esquemas, sin variar el valor de la presión en el contorno. Sin embargo, para valores de CF L
mayores a 1.8 el algoritmo converge a una solución con menor presión a la salida, con la onda
6.1. Casos gasdinámicos 129
de choque ubicada más hacia la derecha, que también es físicamente correcta. Este tipo de
qué ocurrirá con la presión; sin embargo debe tenerse en cuenta que la solución dependerá
buena, aunque la solución fue más fuertemente dependiente del valor de CF L. Para valores
de CF L ≥ 2,0 el algoritmo no llega a la convergencia, y para 1,4 < CF L < 2,0 converge a
soluciones cuya presión a la salida es mayor a la inicial. Sin embargo, si para las primeras
300 iteraciones se emplea otro esquema que preserve la presión de salida para luego cambiar
en el campo lejano (Ec. (5.43)) también produjo muy buenos resultados en cuanto a su
Es importante notar que todos los esquemas muestran curvas de convergencia muy si-
milares para las primeras 300 iteraciones. Esto se debe a que como las condiciones iniciales
dieren mucho de la solución nal, debido a que son simples interpolaciones lineales de las
condiciones en el contorno. Por esta razón, el estado del sistema evoluciona rápidamente
dentro del dominio, formándose una onda de choque que se mueve de posición, hasta esta-
bilizarse en su posición nal. Este estado transitorio demora en enviar su información física
hasta los contornos, por lo que el modelo de BCs no inuye en gran medida durante este
transitorio.
Para este caso la tobera se comporta como un difusor, y conociendo las condiciones a
la entrada sólo es posible una solución física, dada por las ecuaciones de ujo isoentrópico
Este caso fue el que más dicultades presentó para obtener una solución que converja,
130 Capítulo 6. Resultados
anterior, existen diferentes combinaciones posibles para dar valores a los operadores L(2)
y L(3) . Sin embargo, como el ujo que ingresa al dominio es isoentrópico la única opción
posible para L(2) es la de no reexión. Si se elige otra opción, el algoritmo produce resultados
no físicos. Si imponemos condiciones de presión constante o fuerza nula en el contorno el
el tiempo a través del operador L(3) (Ec. (5.52)) obtuvimos buenos resultados, con diferentes
propiedades de convergencia dependiendo del esquema de condiciones de contorno utilizado
para la salida. Los resultados obtenidos con dichas condiciones se muestran en la Figura 6.3.
También probamos establecer una entrada de velocidad constante, sin embargo no logramos
Nuevamente la condición de fuerza nula resultó ser la que mejores propiedades de con-
presión en el contorno en caso de que se utilice una condición inicial cuyo valor no sea
correcto.
de manera satisfactoria, siendo ésta última más propensa a inestabilidades para valores
necesario denir todas las variables a la entrada, como si se tratara de ujo supersónico. En
este caso es posible gracias a que se conoce la solución analítica del problema. Sin embargo,
(Yee, 1981), a menos que se conozca la solución exacta en el contorno (como en este caso).
Todos los esquemas probados muestran a grandes rasgos curvas de convergencia simila-
res para las primeras 800 iteraciones; a partir de ese momento algunos esquemas muestran
que el sistema evoluciona rápidamente dentro del dominio debido a que las condiciones
iniciales son muy diferentes a la solución analítica, como en el caso anterior. En el caso del
dentro del dominio y por la condición de presión impuesta a la salida. Por otro lado, para
1.3 1.05
1
1.25
0.95
1.2
0.9
1.15
0.85
p
ρ
1.1
0.8
1.05 0.75
1 0.7
0.95 0.65
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
−2
10
0.8
0.7 −4
10
0.6
u
eρ
−6
0.5 10
0.4
−8
10
0 0.2 0.4 0.6 0.8 1 500 1000 1500 2000 2500
x n
p=cte force free T = cte p∞ no ref extrap
Esto tiene una gran inuencia en el proceso de convergencia: cuando la solución cambia
dentro del dominio, hasta que alcanza una condición de pseudo equilibrio en la entrada.
que una onda viaje desde la salida hacia el interior del dominio, provocando cambios en la
solución. Una vez que dicha onda llega a la entrada, las variables en la misma evolucionan
para satisfacer la condición de caudal másico constante, y una onda se propaga hacia el
interior. Esta nueva onda produce cambios en el dominio hasta llegar a la salida, haciendo
que el proceso se repita. Esto ocurre muchas veces hasta que se logra la convergencia.
132 Capítulo 6. Resultados
un arco de la corona solar, propuesto por (Cargill y Priest, 1980). Dicho modelo consiste de
un arco circular con hidrógeno monoatómico, en cuyas bases la presión es igual a la presión
arco y despreciando su curvatura, el problema se idealiza con las Ecuaciones de Euler cuasi
se esquematiza el modelo.
Figura 6.4: Esquema del modelo del arco de (Cargill y Priest, 1980)
∂A πs πs T
S = 0, p − ρAg cos , −ρAug cos (6.2)
∂x 2L 2L
o isotérmico) y ujo estacionario, lo que les permitó obtener una relación analítica a través
a2 πs a2 ∂A
∂u
u− = −g cos + (6.3)
u ∂s 2L A ∂s
6.1. Casos gasdinámicos 133
variables.
presión. La discretización del dominio fue de 400 celdas, más dos en cada uno de los con-
tornos. Para este modelo tanto la entrada como la salida serán subsónicas, por lo tanto
empleamos la condición de caudal másico constante para la entrada y probamos las condi-
del arco, ya que la gravedad solar produce un gradiente en la densidad, cuyo mímino se
la descendiente.
Elegimos para simular el caso de un arco de longitud L = 140 · 106 m, cuya densidad y
h πs i
2
A(x) = A0 1 + (k − 1) sin ; L = 140 · 106 m; k=5
2L
del arco A(s = L). Teniendo en cuenta que la solución adimensional resulta independiente
entrada vale u0 = 0,27a. Para estas condiciones, la solución analítica predice ujo subsónico
en ambas bases (entrante a la izquierda y saliente a la derecha). Si la velocidad en la entrada
el ujo reduzca su velocidad (como en un difusor subsónico); y por otro lado la gravedad
variable, cuyo efecto decae con la altura, posibilitando que el ujo se acelere (la aceleración
soluciones posibles pueden tener varios extremos locales. De esta manera, si la contracción
segunda mitad del arco L < s < 2L. Para satisfacer la condición de presión en la cromósfera
para s = 2L necesariamente debe aparecer una onda de choque en dicho segmento del arco.
Sin embargo, de acuerdo a la solución analítica dicha onda puede aparecer en cualquier parte
134 Capítulo 6. Resultados
de la semilongitud del arco. Esto presenta un problema desde el punto de vista numérico: es
imposible que el algoritmo converja a una solución estacionaria si existen innitas soluciones
posibles entre L < s < 2L; cualquier solución obtenida cambiaría hacia otra solución con
se ve claramente que la distribución de velocidades a lo largo del arco posee tres máximos
locales
Empleamos la condición de caudal másico constante a la entrada, ya que dio buenos re-
sultados para el caso del difusor subsónico, y la solución también implica ujo isoentrópico.
presión impuesta en el campo lejano (Ec. (5.43)) y fuerza nula en el contorno (Ec. (5.40)).
Los resultados obtenidos se muestran en la Figura 6.5, donde se presentan los perles
Densidad Velocidad
2 0.8
1.8
0.7
1.6
0.6
1.4
0.5
1.2
1 0.4
0.8
0.3
0.6
0.2
0.4
0.1
0.2
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
regiones donde el ujo es supersónico, que eventualmente pueden cambiar las condiciones
en la salida (dependiendo del tipo de esquema de BCs utilizado), haciendo que el sistema
6.1. Casos gasdinámicos 135
converja a una solución no física (con una presión en la cromósfera diferente de la real).
El sistema también resulta muy sensible a las condiciones iniciales, pudiendo obtenerse
presiones negativas en el caso de que el paso de tiempo sea sucientemente grande. Esto
velocidad que satisfacen las ecuaciones de conservación. En la Figura 6.6 se muestran las
Residuo de la densidad
−2
10
−3
10
−4
10
−5
10
E xtrap .
p=cons
T=cons
−6
10 v=cons
fza nula
no ref
p∞
−7
10
0 1000 2000 3000 4000 5000 6000
n
Figura 6.6: Curvas de convergencia para distintos esquemas de BCs para el arco
coronal de (Cargill y Priest, 1980)
Para evaluar la inuencia de los términos fuente relativa a la de los efectos convectivos,
calculamos la relación entre las componentes del vector de términos fuente Ci0 y las com-
ponentes del vector de ujos numérico Fi para las ecuaciones de cantidad de movimiento
C20 (x = 0) C30 (x = 0)
= 2,50 = 1,093 (6.4)
F2 (x = 0) F3 (x = 0)
De esta manera, los efectos gravitatorios tienen un rol importante en el modelo de BCs,
ya que sus efectos son mayores a los de la convección, y a que los efectos de variación de
constante en el contorno, aunque vimos que la convergencia fue bastante errática, y para
casos más severos (como para k = 20 y L = 175 · 106 m) una vez alcanzada la solución
analítica, ésta no era estable y después de un tiempo el algoritmo convergía a otra solución.
de convergencia, pero ocurrió que resultó muy sensible a las condiciones iniciales, debido a
que el ujo se volvía supersónico durante el transitorio, provocando que la salida se volviera
también supersónica y cambiara el valor en la presión. Una estrategia que empleamos para
subsanar este problema fue utilizar la condición de salida con presión constante para las
primeras 3000 iteraciones, y luego cambiar a las condiciones de no reexión. Sin embargo,
este problema quedó resuelto con el esquema de presión impuesta en el campo lejano, que
resultó ser una buena solución de compromiso entre ambos esquemas. A pesar de ello,
la tasa de convergencia no es tan buena como hubiéramos esperado. Podría ser mejorada
quizás ajustando el coeciente σ en la Ec. (5.43) a un valor mas apropiado para este tipo
∂pA dA
+ ρAg − p =0
∂x dx
El algoritmo converge a una solución con expansión supersónica, que verica la relación
anterior. Finalmente, debemos notar que el esquema de extrapolación para BCs posee
una curva de convergencia de mayor velocidad y con menor cantidad de oscilaciones. Sin
embargo, como en los casos anteriores, esta solución requirió denir todas las variables en
cuando tienen que analizarse casos cuya solución analítica no se conoce en el contorno. Si
bien en la Figura 6.5 se observa que esta solución es la que más se acerca a la analítica, lo que
ocurre es que se mantiene oscilando. Sin embargo, las condiciones extrapoladas probaron
Para validar los efectos de los términos convectivos, de conducción de calor y de fuentes
viscosidad nula
El sistema de ecuaciones del modelo de Euler para el caso 1D estacionario para un dominio
d dρ
dx (ρu) =0→ dx u + ρ du
dx =0
d dρ
2
dx (ρu + p) =0→ ρu du
dx + dx RT + ρR dT
dx =0 (6.6)
d2 T d2 T
d
dx [u (E + p)] = κ dx2
+ Q˙cal → ρuCp dT
dx + ρu du
2
dx =κ dx2
+ Q˙cal
ρu = ρ0 u0
(6.7)
ρu2 + ρRT = ρ0 u20 + ρ0 RT0
Mediante manipulación algebraica de las Leyes de conservación (6.6) y usando las Relaciones
RT du dT R dT du
u− +R = 0 → RT = (6.8)
u dx dx u −u
dx dx
dT R dT d2 T
ρuCp + ρu2 RT = κ 2 + Q˙cal (6.9)
dx u −u
dx dx
luego
R2 T d2 T
dT
ρu Cv + = κ + Q˙cal (6.10)
RT − u2 dx dx2
138 Capítulo 6. Resultados
miento (6.7):
ρ0 2 ρ0
u2 = u + RT0 − RT (6.11)
ρ 0 ρ
RT0
puede resolverse esta ecuación cuadrática para u, deniendo el término auxiliar Γ = u0 + u0
Γ 1p 2
u= ± Γ − 4RT (6.12)
2 2
integral (6.7), se obtiene una ecuación diferencial ordinaria no lineal de segundo orden que
" #
R2 T dT d2 T
ρu Cv + Γ2
√ = κ 2 + Q˙cal (6.13)
2RT − ∓ Γ2 Γ2 − 4RT dx dx
2
Esta ecuación tiene solución analítica con una integral impropia. Por esta razón resulta
más práctico resolverla mediante un esquema de diferencias nitas. Debe tenerse en cuenta
Γ2
T ≤
4R
Dicho problema puede resolverse mediante diferencias nitas discretizando la derivada pri-
Para que el esquema con derivada primera centrada sea estable es necesario agregar
u∆x
P e∆x = (6.16)
2α
ρuCp ∆x ρ0 u0 Cp ∆x
P e∆x = =
2κ 2κ
La viscosidad numérica necesaria para establizar el problema para P e∆x > 1 y que a su
Existen distintas opciones para elegir la función f (P e∆x ). En este caso se eligió:
1 1
f (P e∆x ) = − (6.17)
tanh(P e∆x ) P e∆x
y que con ella la solución aproximada converge con mayor precisión a la solución analítica
e
2 P eL T0 − TL 2 P ex
e (T0 − TL )
T (x) = − (6.18)
e2 P eL − 1 e2 P eL −1
donde
ρ0 u0 L Cp ρ0 u0 x Cp
P eL = ; P ex = (6.19)
κ κ
de ecuaciones no lineales
donde
R2 Ti
∂Fi h ρ0 u0
Ji,i−1 = = −1 − Cv +
∂Ti−1 2 κ + κnum Γ2 Γ √ 2
2RTi − ∓ Γ − 4RTi
2 2
∂Fi h ρ0 u0 R 2
Ji,i = =2+ { − ...
∂Ti 2 κ + κnum Γ2 Γ√ 2
2RTi − ∓ Γ − 4RTi
2 2
ΓR (6.22)
R2 Ti 2R ∓ √ 2
Γ − 4RTi
2 } (Ti+1 − Ti−1 )
2 Γ√ 2
Γ
2RTi − ∓ Γ − 4RTi
2 2
R2 Ti
∂Fi h ρ0 u0
Ji,i+1 = = −1 + Cv +
∂Ti+1 2 κ + κnum Γ2 Γ √ 2
2RTi − ∓ Γ − 4RTi
2 2
J · Y = −F
T(k+1) = T(k) + Y
tomando como condiciones iniciales una interpolación lineal entre los valores en los contor-
6.1. Casos gasdinámicos 141
nos.
ρ0 = 1,225 kg ; u0 = 100 ms ;
m3 T0 = 400 K; en la entrada
T = 500 K en la salida
L
Se observa que las curvas de temperatura de las soluciones por diferencias nitas coinciden
500 130
490
125
480
470
120
460
u[m/ s]
T [K ]
450 115
440
110
430
420
105
410
400 100
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
Incompr. diferencia centr. diferencia upwind.
Figura 6.7: Solución incompresible y soluciones por diferencias nitas de la Ec. (6.13)
y que existe una diferencia inferior al 2 % respecto a la curva que representa la solución
analítica de la ecuación de la advección difusión incompresible. Esto se debe a que para las
condiciones iniciales planteadas el número de Mach es menor a 0,3 en casi todo el dominio, lo
cual permitiría asumir una solución de ujo incompresible. Por otro lado, las distribuciones
solución asume u = cte y ρ = cte. Sin embargo, las soluciones asociadas a ambos esquemas
lineal entre las variables en los contornos. Para llegar a convergencia fueron necesarias 8000
nitas centradas. Para calcular las matrices jacobianas de los términso fuente expresamos
la temperatura en función de las variables conservativas empleando las Ecs. (3.24) y (3.27).
(Poinsot y Lele, 1992) proponen que para modelar condiciones de contorno para las
∂q ∂2T
=0→ =0 (6.24)
∂x ∂x2
Sin embargo, en el presente problema para obtener convergencia fue necesario imponer
imponiendo solamente una condición sobre alguna de las variables (presión, velocidad o
temperatura) más la condición dada por (6.24) a la salida, el algoritmo oscila sin llegar
a un estado estacionario. Creemos que esto puede deberse a la gran importancia relativa
del término de conducción de calor: el sistema deja de comportarse como hiperbólico, para
puesta por el operador L(2) para ondas salientes (5.33b) implica que la entropía del sistema
se conserva, lo cual no es físicamente correcto para la salida derecha, debido al ujo de calor
existente.
Por estas razones, obtuvimos los mejores resultados cuando impusimos los valores de
salida. Otra opción que probamos fue imponer la velocidad y el gradiente de temperatura
a la salida, pero en este caso el sistema converge a una solución con menor temperatura.
una condición de ujo con velocidad y temperatura constantes a lo largo del dominio. Lo
mismo ocurre con las condición de fuerza nula en el contorno. La Figura 6.8 muestra las
probados (dados en la Tabla 6.1), y la Figura 6.9 muestra los respectivos perles de presión
lograr la convergencia fue estableciendo las tres variables en cada uno de los contornos, es
decir que fue necesario imponer dos condiciones adicionales a las físicamente necesarias.
En la Tabla 6.1 se especican las condiciones de contorno que utilizamos en cada uno
500 130
125
480
120
460
115
u[m/s]
T[K]
440
110
420
105
400
100
380 95
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x
extrap 1 2 3 4 5 exact
5
x 10
1.44 1.25
1.42 1.2
1.4 1.15
ρ[kg/m3 ]
1.38 1.1
p [Pa]
1.36 1.05
1.34 1
1.32 0.95
1.3 0.9
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
extrap 1 2 3 4 5 exact
Como se puede ver en la Figura 6.10 en general se llega a la convergencia de forma más
efectiva cuando se establecen las tres variables a la entrada. Cuando se establece sólo la
velocidad, o el caudal másico la solución oscila cuando se encuentra a una solución próxima
cidad (casos 1, 2 y 5) se calcula L(3) para velocidad constante, de acuerdo a la Ec. (5.36),
144 Capítulo 6. Resultados
1 ρ , v T ,v
2 ρ , v, p T ,v
3 ρ , v, p T , ρv
4 ρ , v, p ρ,v
5 v T ,v
−5
10
−6
10
eρ
−7 extrap
10 1
2
3
4
5
−8
10
0 1 2 3 4 5 6 7 8
n 4
x 10
Figura 6.10: Curvas de convergencia para los diferentes modelos de BCs dados en
6.1.
y luego se reemplaza este valor en la expresión de L(2) para temperatura constante, dada
3) es necesario resolver un sistema de dos Ecuaciones (la (5.47) y la (5.51)) cuyas incógnitas
a(γ−1)L(3)
L(2) = γu−a
(3) (6.25)
L(1) = − (γu+a)L
γu−a
6.1. Casos gasdinámicos 145
ρ0 u0 CP L
∂T ρ0 u0 Cp e κ (T0 − TL )
= ρ0 u0 Cp L (6.26)
∂x x=L
κ e κ −1
Luego puede calcularse el valor del operador L(1) en función de L(2) y L(3) mediante la
Relación (5.46). Resultó más práctico imponer la condición de velocidad constante mediante
en un sistema de ecuaciones que tiene como soluciones expresiones complejas para L(1) y
L(2) . Por esta razón, y porque la otra versión de implementación de condiciones basada en
esta versión.
El tubo de choque de Sod es un caso de prueba clásico para validar esquemas numéricos
para ujo compresible. Consiste en un problema de Riemann para gases perfectos con
onda de choque que viaja hacia la derecha seguida por una discontinuidad de contacto. La
onda de choque produce un aumento de presión y temperatura en el uido que atraviesa. Por
otro lado, un abanico de expansión viaja hacia la izquierda, que produce una disminución
En este caso decidimos probar el caso con dos tipos de condiciones de contorno:
derecha
146 Capítulo 6. Resultados
Pared sólida
Para el caso de la pared sólida la solución exacta se obtiene mediante la teoría de las
características y no presenta mayores problemas desde el punto de vista numérico. Las ondas
genuinamente no lineales rebotan contra las paredes generando ondas del mismo tipo (la
onda de choque rebota como onda de choque, y la de expansión rebota como expansión). La
solución analítica puede encontrarse en el trabajo nal de grado de (Acosta Lusa y Tamagno,
2004), así como soluciones numéricas para ambos tipos de condiciones de contorno. Para
modelar una pared sólida con el esquema de BCs extrapoladas, se extrapolan todas las
es reejada:
U(2) = −U(2)
m
m+1
(6.27)
U(l) = U(l) para l 6= 2
m+1 m
choque como el abanico de expansión son correctamente reejados por ambos esquemas.
En la Figura 6.11 se comparan los contornos de densidad obtenidos con ambos esquemas,
y en la Figura 6.12 los contornos de velocidad para 600 pasos de tiempo. Sin embargo,
0.6 0.6
0.8 0.8
0.5 0.5
0.3 0.3
0.4 0.4
0.2 0.2
0.3 0.3
0.1 0.1
0.2 0.2
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
X X
0.8
0.4
0.5 0.5
0.6
0.2
0.4 0.4 0.5
t
t
0
0.4
0.3 0.3
−0.2 0.3
0.2 0.2
−0.4 0.2
0.1 0.1
−0.6 0.1
0 0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
X X
contorno (para garantizar ujo igual a cero) cuando la onda de choque se reeja (de acuerdo
a la Ec. (6.27)). Los valores negativos en la escala del gráco se dejaron para hacer énfasis
en este punto. Por otro lado, el esquema basado en características produce variaciones en
la presión y la densidad para mantener la velocidad igual a cer0. Esto no ocurre para el
comportó de manera muy satisfactoria, conservando la energía del sistema. Cuando la onda
comparamos la evolución de la energía total del sistema ρet , integrada en la longitud del
tubo,
Z L
Et = ρet dx (6.28)
0
La interacción del abanico de expansión con la pared también genera cambios de energía
en el sistema, en este caso disminuyéndola de manera progresiva a medida que rebotan las
ondas del abanico. La cantidad de energía perdida en este caso es también del orden del
148 Capítulo 6. Resultados
279
Extrap.
278.5 Char.
278
277.5
Et
277
276.5
276
275.5
275
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
t
Figura 6.13: Evolución de la energía total del sistema con BC de pared cerrada para
ambos esquemas
1 %. Otra diferencia que notamos entre ambos esquemas es un incremento muy pequeño
en el tamaño del paso de tiempo cuando la onda de choque llega a la pared. El paso de
tiempo para el esquema extrapolado es casi imperceptiblemente mayor (menos del 0,1 %),
el efecto acumulado en el tiempo de cálculo a lo largo de varios pasos da una diferencia de
Para el caso del conducto abierto a la atmósfera la física del problema es más compleja.
Como primera medida deben distinguirse dos casos posibles: cuando el ujo detrás de la
onda de choque es supersónico, o cuando éste es subsónico. En el caso de que el ujo sea
subsónico (como en el test de Sod) cuando una onda característica llega a un extremo
abierto que conecta a una región de presión constante, dicha onda debe reejarse hacia
adentro como una onda de la familia opuesta. De esta manera se mantiene la condición de
presión constante e igual a la presión de salida. Esto implica que si la onda incidente es
de expansión, ésta debe reejarse como onda de compresión, y viceversa. Este resultado se
extiende para choques y expansiones: si una onda de choque llega a un extremo abierto,
debe reejarse como un abanico de expansión. Una discusión completa puede encontrarse
6.1. Casos gasdinámicos 149
Sin embargo, de acuerdo a fotografías Schlieren, lo que ocurre en realidad es que cuando
la onda de choque llega al contorno derecho y sale al exterior se forma en la salida una
estructura de ujo tridimensional inestacionario. Esto hace que la presión aumente por un
1955), si la velocidad detrás del choque es tal que el ujo es subsónico, una vez que la
de expansión, que luego se propaga hacia adentro del tubo como una serie de ondas de
en la Figura 6.14.
Figura 6.14: Reexión de una onda de choque como ondas de expansión discretas,
tomado de (Rudinger, 1955)
El autor además explica que mientras mayor sea la intensidad de la onda de choque,
menor serán los intervalos de tiempo entre las ondas de expansión reejadas. Por esta
razón, para choques de mucha intensidad puede modelarse la reexión de los mismos como
Por otro lado, si el ujo detrás de la onda de choque es supersónico, las condiciones
de contorno analíticas para un estado estacionario exigen que la velocidad en la salida sea
se comportaba mejor, y así poder determinar cuál era sucientemente simple y preciso
para emplearse en casos más complejos con salidas inestacionarias a extremos abiertos.
M = 0,9, de manera que alguna condición debe imponerse sobre la onda que se reeja
hacia adentro. Naturalmente, una condición de no reexión resulta inapropiada. Una opción
posible sería una salida de presión constante. Sin embargo, este tipo de condición tiene dos
todo de con cuánta precisión resuelve el ujo numérico de Harten-Yee a la onda de choque
dentro del dominio. Si el choque es capturado en tres celdas, se reejarán tres ondas desde
el contorno. Segundo: si el ujo se vuelve sónico en la salida, la presión en los nodos del
hacen que el ujo a la salida vuelva a ser subsónico, el esquema de condiciones de contorno
mantendrá constante este nuevo de valor de presión en lugar de la presión física en la salida,
∂p
dada por las condiciones iniciales. Esto es porque el esquema numérico impone
∂t = 0, en
valor de la presión en las celdas de la salida una vez que el ujo se volviera subsónico
nuevamente.
Otra opción que evaluamos fue emplear el esquema de la presión impuesta en el far
eld p∞ . Este esquema, para el valor de σ = 0,25 utilizado por (Poinsot y Lele, 1992), se
del dominio prácticamente no inuye en la presión aguas arriba y muestra una interacción
Para obtener resultados similares a los obtenidos con una salida de presión constante
fue necesario usar un parámetro σ 100 veces mayor, es decir σ = 25. Además, emplear en la
Expresión (5.43) el valor del número de Mach máximo del ujo produjo menores reexiones
y algunas oscilaciones. Optamos por emplear el valor del número de Mach en el contorno
en su lugar, y obtuvimos mejores resultados. En las Figuras 6.15 y 6.16 mostramos una
salida de presión constante y para una salida con presión impuesta en el far eld, con tres
valores de σ.
A pesar de que el abanico de expansión reejado no se observa con tanta claridad como
t
0.4 0.5 0.4 0.5
0.4 0.4
0.2 0.3 0.2 0.3
0.2 0.2
0 0
0 0.5 1 0 0.5 1
X X
0.9 0.9
0.8
0.8
0.8 0.8
0.7 0.6 0.7
0.6
0.6 0.6
t
Figura 6.15: Contornos de densidad obtenidos para distintos tipos de salidas con
presión impuesta.
que una onda se reeja hacia adentro del dominio, pero la intensidad de dicha onda es
signicativamente menor que la que se observa para el caso de salida de presión constante.
En la Figura 6.17 se observan los resultados obtenidos para el esquema de fuerza nula.
RL
Debe notarse además que la energía total del sistema
0 ρet dx no se ve signicativa-
mente afectada por el esquema de condiciones de contorno utilizado. Las curvas asociadas
El siguiente paso fue evaluar lo que ocurre después de que el ujo se vuelve supersónico
en la región vecina a la salida. En este caso probamos con dos condiciones posibles para la
152 Capítulo 6. Resultados
t
0.4 0.4
0.4
0.5
0.2 0.2 0.2
0 0 0 0
0 0.5 1 0 0.5 1
X X
0.8
0.8 0.8 1
0.6 0.8
0.6 0.6
0.6
t
0 0 0 0
0 0.5 1 0 0.5 1
X X
Figura 6.16: Contornos de velocidad obtenidos para distintos tipos de salidas con
presión impuesta.
salida:
t
t
0.2
0 0 0
0 0.5 1 0 0.2 0.4 0.6 0.8 1
X X
280
260
240
220
Et
200
180
p=cons.
σ = 0.25
160 σ = 2.5
σ = 25
force free
140
120
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t
Figura 6.18: Evolución de la energía total del sistema para el caso de extremo abierto
que ocurran aguas arriba, el esquema de BCs de salida de presión constante preservaría
154 Capítulo 6. Resultados
el valor de la presión actual. Dicho valor es distinto del valor de la presión en la salida,
determinado por las condiciones iniciales. Para salvar este inconveniente es necesario res-
tablecer la presión en las celdas del contorno, de la misma manera que para el esquema de
sónicas (u =cte o M =cte) junto con el esquema de presión constante para salidas subsóni-
cas produjo inestabilidades, como se muestra en la Figura 6.19, que llevaron a un colapso
Contornos de M=cte para salida sónica Contornos de M=cte para salida sónica
de M=cte y subsónica de p=cte de M=cte y subsónica de p=cte
0.32 0.32
1.5 3
2.5
0.3 0.3
1 2
t
t
0.5 1
0.26 0.26 0.5
0 0
0.24 0.24
0.8 0.85 0.9 0.95 1 0.8 0.85 0.9 0.95 1
X X
Figura 6.19: Comparación de esquemas de salida sónica con p =cte para ondas
subsónicas
Esto se debe a que cuando el ujo llega a condición sónica en la salida la presión
aumenta. Pero cuando por la propia evolución del sistema el ujo se vuelve subsónico
en la velocidad del sonido, que a su vez incrementa el número de Mach local. Este proceso
se repite varias veces, siendo más acentuado cuando se emplea la condición de velocidad
constante.
Este tipo de inestabilidad no fue tan severa cuando empleamos el esquema de presión
de tiempo cuando la salida se volvía sónica. Sin embargo, observamos que aguas arriba de
la salida el ujo se vuelve siempre supersónico, sin importar el esquema utilizado. Esto nos
lleva a pensar que la solución físicamente correcta es una salida supersónica inestacionaria,
Por último comparamos los resultados obtenidos con el esquema extrapolado. Lo im-
través de la energía total y extrapola los valores de las otras variables conservadas.
La Figura 6.20 muestra los contornos de Mach constante obtenidos para la condición
Contornos de M=cte para salida sónica Contornos de M=cte para salida sónica
de v=cte y subsónica de p impuesta de M=cte y subsónica de p impuesta
0.7
1.2 0.7
0.6 1
1 0.6
0.5 0.5 0.8
0.8
0.4 0.4 0.6
t
t
0.6
0.3 0.3
0.4 0.4
0.2 0.2
0.1 0.2 0.1 0.2
0 0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
Figura 6.20: Comparación de esquemas de salida sónica con p∞ impuesta para ondas
subsónicas
que viaja hacia la derecha reejado desde la pared izquierda con el abanico de expansión
cerca de la salida. Por otro lado, el esquema extrapolado produjo pequeñas oscilaciones en
la salida, sobre todo luego de que la onda de choque llegara a la salida y cuando el abanico
t
0.2
1 0.4 1
0
0.5 0.5
0.2 −0.2
0 0 −0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
2 2 0.6
0.6
t
0.5 0.5 0
0.2
0 0 −0.2
0 0.2 0.4 0.6 0.8 1 0 0.5 1
X X
2 2
0.6 0.5
1.5 1.5
t
t
0.4
1 1
0
0.5 0.2 0.5
0 0 −0.5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
1.5 1.5
0.5
0.4
1 1
0 0
0 0.2 0.4 0.6 0.8 1 0 0.5 1
X X
(s − s0 )
EH = EH0 exp −
sH
158 Capítulo 6. Resultados
Para el caso general de un arco con una distribución de temperatura, sección y una fun-
T (s) = Tmax [1 − z a ]b para
sH
> 0,3
L
h i (6.30)
T (s) = Tmax [1 − z a ]b 1 + 0,5 · log10 L
(1 − z) z 5 para
sH
≤ 0,3
sH L
h(s)−h0
−
p(s) = p0 e λp (s)
T (s) h(s)
λp (s) = λ0 1+ qλ (L, sH )
106 Rsun
2L πs
h= sin
π 2L
a2
L0
a(L, sH ) = a0 + a1 sH
b 2
L0
b(L, sH ) = b0 + b1
sH c2
qλ (L, sH ) = c0 + c1 sLH0
Este modelo sencillo permite validar los términos fuente para arcos de diferentes longi-
A continuación se presentan los resultados para un arco de sección constante con semilongi-
tud L = 40 · 106 m. Para este caso además existen grácos de resultados de una simulación
se emplearon para comparar resultados. Los datos de dicha simulación se levantaron co-
mo puntos discretos de las guras del trabajo citado mediante un programa para tomar
Una vez obtenidas las soluciones hidrostáticas aproximadas para la presión y densidad,
se emplearon éstas como condición inicial para el código de FVM desarrollado. Con dichas
las empleadas en la Sección 6.1.3, corrimos 4000 pasos de tiempo con CF L = 0,9. En ese
a una solución hidrostática idéntica a la obtenida por (Aschwanden y Schrijver, 2002) con
Por otro lado,para la condición inicial planteada el código de volumenes nitos dio como
resultado velocidades casi nulas en todo el dominio (salvo pequeñas oscilaciones cerca de
Figura 6.25 se presentan los perles de densidad y temperatura obtenidos, comparados con
por radiación para las tres soluciones mencionadas. De esta forma, comparamos los términos
fuente calculados por el código FVM para dichas condiciones como los calculados para
casos se obtuvo una muy buena correlación. Para evaluar los ujos de calor por conducción
nitos, Ec. (4.37). Ambos esquemas funcionaron adecuadamente, existiendo diferencias del
orden del 1 % entre los resultados obtenidos con el método de volúmenes nitos y los con
5
x 10 11
12 10
Sol Aschwanden
Sol hidrost
10 Sol FVM
10
8 10
log10 ρ
T[K]
9
4 10
8
0 10
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
9 9
x 10 x 10
Por otro lado, evaluamos el efecto de considerar una limitación por saturación del ujo
de calor por conducción, de acuerdo a la Ec. (3.12). Dicha condición no fue tenida en cuenta
en la solución obtenida por los autores. Sin embargo, se observa que la limitación del ujo
de calor posee escasa inuencia en la cantidad de calor transportada. Creemos que esto
se debe a que en las zonas de las bases (donde existe mayor gradiente de temperatura) la
densidad es más elevada, lo que implica un valor del ujo de saturación qsat mayor.
Por otro lado, derivamos el ujo de calor por conducción de la solución aproximada
por la Ec.(6.30). Se observa que la solución aproximada para el ujo de calor (obtenida
diferente de la solución numérica obtenida por Aschwanden. Creemos que esto se debe a que
la solución aproximada es menos exacta en las regiones cercanas a las bases, y estas pequeñas
desarrollado muestra muy buena coincidencia con los valores obtenidos por la solución
−4 −4
x 10 x 10
0.5
0 L Asch
rad
Lrad Hid
−0.5 Lrad FVM
−1 2
Qcond [erg/cm3 s]
[erg/cm3 s]
−1.5
−2
Rad
L
−2.5 Cond. Cal. Asch. 1
Cond. Cal. An. hid
Cond. cal. ap. hid
−3 Cond. Cal. FDM
Cond. Cal. FVM
Cond. Cal. FVM Lim.
−3.5
−4 0
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x[cm] 9
x 10 x[cm] 9
x 10
Como primera medida, corrimos el test de Sod con el algoritmo MHD para comparar los
resultados con el algoritmo que resuelve las ecuaciones de Euler, con el objetivo de vericar
En este caso el dominio computacional fue discretizado en 200 celdas y ambos contornos
con dos celdas fantasmas. Corrimos 2000 pasos de tiempo con CF L = 0,5.
Al igual que para el test de Sod analizado con el modelo gasdinámico, resolvimos el
problema con dos tipos de condiciones de contorno: ambos extremos con paredes cerradas,
Pared cerrada
Comparamos el problema del test de Sod con condiciones de contorno de pared cerrada
MHD usamos condiciones de no reexión para las ondas magnetosónicas lentas y de Alfvén.
a ondas magnetosónicas rápidas. En las Figuras 6.27 y 6.28 se muestran los resultados para
2 0.8 2
0.6
1.5 0.6 1.5
t
t
0.4
1 1
0.4
0.5 0.5 0.2
0.2
0 0 0
0 0.5 1 0 0.5 1
X X
2.5 2.5
0.5 0.5
2 2
1.5 0 1.5
t
0
1 1
−0.5 −0.5
0.5 0.5
0 0
0 0.5 1 0 0.5 1
X X
ondas en ambos casos. Tanto la onda de choque que viaja hacia la derecha como el aba-
nico de expansión son correctamente reejadas en ambas paredes sólidas. Estos resultados
6.3. Casos Magnetohidrodinámicos 163
características.
Extremo abierto
Por otro lado, repetimos la comparación entre el modelo MHD y gasdinámico para
dos variantes de condiciones de contorno para extremos abiertos: con presión constante
∂p
implementada a través de
∂t =0 y de presión impuesta en el far eld. Hacemos esto sólo
a nes de comparación, más allá de una precisa descripción de la física del problema, como
vimos en la sección anterior. En las Figuras 6.29 y 6.30 mostramos la comparación de los
∂p
contornos de densidad y velocidad obtenidos para el caso de
∂t = 0.
0.9
2.5 2.5 0.8
0.8
2 0.7 2
0.6
0.6
1.5 1.5
t
0.5 0.4
1 0.4 1
Se observa que en todos los casos los resultados comparados entre ambos modelos
coinciden, y se observa que las ondas poseen las mismas velocidades de propagación. La
∂p
diferencia entre el contorno con p∞ impuesta y con
∂t =0 se debe al fenómeno que ocurre
cuando el ujo se vuelve supersónico a la salida, que produce una evolución temporal del
2.5 2.5
1 1
2 2
0.5 0.5
1.5 1.5
t
t
1 0 1
0
0.5 0.5
−0.5 −0.5
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
2 2 0.6
0.6
t
1.5 1.5
0.4
1 0.4
1
0.2
0.5 0.5
0.2
0 0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
rápidas.
6.3. Casos Magnetohidrodinámicos 165
t
1.5 1.5
0.4 0.4
1 1
0.2 0.2
0.5 0.5
0 0
0 −0.2 0
0 0.5 1 0 0.2 0.4 0.6 0.8 1
X X
Solución en el dominio
Para validar el correcto funcionamiento del solver de Roe para MHD implementado,
corrimos el caso del tubo de choque MHD de (Brio y Wu, 1988). Este tubo de choque se
basa en el test de Sod, con la diferencia de que se agrega una discontinuidad en una de las
componentes transversales del campo magnético B. Las condiciones iniciales para este caso
son:
[ρ, u, v, w, By , Bz , p]l = (1, 0, 0, 0, 1, 0, 1)
1 1
[ρ, u, v, w, By , Bz , p]r = 8, 0, 0, 0, −1, 0, 10
Gallice, 1997). La solución consiste de tres ondas que viajan hacia la derecha y dos que
expansión.
puesta anteriormente mencionada. Ésta es una estructura de onda propia del sistema MHD;
aparece en esta condición porque el sistema MHD es no convexo. Cuando By se hace nu-
lo la onda magnetosónica lenta cambia de carácter, como lo explican (Brio y Wu, 1988;
Goedbloed y Poedts, 2004). Para una discusión más profunda de las distintas estructuras
de ondas posibles en MHD puede consultarse la tesis de (Yalim, 2008) y sobre las posibles
Como primera medida, vericamos que el esquema resuelva correctamente las ondas
que ocurren dentro del dominio. Corrimos un primer caso empleando 400 celdas de longitud
L=2 y número CF L = 0,5. En la Figura 6.33 se muestran los resultados obtenidos para
distintos intervalos de tiempo. Los resultados se comparan con la solución para t = 80 con
Se observa que el código captura las cinco ondas del sistema, por lo que se compor-
presibilidad articial.
discontinuidades. Esto hace que no aparezcan resultados espurios para la velocidad w y para
el campo magnético Bz (que deben mantenerse iguales a cero durante toda la simulación).
Esto no ocurre con solvers de Riemann aproximados que usan promedios aritméticos.
Condiciones de contorno
Una vez que vericamos que el algoritmo resuelve correctamente la solución dentro del
dominio, analizamos lo que ocurre cuando las ondas interactúan con distintos tipos de con-
torno. Para ello comparamos las soluciones obtenidas con el esquema de BCs basado en
características con las obtenidas con condiciones de extrapolación de orden cero. Nueva-
mente, comparamos los resultados para un contorno de pared cerrada y un contorno abierto
0.8 0.8
0.6 0.6
ρ
p
0.4 0.4
0.2 0.2
0 0
0 200 400 600 800 1000 0 200 400 600 800
x x
Velocidad adimensionalizada Campo Magnetico en direccion Y
1 1.5
0.5 0.5
u
By
0 −0.5
−1
−0.5 −1.5
0 200 400 600 800 0 200 400 600 800 1000
x x
1 51 101 151 201 sol C−G
manera que para el modelo gasdinámico (mediante la Ec. (6.27)), con las del esquema
Ecs. (5.72) y (5.73). Para las demás ondas implementamos condiciones de no reexión, de
La Figura 6.34 muestra los perles de densidad obtenidos para distintos instantes:
que viaja hacia la derecha llega a la pared, y es reejado de manera correcta por ambos
esquemas. Se observa que para t = 0,2 existe una onda de expansión que viaja desde la
Sin embargo, en los instantes de tiempo siguientes aparecen diferencias entre los resul-
tados asociados a cada esquema. Encontramos que las mencionadas discrepancias ocurren
168 Capítulo 6. Resultados
0.9
0.8
0.7
0.6
ρ
0.5
0.4
0.3
Ex t=0.2
0.2 Ex t=0.4
Ex t=0.6
Char t=0.2
Char t=0.4
0.1 Char t=0.6
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x
cuando ondas que transportan perturbaciones en la componente transversal del campo mag-
nético By interactúan con las paredes. Esto puede observarse claramente en la Figura 6.35,
velocidad v.
Creemos que las diferencias mencionadas pueden explicarse en virtud de una restricción
primitivas como:
∂v + u ∂v + Bx ∂By = 0
∂t ∂x ρ ∂x
(6.31)
∂By + B ∂u − B ∂v + u ∂By = 0
∂t y ∂x x ∂x ∂x
Para imponer una pared sólida es necesario que um+1 = −um , pero todas las demás
variables se extrapolan, o sea vm+1 = vm , Bym+1 = Bym . Esto implica que todas las
6.3. Casos Magnetohidrodinámicos 169
Ex t=0.2
0.8
Ex t=0.4
Ex t=0.6
Char t=0.2
Char t=0.4
0.6
Char t=0.6
0.4
0.2
By
−0.2
−0.4
−0.6
−0.8
−1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
X
respecto a x en la interfaz entre la última celda del dominio (celda m) y la primera del
∂u um
≈2 (6.32)
∂x m ∆x
∂By 2um
≈ −By (6.33)
∂t ∆x
Esto demuestra que la extrapolación impone condiciones sobre By , que en general son
diferentes por la evolución temporal obtenida en función de los operadores L, dada por la Ec.
(5.67e). Creemos que la restricción dada por esta última relación es más consistente porque
ρet del sistema. Por esta razón, una diferencia en By produce una diferencia signicativa en
la densidad y la presión. La Figura 6.36 muestra los distintos patrones de ondas obtenidos
0 0
1 1
t
t
−0.5 −0.5
0.5 0.5
−1 −1
0 0
0 0.5 1 0 0.5 1
X X
Extremo abierto
En este caso consideramos que el plasma sale a una región de estancamiento, donde la
presión p y la componente normal del campo manético Bx son constantes en todas partes.
Este caso presentó varias dicultades desde el punto de vista numérico, tanto para el
tación del esquema extrapolado debe hacerse una distinción a priori: si el ujo a la salida es
súper magnetosónico rápido, o sea u > cf , entonces todas las variables conservadas deben
ser extrapoladas. Pero si el ujo es sub magnetosónico rápido debe imponerse al menos una
de otras variables del sistema llevó a que el ujo se vuelva súper magnetosónico rápido
en las celdas del contorno. Una vez que el ujo se volvió súper magnetosónico rápido la
presión comenzó a incrementarse en el tiempo, hasta que el ujo se volvió una vez más
∂p
submagnetosónico rápido. Cuanbdo esto ocurre, el esquema
∂t =0 vuelve a imponer una
condición sobre la conservación del valor de la presión, pero conservando el nuevo valor.
dinámico, a pesar de que en este caso se trata de un choque magnetosónico lento. Este
sea la más apropiada. Sin embargo, no es sencillo determinar cuál es la solución física-
queda claro sobre qué variable debe imponerse qué condicion, debido a la complejidad del
problema.
Para salvar esta dicultad empleamos también el esquema que impone presión en el
campo lejano, dado por la Ec. (5.43). Sin embargo, este esquema debió ser nuevamente
modicado para representar los mismos patrones de ondas antes de que el ujo se torne
Para ello fue necesario establecer un valor para σ = 25000. En las Figuras 6.37 y 6.38
expansión magnetosónica rápida que viaja hacia la derecha casi no se percibe por su pequeña
intensidad.
Por otro lado, comparamos el equema basado en características para p∞ con el esquema
de BCs extrapoladas. En las Figuras 6.39 y 6.40 comparamos los contornos de densidad y
velocidad obtenidos.
pleando este esquema también obtuvimos ujo magnetosónico rápido una vez que el choque
esquema extrapolado impone de facto el valor inicial de la presión en el contorno, una vez
Nuevamente, los patrones de ondas obtenidos con ambos esquemas de BCs dieren,
172 Capítulo 6. Resultados
0.9 0.9
0.35 0.35
0.8 0.8
0.3 0.3
0.7 0.7
0.25 0.25
0.6 0.6
t
t
0.2 0.2
0.5 0.5
0.15 0.15
0.4 0.4
0.1 0.1
0.3 0.3
0.05 0.05
0.2 0.2
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1
X X
Figura 6.37: Comparación de contornos de densidad del test de Brio-Wu con extremo
abierto para esquemas de BCs basados en características
0.4 0.4
0.3 0.3
0.2 0.2
0 0
0.15 0.15
−0.2
−0.2
0.1 0.1
−0.4
0.05 0.05 −0.4
0 −0.6 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
Figura 6.38: Comparación de contornos de velocidad del test de Brio-Wu con extremo
abierto para esquemas de BCs basados en características
especialmente en la región afectada por la reexión del choque magnetosónico lento. Sin
embargo, esta imposición abrupta de la presión por parte del esquema extrapolado hace
que el algoritmo interprete que existe una discontinuidad en el contorno, lo cual creemos no
soluciones más suaves y estables en el tiempo, con cambios que ocurren más gradualmente.
Podemos observar estos fenómenos en la Figura 6.41. Además, el esquema basado en carac-
terísticas muestra más reexiones de ondas desde el contorno derecho, sobre todo después
6.3. Casos Magnetohidrodinámicos 173
t
0.6 0.5 0.5
0.4 0.4 0.4
0.4
0.3 0.3
0.2 0.2
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.5 1
X X
0.6
0.4
0.4 0
0
0.2 0.2
−0.5
0 0
0 0.2 0.4 0.6 0.8 1 0 0.5 1
X X
t
t
0.6
0.4 0.4 0.4
0.4
0.2 0.2
0.2 0.2
0 0
0 0.5 1 0 0.5 1
X X
la información que se posee sobre la dinámica de los arcos coronales se obtiene mediante
telescopios que miden en el espectro EUV SXR. Ellos permiten determinar indirectamente
existen ondas que se propagan dentro del arco, conocidas como abrillantamientos, que se
Alfvén. Estas perturbaciones se propagan desde la base hacia arriba, con períodos de entre
5 y 20 minutos, con intensidades que van decayendo a medida que la onda se acerca al ápice.
Las ondas que descienden desde el ápice hacia abajo son de menor intensidad, debido a los
han observados eventos en espectro EUV donde los abrillantamientos se propagan desde
arriba hacia las bases, para luego rebotar hacia arriba con una intensidad aún mayor. Los
podría esperarse que para cambios de densidad del orden de un 40 % a través de una onda,
sieron que los abrillantamientos pueden explicarse como ondas de choque que se propagan
asumiendo que las pérdidas de energía por radiación eran compensadas por el transporte
de calor por conducción y el efecto del calentamiento, tanto local como globalmente. Sin
embargo, no pudo lograrse que el modelo reprodujera el rebote de las ondas desde la base,
ciones iniciales más realistas basadas en soluciones hidrostáticas con términos fuente en
génea.
bases o en el ápice.
tica, y se emplearon tanto el modelo gasdinámico como el modelo MHD. Para el modelo
transversales nulas. Dichas componentes del campo magnético son nulas ya que los arcos
coronales son estructuras fuertemente dominadas por el campo magnético, donde las líneas
Ec. (3.155). En este caso la onda magnetosónica lenta se transforma en una onda sónica,
y la onda magnetosónica rápida en una onda de Alfvén trivialmente nula. Por esta razón,
a que se calcula con la velocidad característica mayor |u| + |cf |, que en este caso está
asociada a una onda trivial. En las Figuras 6.42 y 6.43 analizamos un tubo de choque con
ρ = 1, u = 0 ∀x
p = 0,1 0 ≤ x ≤ 0,4 ∧ 0,6 ≤ x ≤ 1 (6.36)
p = 1 0,4 ≤ x ≤ 0,6
Los valores adimensionalizados que empleamos para el campo magnético longitudinal fueron
de Bx = 0, 1, 2. La solución consiste en dos ondas de choque que viajan hacia los extremos
del tubo, con dos ondas de expansión que viajan hacia el centro a medida que las ondas de
3 0.6
0.4
2.5
0.2
2
1.5
u
ρ
−0.2
1
−0.4
0.5
−0.6
0 −0.8
0 2 4 6 8 10 0 2 4 6 8 10
Euler t=0,26 MHD B=0 t=0,26 MHD B=1 t=0,26 MHD B=2 t=0,26
Figura 6.42: Perles de densidad y velocidad normal para el tubo de choque (6.36)
con distintos valores de Bx
Dichas guras ponen en evidencia lo explicado anteriormente, esto justica que muchos
autores, como (Müller et al., 2003) empleen modelos gasdinámicos para llevar a cambo
una dimensión ya es necesario incluir los efectos de los campos magnéticos, para contemplar
0.55 4
0.5
3.5
0.45
3
0.4
2.5
0.35
ρet
0.3
p
0.25
1.5
0.2
1
0.15
0.5
0.1
0.05 0
0 1 2 3 4 5 6 7 8 9 10
0 2 4 6 8 10
Euler t=0,26 MHD B=0 t=0,26 MHD B=1 t=0,26 MHD B=2 t=0,26
Figura 6.43: Perles de presión y energía para el tubo de choque (6.36) con distintos
valores de Bx
un factor arbitrario Kcal . El factor Kcal toma valores de 2, 5, 10, 15 y 20. En todos los casos
empleamos condiciones de contorno basadas en características de entrada con velocidad nula
de pasos de tiempo varió dependiendo del caso analizado, siendo el valor máximo para
Kcal = 20 de 300.000 pasos, para lograr una simulacion de unos 3500 segundos.
desde las bases hacia el ápice. Dicho incremento de velocidad produce ondas que a su
simétricas. En la Figura 6.45 los perles de velocidad para distintos instantes de tiempo en
de velocidad.
5
x 10 9
18 10
16
14
12
10
log (ρ)
T[K]
E =1
6 H0
EH0=2
EH0=5
4
E =1 F no lim
H0 C
EH0=2 FC no lim
2
EH0=5 FC no lim
8
0 10
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x[cm] x 10
9 x[cm] 9
x 10
Se observa que para estos casos el modelo de conducción térmica (limitado o no limitado)
no tiene mucha inuencia en la solución de equilibrio estacionario. Creemos que esto se debe
cantidad de calor aportada por EH no llega a ser suciente para generar temperaturas
mucho más elevadas, por lo que el valor de Qcond en la mayor parte del dominio es mucho
de recorrer varias veces el arco toman una forma similar a una sinusoide. Creemos que la
difusión juega un papel fundamental para hacer que las ondas mencionadas posean esa forma
más suave. Su amplitud se incrementa durante los primeros 200 segundos hasta alcanzar un
valor máximo, para luego amortiguarse. Resulta llamativo además que la velocidad crece
sólo en el sentido desde la base hacia el ápice, existiendo velocidades mucho menores en la
dirección descendente para Kcal = 2 y siendo las velocidades ascendentes mucho mayores
que las descendentes para Kcal = 5. Creemos que esto puede explicarse debido al efecto de
más eciente por conducción térmica que por efectos convectivos. Esta idea es apoyada por
la Figura 6.45, donde se observa que las velocidades descendentes para el ujo de calor
6.4. Arcos magnéticos con perturbaciones 179
no limitado son mayores que aquellas para el modelo con ujo limitado (t = 600 s). Este
resultado además estaría de acuerdo con lo que discuten (Nakariakov y Verwichte, 2005)
EH= 2 EH =5
4000 8000
t=60s t=60s
t=200s t=200s
t=400s t=400s
3000 t=800s 6000 t=600s
t=60s FC no lim t=60s FC no lim
t=200s FC no lim t=200s FC no lim
t=400s FC no lim t=400s FC no lim
2000 4000
t=800s FC no lim t=600s FC no lim
1000 2000
u [m/s]
u [m/s]
0 0
−1000 −2000
−2000 −4000
−3000 −6000
−4000 −8000
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
Si bien el proceso de evolución del sistema es complejo, pueden verse algunas caracte-
punto representativo. Para ello gracamos u en x = L/2, donde se produjeron los máximos
de las velocidades, y la presión, densidad y temperatura en el centro del arco, esto es x = L.
Dichos resultados se muestran en la Figura 6.46. En esta gura también se evidencian los
cadas cuanto menor sea Kcal ) y su eventual amortiguamiento. Suponemos que el principal
efecto causante del amortiguamiento es la fuga de calor por las bases, que disipa la energía
de las ondas mecánicas que se propagan hacia arriba cuando éstas retornan hacia las bases.
ondulatorio observado, y los períodos de las ondas están dentro del margen posible regis-
de conducción térmica elegido, quizás por las mismas razones por las que que la solución
180 Capítulo 6. Resultados
cociente entre L/2 y el tiempo en el que ocurre el valor máximo de velocidad para x = L/2
(al que asociamos con el paso de la onda). Al período τ lo estimamos determinando el
tiempo en el que la velocidad en x = L/2 hace un ciclo completo, es decir, aumenta hasta
llegar a un máximo para luego disminuir hasta alcanzar un mínimo, y volver a cero.
L/2
Vonda ≈ (6.37)
tumax
5
x 10
0.09 10
0.085
8
0.08
0.075 6
p[dyn/cm2]
0.07
4
u[cm/s]
0.065
2
0.06
0.055 0
0.05
−2
0.045
0.04 −4
0 500 1000 1500 2000 2500 3000 3500 0 500 1000 1500 2000 2500 3000 3500
t[seg] 6 t[seg]
x 10
8 x 10
1.8 1.8
1.7
1.75
1.6
1.7 1.5
dens[part/cm 3]
1.4
T[K]
1.65
1.3
1.6 1.2
1.1
1.55
1
1.5 0.9
0 500 1000 1500 2000 2500 3000 3500 0 500 1000 1500 2000 2500 3000 3500
t[seg] t[seg]
EH0=1 EH0=2 EH0=5 EH0=1 FC no lim EH0=2 FC no lim EH0=5 FC no lim
Sin embargo, notamos que los cambios en la densidad a través de la onda para Kcal = 5
son del orden del 5 % y en la temperatura no se ven cambios apreciables, luego no creemos
que dichas ondas puedan asociarse a abrillantamientos. Otra particularidad que notamos
la presión y la densidad. Esta última presenta un incremento constante con muy baja
a las pérdidas por radiación. Se observa que las pérdidas radiativas no son prácticamente
afectadas por el incremento de temperatura del sistema. Esto se debe a que la función
Λ(T ) (Ec. (3.2)) es constante hasta T ≈ 2 · 106 K, y a partir de ese valor decrece con la
temperatura. Además, como la densidad no varía en gran medida, concluimos que cualquier
exceso de calor que posea el sistema debe equilibrarse mediante conducción térmica o a
−4 −4
x 10 x 10
2.5 1
EH0=1
2 0.9 EH0=2
E =5
H0
0.8
1.5 E =1 F no lim
H0 C
EH0=2 FC no lim
0.7
1 EH0=5 FC no lim
QCond [erg/cm /s]
0.6
3
0.5
0.5
0
0.4
−0.5
0.3
−1
0.2
−1.5 0.1
−2 0
0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8
x[cm] 9
x 10 x[cm] 9
x 10
Kcal mayor a 5
Cuando el parámetro Kcal toma valores mayores a 5 las oscilaciones en el arco presen-
tan amplitudes mucho mayores, y las velocidades en la dirección descendente toman valores
más cercanos a las velocidades en la dirección ascendente. Se observa que los efectos convec-
similar a una solución hidrostática. Creemos que el sistema no está en un estado estaciona-
rio porque los grácos de convergencia siguen mostrando cambios cada ciertos intervalos de
tiempo. Sin embargo, es difícil aseverar si las pequeñas velocidades de las últimas oscilacio-
nes son amortiguadas por la viscosidad numérica del esquema de volúmenes nitos o por
182 Capítulo 6. Resultados
algún proceso físico. En la Figura 6.48 se muestran los perles de temperatura y densidad
6
x 10 9
3 10
2.5
log(ρ)
T [K]
1.5
1 EH0=10
E =15
H0
EH0=20
EH0=10 FC no lim
0.5
E =15 F no lim
H0 C
E =20 F no lim
H0 C
8
0 10
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x[cm] 9 x [cm] 9
x 10 x 10
peratura en el centro del arco. Esto se explica debido a que para el modelo de ujo de
calor limitado la cantidad de calor transportada hacia las bases es menor. Luego, el sistema
posee más energía en la zona central del arco, que no puede ser evacuada hacia las bases.
en las ondas observadas. Los períodos son más cortos y las velocidades mayores debido a
la mayor energía involucrada. Una diferencia notable es que, para valores elevados de Kcal
a partir de la segunda oscilación comienza a haber un desfasaje entre las ondas asociadas
períodos más cortos pero amplitudes menores. Sin embargo, resulta llamativo que la distri-
la cantidad de calor depositada. A pesar de ello, a través de las ondas ocurren variaciones
6.4. Arcos magnéticos con perturbaciones 183
6
x 10
0.15 2
1
p[dyn/cm2]
u[cm/s]
0.1
0
−1
0.05
−2
0 500 1000 1500 0 6
500 1000 1500
8
x 10 t[seg] x 10 t[seg]
2.1
2 2.5
dens[part/cm 3]
1.9
2
T[K]
1.8
1.7
1.5
1.6
1.5 1
0 500 1000 1500 0 500 1000 1500
t[seg] t[seg]
EH0=10 EH0=15 EH0=20 EH0=10 FC no lim EH0=15 FC no lim EH0=20 FC no lim
En estos casos observamos que las velocidades máximas se incrementan cada vez menos
si Kcal toma valores mayores a 10, obteniéndose períodos cada vez más cercanos a un
hipotético valor asintótico. Esto permite conjeturar que el sistema llega a un estado de
Tabla 6.2: Períodos y velocidades de ondas estimadas para distintos valores de Kcal
con sH = ∞
modelo de conducción de calor no tiene mayor inuencia en las velocidades del sistema.
Creemos que esto ocurre porque para estos casos la energía depositada es tan elevada que
De la misma manera que para los valores de Kcal menores a 5, las pérdidas por radiación
184 Capítulo 6. Resultados
Lrad no se vieron inuenciadas por la cantidad de calor depositada. Las curvas asociadas
al ujo de calor por conducción fueron análogas a las del caso anterior, pero con un valor
Por otro lado, comparamos los resultados obtenidos con los que obtuvimos con el esque-
ma de condiciones extrapoladas análogas a las empleadas para el tubo de choque con pared
sólida (Ec. (6.27)), para evaluar el funcionamiento del esquema característico desarrollado
para esta tesis. Los resultados fueron llamativamente diferentes, tanto en la conguración
de estado estacionario del sistema, como en las ondas que se formaron durante el transi-
torio. En la Figura 6.50 se muestran los perles de temperatura y densidad para el estado
estacionario nal.
6
x 10 11
2.5 10
10
10
1.5
logρ
T[K]
1
9
10
0.5
8
0 10
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x[cm] x[cm] 9
9
x 10 x 10
Nuevamente, para lograr que el sistema converja a una solución fue necesario establecer
interfaz de las celdas del dominio con el contorno (mediante la Ec. (5.2)). Esto, como se
velocidad nula en las caras, produce una discontinuidad de velocidades en los contornos.
tratara de una salida supersónica. Este resultado no puede ser físicamente correcto, ya
que en todos los casos el estado del sistema en la cromósfera implica condiciones de ujo
subsónico. Además, los valores de densidad y presión en el interior del dominio son menores
a los obtenidos por la solución paramétrica, contrario a lo que se esperaría debido a que se
le está introduciendo energía al sistema y a que la masa debe conservarse. En la Figura 6.51
u(x = L/2).
8 6
x 10 x 10
2 2
1.9
1.5
1.8
1
1.7
0.5
dens[part/cm 3]
1.6
u[cm/s]
1.5
−0.5
1.4
−1
1.3
−1.5
1.2
1.1 −2
0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200
t[seg] t[seg]
EH0=2 ex EH0=5 ex EH0=10 ex EH0=2 ch EH0=5 ch EH0=10 ch
presión sin mostrar el comportamiento cuasi periódico obtenido con el modelo basado en
dirección opuesta a la del otro modelo de BCs debido a que la discontinuidad de velocidades
de velocidades para distintos instantes de tiempo. Para valores reducidos del parámetro de
E =2 4 EH =5
H
x 10
4000 1
3000
2000 0.5
1000
u [m/s]
u [m/s]
0 0
−1000
−2000 −0.5
−3000
−4000 −1
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
EH= 10 4 EH =15
4
x 10 x 10
2 3
1.5
2
1
1
0.5
u [m/s]
u [m/s]
0 0
−0.5
−1
−1
−2
−1.5
−2 −3
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
Figura 6.52: Perles de velocidad para BCs extrapoladas (línea contínua) y caracte-
rísticas (puntos) con Kcal = 2, 5, 10, 15 y sH = ∞
que la descendente. Para valores mayores de Kcal los perles de velocidad resultaron más
simétricos. Mientras mayor sea dicho parámetro, más sinusoidal es la forma de la distribu-
ción. Incluso, con este modelo de condiciones se obtuvieron velocidades mayores, a pesar de
existir gradientes de presión menores. Creemos que esto puede explicarse porque el modelo
la densidad en las celdas del contorno y al poseer ésta una velocidad igual y opuesta a
las celdas del dominio, existirá un ujo neto de cantidad de movimiento de los contornos
hacia el interior del dominio. Este inconveniente no existe con el modelo de BCs basadas
en características, que justamente satisface las leyes de conservación. Por estas razones,
La siguiente hipótesis que probamos fue considerar un arco de las mismas dimensiones,
con idénticas condiciones de contorno, y con una función de calentamiento que deposita
6.4. Arcos magnéticos con perturbaciones 187
calentamiento dado por la Ec. (3.15), e igualando a la cantidad de calor depositada por una
paramétrica explicada en la Sección 6.2. Si se emplea dicha solución como condición inicial
con una temperatura máxima un 20 % mayor, pero con una distribución espacial idéntica.
casos extremos. Para este modelo utilizamos el modelo de ujo de calor limitado, ya que
produjo períodos más cortos en los casos anteriores (lo consideramos el caso más desfavo-
rable). Además, de acuerdo a la literatura sería el caso físicamente más realista (Petralia
et al., 2014).
En todos los casos el sistema llega a una condición de pseudo equlibrio, donde la con-
el sistema evoluciona por escalones, pasando por distintos estados de pseudo equilibrio y
le toma bastante más tiempo alcanzarlo. En la Figura 6.53 se muestran las distribuciones
tura casi isotérmicas en la parte central del arco, con fuertes gradientes cerca de las bases.
Como es de esperarse, las temperaturas máximas son menores que las obtenidas para el
caso uniforme a igualdad de energía depositada. Las curvas de densidades para el nal de
la simulación son similares para todos los casos, presentando algunas ondas que aún se
propagan.
de propagación de ondas mucho mayores y períodos aún más cortos. Aparecieron pertur-
baciones en las velocidades que parten desde las bases hacia el ápice del arco, pero para
este caso fueron de mayor amplitud y con forma de diente de sierra. A medida que el
parámetro Kcal fue aumentando, incrementó el valor de las velocidades locales hasta llegar
188 Capítulo 6. Resultados
6
x 10 9
2 10
1.8
1.6
1.4
1.2
logρ
T [K]
8
1 10
0.8
E =1
0.6 H0
EH0=2
0.4 EH0=5
EH0=10
0.2 EH0=15
E =20
H0
7
0 10
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
x[cm] 9 x[cm] 9
x 10 x 10
a un valor de saturación en la velocidad local del ujo para Kcal = 15. Sin embargo, la velo-
cidad de propagación de las ondas sigue aumentando ligeramente, aunque presumiblemente
tendiendo a un valor asintótico como en los casos anteriores. En los primeros instantes de
tiempo las ondas se propagaron con la misma velocidad, pero a medida que se produjeron
rebotes aparecieron defasajes en las velocidades de propagación. Sin embargo, los valores
máximos de las velocidades y presiones fueron inferiores a los del caso con calentamiento
peratura en x=L y la velocidad en x = L/2. Allí pueden observarse los desfasajes ante-
de una onda. Resulta llamativo que mientras mayor sea la cantidad de calor depositada,
más lenta se vuelve la propagación de las ondas luego de los rebotes. Se observa que en
estos casos existen cambios signicativos en la densidad (del orden de hasta un 100 %) y
en la temperatura a través de las ondas, por lo que es posible considerar estas ondas como
con la solución hidrostática a la que converge para Kcal = 1. Los períodos y velocidades
obtenidas están cerca del límite inferior de lo que normalmente se acepta en la literatura,
EH= 2 4 EH= 5
4
x 10 x 10
4
5
2
u [m/s]
u [m/s]
0 0
−2
−5
−4
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
5 EH=10 5 EH =15
x 10 x 10
1 1
0.5 0.5
u [m/s]
u [m/s]
0 0
−0.5 −0.5
−1 −1
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L t=60s t=200s t=400s t=600s x/L
Figura 6.54: Perles de velocidad para arcos coronales con Kcal = 2, 5, 10, 15 y sH =
0,1L
6
x 10
0.08 10
0.06 5
p[dyn/cm2]
u[cm/s]
0.04
0
0.02
−5
0
0 8 500 1000 1500 2000 0 6 500 1000 1500 2000
x 10 t[seg] x 10 t[seg]
2.5
1.5 2
dens[part/cm 3]
1.5
T[K]
1
1
0.5
0.5
0
0 500 1000 1500 2000 0 500 1000 1500 2000
t[seg] EH0=1 EH0=2 EH0=5 EH0=10 EH0=15 t[seg]
E =20
H0
Figura 6.55: Evolución temporal de p(x = L), ρ(x = L), T (x = L) y u(x = L/2)
para Kcal = 2, 5, 10, 15 y sH = 0,1L
190 Capítulo 6. Resultados
Tabla 6.3: Períodos y velocidades de ondas estimadas para distintos valores de Kcal
con sH = 0,1L
La siguiente hipoótesis que probamos fue asumir que el arco sufre una deposición brusca
de energía en la zona del ápice. Dicho proceso fue modelado como un incremento discon-
tínuo en la presión para la solución hidrostática de equilibrio, mediante la Ec. (6.39). Las
condiciones iniciales para la densidad y la velocidad son las mismas que para el problema
hidrostático
p = p
hidr , 0 ≤ x ≤ 0,4 ∧ 0,6 ≤ x ≤ 1
(6.39)
p = K p
p hidr , 0,4 ≤ x ≤ 0,6
(Fernandez et al., 2009). Como en los casos anteriores, modelamos ambos contornos como
entradas con velocidad (nula) impuesta y temperatura impuesta. Corrimos 6 · 105 pasos de
tiempo con CF L = 0,5 y F O = 0,2. También comparamos los dos modelos presentados
para el ujo de calor por conducción: el modelo clásico y el modelo limitado, obteniendo
Como en los casos anteriores, se obtuvieron soluciones simétricas para la presión, tem-
peratura y densidad para amobos modelos de conducción térmica. Sin embargo, puede verse
claramente que la estructura de las ondas obtenidas son muy diferentes en ambos casos.
Para el caso del ujo de calor no limitado, la conduccción térmica rápidamente redistribu-
presión impuesto por la gravedad), hasta llegar a los contornos, donde el sistema elimina
En la Figura 6.57 se muestran los perles de velocidad para los mismos instantes de
tiempo. Para el caso de ujo de calor no limitado, se propagan dos ondas de compresión
6.4. Arcos magnéticos con perturbaciones 191
Fc no limitado Fc limitado
0.025
0.02
0.02
0.018
0.016 0.015
p [Pa]
p [Pa]
0.014
0.012 0.01
0.01
0.005
0.008
0.006 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
t=5 s t=20 s t=40 s t=60 s t=80 s
0.02 0.02
0.018
0.016 0.015
0.014
p [Pa]
p [Pa]
0.012
0.01
0.01
0.008
0.005
0.006
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
t=100 s t=140 s t=180 s t=240 s
hacia los extremos, que no llegan a ser ondas de choque por los efectos de la difusión.
El paso de tiempo se encuentra en todos los casos más limitado por la difusión que por
debido al gradiente favorable de presión asociado a la gravedad. Por otro lado, para el
el instante inicial la conducción térmica combinada con los efectos convectivos genera una
compresión simétricas: dos que se propagan hacia los extremos, y dos que se propagan
hacia el centro. En la región central existen dos efectos que compiten: la difusión térmica
que disminuye la presión, y las dos ondas de compresión que viajan hacia adentro. Como
en los casos anteriores, la difusión térmica tuvo el rol más importante. Sin embargo, dicha
discontinuidad. Esto lleva a que la convección tenga un rol más importante, apareciendo
192 Capítulo 6. Resultados
FC no limitado 5 FC limitado
x 10
4 x 10
3 1.5
2 1
1 0.5
u[m/s]
u[m/s]
0 0
−1 −0.5
−2 −1
−3 −1.5
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
t=5 s t=20 s t=40 s t=60 s t=80 s
4 4
x 10 x 10
4 6
4
2
2
u[m/s]
u[m/s]
0 0
−2
−2
−4
−4 −6
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
t=100 s t=140 s t=180 s t=240 s
velocidades con valores del doble que para el caso del ujo no limitado.
Para descartar que los resultados asociados al ujo limitado sean espurios debido a un
problema del esquema numérico analizamos el efecto de emplear una discretización espacial
Figura 6.58.
Se observa que para ambas mallas los resultados son casi idénticos, sobre todo la velo-
cidad de propagación del frente de onda. Luego, concluimos que si existe un problema con
Como en los casos anteriores, eventualmente las ondas del sistema terminan siendo
amortiguadas y el sistema converge a una solución similar a una hidrostática. Con las tres
6.4. Arcos magnéticos con perturbaciones 193
5
x 10
1.5 0.025
1
0.02
0.5
0.015
u [m/s]
p [pa]
0
0.01
−0.5
0.005
−1
−1.5 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
x/L x/L
t=5s m=400 t=20s m=400 t=40s m=400 t=60s m=400 t=5s m=200 t=20s m=200 t=40s m=200 t=60s m=200
de (Fernandez et al., 2009), por lo que podemos asumir que será necesario analizar arcos
7.1. Conclusiones
En la presente tesis nos propusimos desarrollar un modelo numérico unidimensional
para la simulación numérica de la dinámica de los arcos de la corona solar, que a su vez
Para poder incluir dichos efectos y lograr que el algoritmo numérico fuera estable fue
de (Yee et al., 1985), junto con el ujo numérico de (Yee, 1989). Dichas herramientas
son bajas.
tas implementadas, para luego comparar los resultados con los del modelo MHD. Para el
modelo MHD fue necesario normalizar los vectores y valores propios del sistema, para evi-
tar indeterminaciones. Además, utilizamos un solver de Roe exacto para MHD que no se
había empleado anteriormente en el grupo de trabajo. Dicho solver agregó robustez y pre-
cisión al algoritmo. Por otro lado, todos los algoritmos fueron implementados en variables
adimensionales. Esto se hizo para evitar problemas numéricos por diferencias de orden de
magnitud entre las diferentes ecuaciones que componen el sistema de leyes de conservación,
195
196 Capítulo 7. Conclusiones y trabajos futuros
y para poder comparar el efecto de los términos fuente con el de los términos convectivos.
Probamos dichos esquemas con problemas muy exigentes desde el punto de vista nu-
mérico, como el arco de (Cargill y Priest, 1980) (Sección 6.1.2) o el modelo de advección
difusión compresible (Sección 6.1.3). Ambos casos poseen soluciones subsónicas con una
importancia relativa de los términos fuente muy severa, que exigen al integrador tempo-
ral. Sin embargo el algoritmo convergió a la solución estacionaria, incluso con valores de
CF L mayores a uno y en condiciones próximas a las que no existe una solución única
del problema, para el caso del arco de Cargill. Podemos asegurar entonces que el esquema
complejidad.
Otro desafío que encontramos en el desarrollo de la tesis fue obtener una forma consis-
cercanas a las bases de los arcos coronales existen fuertes gradientes de densidad y tem-
peratura. Además, el arco intercambia energía y masa con la cromósfera; existen ujos de
contorno resultó inadecuado para modelar estos fenómenos. Para ello utilizamos un esque-
posibilitar imponer diversas restricciones físicas a los contornos, tanto para entradas como
salidas de ujo. Evaluamos diferentes condiciones posibles para benchmarks conocidos o so-
luciones analíticas especialmente seleccionadas para validar los modelos físicos necesarios.
En todos los casos el esquema basado en características convergió con la cantidad de condi-
las curvas de convergencia para problemas estacionarios. Dicho esquema también probó ser
exitoso en casos donde la conducción de calor y otros términos fuente tienen valores com-
parables a los de los ujos convectivos, empleando condiciones sobre la onda característica
que transporta el ujo de entropía. Esto permitió mantener la consistencia con la física en
todos los casos. Además, dicho esquema logró acoplarse de manera relativamente sencilla
Este esquema además resultó ser más exible y físicamente consistente para problemas
choque a regiones de uido sin perturbar. Probamos diferentes condiciones posibles, ha-
ciendo ajuste no en algunas de ellas, lo que permitió adquirir experiencia en este tipo de
fenómenos transitorios.
7.2. Contribuciones originales 197
magnetohidrodinámico con todas sus variantes. En las pruebas que realizamos este esquema
demostró tener ventajas claras respecto de la consistencia con la física comparado con el
esquema extrapolado. Sin embargo, quedaron ciertos interrogantes sobre qué condiciones
se pueden imponer sobre cuáles de las ondas que componen el sistema. A pesar de ello,
estamos conados en que los fundamentos físicos de nuestra elección son sólidos, y como
del orden de los observados. Las soluciones obtenidas nos permitieron sacar conclusiones
importantes sobre la dinámica de la energía del sistema y el balance de los términos fuente.
al esquema tradicional extrapolado. Sin embargo, resta aún hacer un barrido de parámetros
del problema para vericar si alguno de los casos obtenidos coincide con observaciones de
A pesar de ello, el caso que modela una deposición de energía como un incremento brus-
dependientes del modelo de conducción de calor empleado. En ambos casos dicho fenómeno
tuvo gran inuencia en la dinámica del sistema, siendo el medio más efectivo de trans-
porte de energía. El modelo de ujo de calor limitado produjo resultados inesperados que
ría identicar casos donde los términos fuente y difusivos tengan menor importancia. Por
a la longitud del arco, luego para longitudes mayores éstos tendrían menos importancia.
de forma aproximada permitió estudiar las propiedades de convergencia de cada uno de los
esquemas descritos y qué esquemas fueron más apropiados para distintos tipos de problema.
Las conclusiones de los diferentes casos gasdinámicos que probamos (tanto estacionarios
como inestacionarios) fueron vertidas en el trabajo (Cimino et al., 2015b), aún en revisión.
Dichas herramientas constituyen una base para continuar con una línea de investigación en
MHD es una innovación en sí misma, puesto que no existen muchos trabajos en esta área, y
es un tema muy vigente en la comunidad cientíca. El desarrollo del esquema y las pruebas
modelo tiene muchas aplicaciones a problemas astrofísicos (como los arcos coronales, las
contorno.
Finalmente, los resultados obtenidos para los arcos coronales con perturbaciones cons-
tituyen una extensión y una mejora del modelo utilizado por el grupo de trabajo. A su
vez permitieron plantear como posible explicación de los fenómenos observados efectos de
acoplamiento entre transferencia de calor, ondas del tipo convectivo y deposiciones de ener-
gía. También resulta novedoso el efecto que tiene el modelo de conducción térmica en las
pormenorizado pasaremos en limpio las conclusiones, que serán vertidas en una publicación
que se encuentra en proceso de escritura, ya que los resutados que obtuvimos con valores
a cada uno de ellos, y poder encontrar correlaciones posibles entre las velocidades de las
ondas y los mismos. Además, resta determinar si alguno de los casos obtenidos se asemeja
A pesar del éxito obtenido, el modelo unidimensional tiene algunas limitaciones, so-
7.3. Trabajo futuro 199
bre todo cuando se trata del sistema MHD. En este caso sólo existen ondas transversales
asociadas a las componentes transversales del campo magnético, que no pueden estar pre-
sentes en este tipo de modelo para los arcos coronales. Además, el hecho de que exista
de propagación de ondas no permite obtener todas las ondas posibles del sistema (switch
reconexión. El paso lógico siguiente sería extender el esquema a dos dimensiones. Para ello
será necesario hacer alguna transformación sobre el solver de Roe empleado para extenderlo
al caso bidimensional.
sea capaz de detectar las direcciones de propagación de las ondas en el contorno, para
imponer restricciones en las direcciones de las mismas. Esta extensión permitiría a su vez
vez podrían resolver los interrogantes sobre qué condiciones pueden imponerse sobre cuáles
ondas características.
Por otro lado, como en el grupo de trabajo también se emplea el programa FLASH para
mallado adaptativo que posee FLASH, permitiendo correr casos más grandes y complejos,
Otra posible línea de desarrollo sería la extensión al modelo MHD del código paralelo
integrante del grupo de trabajo. Esto también posibilitaría analizar casos más grandes
Esas derivadas son las asociadas a los operadores Li cuando las ondas son salientes. Si la
onda i-ésima es entrante,el operador asociado será cero en caso de condición de no reexión,
o estará expresado en función de otros L(j) para otro tipo de condición. Luego, su derivada
será tambien 0 (si se trata de una condición de no reexión) o una combinación lineal de
Debe tenerse en cuenta que, al emplearse derivadas descentradas para evaluar los ope-
radores L(j) en el contorno, las expresiones de estas derivadas variarán según se trate del
descentrada a la derecha
∂A Aj+1 − Aj
=
∂x ∆x
201
202 Apéndice A. Derivadas de los operadores Li para el modelo de Euler
∂L(1) (γ − 1) γuj 2 uj
1 aj
= − − [pj+1 − pj − aj ρj (uj+1 − uj )] (A.1a)
∂U1j ∆x 2ρj 4aj ρj ρj
(γ − 1) γuj 2 (γ − 1) uj 2
(uj − aj ) aj
+ − (uj+1 − uj ) − (uj+1 + uj ) −
∆x 4aj 2 2
(1)
∂L 1 (γ − 1) γuj 1
= + [pj+1 − pj − aj ρj (uj+1 − uj )] (A.1b)
∂U2j ∆x 2aj ρj ρj
(uj − aj ) (γ − 1) γuj (uj+1 − uj )
+ + (γ − 1) uj + aj
∆x 2 aj
∂L(1)
(uj − aj ) (γ − 1) γ (uj+1 − uj )
= − −γ+1 (A.1c)
∂U3j ∆x 2aj
(γ − 1) γ [pj+1 − pj − aj ρj (uj+1 − uj )]
−
2ρj aj ∆x
∂L(2) (γ − 1) γ (ρj+1 − ρj ) uj 2 (γ − 1) uj 2
uj 2 ρj+1 2
= − − + aj + aj (A.3a)
∂U1j ∆x 2 ρj 2 ρj
uj pj+1 − pj − a2j (ρj+1 − ρj )
−
∆xρj
∂L(2) pj+1 − pj − a2j (ρj+1 − ρj )
= + (A.3b)
∂U2j ∆x ρj
uj (γ − 1) γ (ρj+1 − ρj ) uj
+ (γ − 1) uj
∆x ρj
∂L(2)
uj (γ − 1) γ (ρj+1 − ρj )
= − −γ+1 (A.3c)
∂U3j ∆x ρj
∂L(2) (γ − 1) uj+1 2
uj
= − a2j (A.4a)
∂U1j+1 ∆x 2
∂L(2) (γ − 1)
=− uj uj+1 (A.4b)
∂U2j+1 ∆x
∂L(2) (γ − 1) uj
= (A.4c)
∂U3j+1 ∆x
204 Apéndice A. Derivadas de los operadores Li para el modelo de Euler
∂L(3) 1 (γ − 1) γuj 2
aj uj
= − − (pj+1 − pj + aj ρj (uj+1 − uj )) (A.5a)
∂U1j ∆x 4aj ρj 2ρj ρj
(γ − 1) γuj 2 (γ − 1) uj 2
(uj + aj ) aj
+ (uj+1 − uj ) + (uj+1 + uj ) −
∆x 4aj 2 2
(3)
∂L 1 1 (γ − 1) γuj
= − [pj+1 − pj + aj ρj (uj+1 − uj )] (A.5b)
∂U2j ∆x ρj 2aj ρj
(uj + aj ) (γ − 1) γuj (uj+1 − uj )
+ − + (γ − 1) uj − aj
∆x 2aj
∂L(3) (γ − 1) γ
= [pj+1 − pj + aj ρj (uj+1 − uj )] (A.5c)
∂U3j 2ρj aj ∆x
(uj + aj ) (γ − 1) γ (uj+1 − uj )
+ −γ+1
∆x 2aj
descentrada a la izquierda
∂A Aj − Aj−1
=
∂x ∆x
ya que las velocidades de ondas salientes irán hacia la derecha.
A.2. Extremo Derecho 205
∂L(1) (γ − 1) γ uj 2
1 aj uj
= − + − (pj − pj−1 − aj ρj (uj − uj−1 )) (A.7a)
∂U1j ∆x 4ρj aj 2ρj ρj
(γ − 1) γ uj 2 (γ − 1) uj 2
(uj − aj ) aj
+ − (uj − uj−1 ) + (uj + uj−1 ) +
∆x 4aj 2 2
(1)
∂L 1 (γ − 1) γuj 1
= + (pj − pj−1 − aj ρj (uj − uj−1 )) (A.7b)
∂U2j ∆x 2aj ρj ρj
(uj − aj ) (γ − 1) γuj (uj − uj−1 )
+ − (γ − 1) uj − aj (A.7c)
∆x 2 aj
∂L(1)
(uj − aj ) (γ − 1) γ (uj − uj−1 )
= − +γ−1 (A.7d)
∂U3j ∆x 2aj
(γ − 1) γ (pj − pj−1 − aj ρj (uj − uj−1 ))
− (A.7e)
2 aj ∆xρj
!
∂L(2) 1 (γ − 1) γ (ρj − ρj−1 ) uj 2 (γ − 1) uj 2 a2j (ρj − ρj−1 ) 2
= uj − + + − aj
∂U1j ∆x 2 ρj 2 ρj
uj pj − pj−1 − a2j (ρj − ρj−1 )
− (A.9a)
∆x ρj
∂L(2) pj − pj−1 − a2j (ρj − ρj−1 )
=
∂U2j ∆x ρj
uj (γ − 1) γ (ρj − ρj−1 ) uj
+ − (γ − 1) uj (A.9b)
∆x ρj
∂L(2)
uj (γ − 1) γ (ρj − ρj−1 )
= − +γ−1 (A.9c)
∂U3j ∆x ρj
(γ−1) uj−1 2
∂L(2) uj aj 2 − 2
= (A.10a)
∂U1j−1 ∆x
∂L(2) (γ − 1) uj uj−1
= (A.10b)
∂U2j−1 ∆x
∂L(2) (1 − γ) uj
= (A.10c)
∂U3j−1 ∆x
A.2. Extremo Derecho 207
∂L(3) (γ − 1) γ uj 2
1 aj uj
= − − (pj − pj−1 + aj ρj (uj − uj−1 )) (A.11a)
∂U1j ∆x 4 aj ρj 2ρj ρj
(γ − 1) γ uj 2 (γ − 1) uj 2
(uj + aj ) aj
+ (uj − uj−1 ) − (uj + uj−1 ) +
∆x 4 aj 2 2
(3)
∂L 1 1 (γ − 1) γ uj
= − (pj − pj−1 + aj ρj (uj − uj−1 )) (A.11b)
∂U2j ∆x ρj 2 aj ρj
(uj + aj ) (γ − 1) γ uj (uj − uj−1 )
+ − − (γ − 1) uj + aj
∆x 2 aj
∂L(3) (γ − 1) γ (pj − pj−1 + aj ρj (uj − uj−1 ))
= (A.11c)
∂U3j 2ρj aj ∆x
(uj + aj ) (γ − 1) γ (uj − uj−1 )
+ +γ−1
∆x 2 aj
!
∂L(3) (uj + aj ) (γ − 1) u2j−1 aj ρj uj−1
= − + (A.12a)
∂U1j−1 ∆x 2 ρj−1
∂L(3)
(uj + aj ) aj ρj
= (γ − 1) uj−1 − (A.12b)
∂U2j−1 ∆x ρj−1
∂L(3) (1 − γ)
= (uj + aj ) (A.12c)
∂U3j−1 ∆x
208 Apéndice A. Derivadas de los operadores Li para el modelo de Euler
Solutions and Scaling Laws of Coronal Loops. The Astrophysical Journal Supplement
Series, 142:269283.
coronal loops ii. a trace case study. Astronomy and Astrophysics, 441:327 335.
Brio, M. y Wu, C. (1988). An upwind dierencing scheme for the equations of ideal mag-
Cargill, P. J. y Priest, E. (1980). Siphon ows in coronal loops. I - Adiabatic ow. Solar
Physics, 65:251269.
Cargo, P. y Gallice, G. (1997). Roe Matrices for Ideal MHD and Systematic Construction
209
210 Bibliografía
Cimino, A. M., Krause, G. J., Elaskar, S. A., y Costa, A. (2015a). Characteristic boundary
conditions for magnetohydrodynamics: the brio-wu shock tube. Computers & Fluids
(en revisión).
Cimino, A. M., Tamagno, J. P., Elaskar, S. A., y Costa, A. (2015b). Some applications of
boundary conditions based on characteristics for gas dynamic ows with source terms.
Coen, O. (2008). The Solar Corona Through Numerical Eyes. PhD thesis, University of
Michigan.
Classical and Saturated Mass Loss Rates . The Astrophysical Journal, 211:135146.
Dedner, A., Kröner, D., Sofronov, I., y Wesemberg, M. (2001). Transparent boundary
Dutt, P. (1988). Stable Boundary Conditions and Dierence Schemes For Navier-Stokes
Einfeldt, B., Muntz, C., Roe, P., y Sjögreen, B. (1991). On godunov type methods near
Engquist, B. y Madja, A. (1977). Absorbing boundary conditions for the numerical simu-
Fernandez, C. A., Costa, A., Elaskar, S., y Schulz, W. (2009). Numerical simulation of the
internal plasma dynamics of post are loops. Monthly Notices of the Royal Astrono-
mical Society.
Handy, B. N., Acton, L. W., y Kankelborg, C. C. (1999). The transition region and coronal
Harten, A. (1983). High resolution schemes for hyperbolic conservation laws. Journal of
Computational Physics, 49:357393.
Bibliografía 211
Harten, A. y Hyman, J. M. (1983). Self adjusting grid methods for onedimensional hy-
Hayashi, K. (2005). Magnetohydrodynamic simulations of the solar corona and solar wind
using a boundary treatment to limit solar wind mass ux. The Astrophysical Journal
Supplement Series, 161.
Inan, U. y Golkowski, M. (2011). Principles of Plasma Physics for Engineers and Scientists.
Cambridge.
Matatsuka, K. (2013). I do like CFD, Vol 1. Governing Equations and Exact Solutions.
CFDbooks.
Müller, D. A. N., Hansteen, V. H., y Peter, H. (2003). Dynamics of solar coronal loops.
i. condensation in cool loops and its eect on transition region lines. Astronomy and
Astrophysics, 411:605613.
Ogawara, Y., Takano, T., Kato, T., Kosugi, T., y Tsuneta, S. (1991). The solar a mission:
Petralia, A., Reale, F., 2, S. O., y Klimchuk, J. A. (2014). MHD modeling of coronal loops:
Peyrard, P. F. y Villedieu, P. (1999). A roe scheme for ideal mhd equations on 2d adaptively
Roe, P. L. (1981). Approximate riemann solvers, parameter vectors, and dierence schemes.
Rosner, R., Tucker, W., y Vaiana, G. (1978). Dynamics of the Quiescent Solar Corona. The
Astrophysical Journal, 220:643655.
Rudinger, G. (1955). On the Reection of Shock Waves from an Open End of a Duct.
Serio, S., Peres, G., y Vaiana, G. (1981). Closed coronal structures ii. generalized hydrostatic
Shapiro, A. (1953). The Dynamics and Thermodynamics of Compressible Fluid Flow. Num-
ber v. 1 in The Dynamics and Thermodynamics of Compressible Fluid Flow. Ronald
Press Co.
namics, based on solution of the well-posed riemann problem. ASP Conference Series,
385.
[Link] (2004). Modeling articial boundary conditions for compressible ow. Ann.
Rev. Fluid. Mech.
Thompson, K. (1987). Time dependent boundary conditions for hyperbolic systems. Journal
of Computational Physics, 68:124.
Thompson, K. (1990). Time dependent boundary conditions for hyperbolic systems, ii.
Toro, E. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer
Verlag.
Yee, H., Warming, R., y Harten, A. (1985). Implicit total variation diminishing (tvd)