Physics">
Nothing Special   »   [go: up one dir, main page]

Deducción de Ecuaciones de Transporte Noeclásico en Un Tokamak Axisimétrico

Descargar como pdf o txt
Descargar como pdf o txt
Está en la página 1de 67

INFORME DEL SERVICIO SOCIAL

Deduccción de las Ecuaciones de Transporte


Neoclásicas en un Tokamak Axisimétrico

Cecilio Jorge Vega Vázquez


Universidad Autónoma Metropolitana
Unidad Iztapalapa
División de Ciencias Básicas e Ingenierı́a

Realizado en el Instituto Nacional de Investigaciones Nucleares


Departamento de Fı́sica
Gerencia de Ciencias Básicas

Asesor: Dr. César R. Gutiérrez Tapia


Proyecto: FI-001

Diciembre, 2014.
Contenido

0.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1 Teorı́a general de transporte 4


1.1 Descripción cinética del plasma . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Ecuaciones cinéticas del plasma . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3 Ecuación de continuidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4 Relajamiento del plasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.5 Descripción hidrodinámica de un plasma . . . . . . . . . . . . . . . . . . . . 20
1.6 Ecuaciones de la Magneto-hidrodinámica . . . . . . . . . . . . . . . . . . . . 22

2 Aproximación de fluido 28
2.1 Teorı́a general de momentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Transporte neoclásico y el código ASTRA 44


3.1 Sistema de coordenadas, promedios de superficie y flujos . . . . . . . . . . . 45
3.2 Ecuaciones de transporte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3 Fuentes y Sumideros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4 Condiciones Iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.5 Condiciones de frontera para la densidad y la temperatura . . . . . . . . . . 52
3.6 Condiciones de frontera para el flujo poloidal . . . . . . . . . . . . . . . . . . 53
3.6.1 Corriente preescrita del plasma . . . . . . . . . . . . . . . . . . . . . . 53
3.6.2 Búcle de voltaje preescrito . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.6.3 Ecuación del circuito externo . . . . . . . . . . . . . . . . . . . . . . . . 54
3.7 Cierre de las ecuaciones de equilibrio y transporte . . . . . . . . . . . . . . . 55
3.8 Parámetros del Tokamak TCABR . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.8.1 Análisis del transporte del Tokamak TCABR con el código ASTRA . . 56

1
Resumen

Se describe la deducción de las ecuaciones de transporte neoclásicas para un Tokamak


axisimétrico en el marco de las actividades del Servicio Social realizado en el Departa-
mento de Fı́sica del Instituto Nacional de Investigaciones Nucleares.
0.1 Introducción

El panorama energético actual, impone una gama de retos para las crisis que derivan
de fenómenos como el incremento desmedido en la población mundial, y el desarrollo
tecnológico. Crisis que surgen del agotamiento de recursos, la destrucción del medio
ambiente, la escasez alimentaria, la demanda masiva de energı́a, etc. La energı́a es un
factor determinante en los procesos de la vida moderna, el ritmo actual de consumo en
las actividades humanas requiere grandes cantidades de energı́a, y las fuentes actuales
se han visto, en algún punto, superadas por la demanda. El agotamiento de las fuentes
no renovables (combustibles fósiles) además del fuerte impacto ambiental que generan,
aunado a las limitantes de las energı́as renovables en la disponibilidad, densidad, y al-
macenamiento, nos obligan a buscar alternativas para la producción de potencia. Algu-
nas de las energı́as renovables requieren grandes extensiones de territorio, aspecto que
en un mundo con tendencia a la sobrepoblación es difı́cil de proveer. La densidad en-
ergética y dispónibilidad de las fuentes renovables limitan su empelo y fiabilidad. Estas
restricciones han conferido en la energı́a nuclear de fusión, la opción más viable para la
producción de energı́a eléctrica continua para los próximos años, debido a la densidad
de energı́a, la reducción en las emisiones y la disponibilidad del combustible. Recorde-
mos que actualmente acuerdos internacionales limitan la cantidad de emisiones que los
paı́ses suscriptores pueden liberar en la atmosfera. La producción de energı́a nuclear a
través del proceso de fisión se ha encontrado con un gran número de detractores, debido
a los problemas que implican los desechos generados durante el proceso. El proceso
de fusión nuclear a pesar de ser estudiado desde los años 50’s, ha tenido poco eco en
paı́ses en vı́as de desarrollo, siendo el más grande reto de la comunidad cientı́fica para
el presente siglo.
En 1979, Estados Unidos, la URSS, la comunidad Europa y Japón se juntaron bajo
el auspicio de la Agencia Internacional de Energı́a Atómica, y con la recomendación del
cientı́fico soviético Evgeniy Velikhov se trabajó en una serie de talleres sobre los reactores
basados en dispositivos tipo Tokamak (INTOR), los talleres se extendieron la década sigu-
iente para constatar la viabilidad del Tokamak. Actualmente, los paı́ses más desarrolla-
dos (U.S.A, UE, Rusia, China, Corea del Sur, India, Japón) se han unido en el Proyecto

1
ITER cuyo objetivo es desarrollar el primer reactor de fusión comercial en un periodo no
mayor a 30 años. El proyecto que actualmente se construye en Francia, ya es consider-
ado el más caro en la historia de la comunidad cientı́fica a la par de la estación espacial
internacional. De lograrse el objetivo, los beneficios que aportará a las naciones desar-
rolladoras tendrán un impacto sin precedentes, en la economı́a y polı́tica mundial, pues
se dispondrá de una fuente de energı́a limpia, abundante, y relativamente barata. Sin
embargo el gran salto tecnológico, cientı́fico y económico que aportará el proyecto, puede
aumentar la dependencia hacia los paı́ses desarrollados en una escala importante. Por
lo que es inminente comenzar a desarrollar todo tipo de recursos en este campo. En los
últimos años la factibilidad cientı́fica del confinamiento magnético ha sido demostrada, y
se espera que la del confinamiento inercial sea establecida en un corto periodo de tiempo.
Los Tokamaks son un dispositivo de confinamiento cerrado en el cual un fuerte campo
toroidal Bφ es producido en el plasma por un conjunto de bobinas que encierran al
toro mediante el ensamble de varios segmentos de estas bobinas, formando un plasma
toroidal. El campo es producido por la corriente que fluye a través de las bobinas, y en
la región del plasma puede describirse como:
Bφ0
Bφ (r, θ) =  
r
1+ R0 cos θ

La teorı́a neoclásica está considerada como uno de los pilares en la descripción de


plasmas contenidos en dispositivos toroidales y es un punto de partida básico para la
comprensión del proceso de transporte de partı́culas y energı́a en plasmas contenidos
dentro de campos magnéticos en equilibrio, donde el plasma se encuentra en un estado
no turbulento. Este punto de partida permite la comprensión de diferentes procesos que
se presentan en sistemas toroidales como el Tokamak. Los fenómenos de transporte en
un fluido como el plasma es un caso muy particular pues al componerse de partı́culas
con carga implica la consideración de fuerzas eléctricas y magnéticas generadas en el
plasma o por una fuente externa. Estas propiedades hacen que el estudio difiera al de
un gas ordinario, por lo que el análisis requiere de herramientas de análisis adicionales
para este tipo de fluido. La geometrı́a del sistema, además de inducir procesos en el
plasma, requiere el empleo de coordenadas de flujo especı́ficas raramente empleadas en
el análisis de lı́quidos o gases. La literatura disponible por lo general trata cada uno de es-
2
tos puntos por separado. El objetivo de este trabajo de servicio es tratar de conjuntar las
herramientas anaı́ticas y numéricas necesarias para la deducción del modelo de trans-
porte neoclásico en un solo documento. Se muesdtra en un ejemplo, un modelo para
simular el transporte en el Tokamak TCBAR mediante el código ASTRA. Este Tokamak
donado por suiza al gobierno brasileño es el experimento más desarrollado que actual-
mente se mantiene en operación en América latina. La importancia de este Tokamak
radica en el sentido que ningún paı́s de Latinoamérica forma parte del proyecto ITER, y
en este sentido, el TCBAR es el experimento más importante para la investigación de la
fusión nuclear dentro de los paı́ses latinoamericanos.

3
Capı́tulo 1

Teorı́a general de transporte

La modelación de plasmas consiste en resolver ecuaciones de movimiento que describen


el estado del plasma. Este conjunto de ecuaciones por lo general son acopladas con
las ecuaciones de Maxwell para los campos electromagnéticos ó con las ecuaciones de
Poisson para los campos electrostáticos. Existen distintos enfoques para la descripción
del plasma, los principales modelos por orden de complejidad son:

• Descripción por partı́culas

• Descripción de fluido

• Descripción cinética

• Descripción hı́brida (cinética-fluido)

• Descripción Giro-cinética

El modelo por partı́cula describe el plasma como cargas individuales moviéndose a


través de un campo eléctrico y magnético impuesto. El movimiento de cada partı́cula
está regulado por la ecuación de Lorentz.
El modelo cinético es uno de los tratamientos fundamentales para describir el plasma,
cuyos resultados derivan en una función de distribución dependiente de la posición, la ve-
locidad y el tiempo. La descripción cinética se logra mediante la resolución de la ecuación
de Boltzmann, o en el caso de requerir una correcta descripción de las interacciones de
Coulomb de largo alcance mediante la resolución de la ecuación de Vlasov que contiene
el campo electromagnético auto-consistente. También es posible utilizar la ecuación de

4
Fokker-Planck cuyas aproximaciones se han empleado para deducir términos de col-
isión operables. Las cargas y las corrientes producidas por las funciones de distribución
auto-consistentes determinan el campo electromagnético a través de las ecuaciones de
Maxwell.
Con el modelo de fluido se busca reducir las complejidades del modelo cinético, en
el se describe al plasma mediante los valores de cantidades macroscópicas. Las ecua-
ciones para estas cantidades se conocen como ecuaciones de fluido. Estas se obtienen
de los momentos respecto de la velocidad de la ecuación de Boltzmann o de Vlasov. Las
ecuaciones de fluido no están completas sin antes obtener la determinación de los co-
eficientes de transporte como: difusión, movilidad, frecuencia de colisiones promedio,
etc. Para determinar estos coeficientes de transporte es necesario asumir una función
de distribución de las velocidades adecuada para la velocidad de las partı́culas.
El modelo cinético da una descripción precisa de la fı́sica del problema, sin embargo es
más complejo que el modelo hidrodinámico y requiere mayor poder de cómputo para las
soluciones numéricas. El modelo hı́brido trata algunos componentes del sistema como
fluido y otros como partı́culas.
El modelo giro-cinético es apropiado para sistemas en presencia de un campo magnético
fuerte. En este modelo se promedian las ecuaciones cinéticas sobre la órbita circular del
giro-radio (cı́rculos de Larmor).
La descripción del plasma se obtiene de una extensión de la teorı́a cinética de gases y la
hidrodinámica de fluidos neutros para partı́culas cargadas. La descripción del fenómeno
tiene un mayor grado de complejidad que un fluido neutro, debido a los campos eléctricos
y los campos magnéticos que surgen del movimiento de las partı́culas con carga que
componen el plasma. Los campos que deben calcularse tienen una naturaleza auto-
consistente. Adicionalmente, un plasma magnetizado presenta anisotropı́a ante las per-
turbaciones.
Los campos eléctricos y magnéticos están descritos mediante las ecuaciones de Maxwell.
La mayorı́a de los cálculos asumen que el plasma se contiene en el vacı́o.
Los efectos del plasma se presentan en las ecuaciones de Maxwell mediante la densi-
dad de carga y fuentes de corriente producidas por la respuesta del plasma a los campos
eléctricos y magnéticos.

5
Carga del plasma:
X
ρq = n s qs
s

Densidad de corriente:
X
J= ns qs Vs
s

donde s es la especie, ns es la densidad, qs es la carga, y Vs es velocidad de flujo de la


partı́cula.
Para el caso en donde las corrientes en el plasma son pequeñas, (e.g. plasmas de
baja presión) y en presencia de un campo magnético estático, el modelo electrostático
apropiado es:
E = −∇φ

ρq
∇·E=
0

ρq
−∇2 φ =
0
En donde solo se requiere conocer la densidad de la carga. La descripción del plasma
busca proveer un procedimiento para calcular la densidad de carga y corriente para un
campo eléctrico E y magnético B.
La descripción termodinámica o mecánico estadı́stica de un plasma son posibles de
aplicar en casos donde el plasma se acerca a un equilibrio de colisiones de Coulomb. Sin
embargo los plasmas usualmente están lejos del equilibrio termodinámico o mecánico
estadı́stico, además la respuesta del plasma que se busca está en una escala de tiempo
menor a la relajación del plasma.

1.1 Descripción cinética del plasma

La descripción cinética del plasma incluye los efectos del movimiento de las partı́culas
cargadas que constituyen el plasma. En esta descripción se puede considerar a las
partı́culas como cargas puntuales, debido a que los efectos de la mecánica cuántica son
en su mayorı́a despreciables en los plasmas. Por lo tanto la distribución espacial de

6
una sola partı́cula moviéndose a lo largo de la trayectoria x(t) puede representarse por
la función delta:
δ [x − x (t)] = δ [x − x (t)] δ [y − y (t)] δ [z − z (t)]

Similarmente, la distribución espacial de velocidades de la partı́cula mientras esta se


desplaza por la trayectoria v(t) es:
δ [v − v (t)]

Aquı́ x y v son coordenadas fijas Eulerianas de un espacio-fase de seis dimensiones:

(x, y, z, vx , vy , vz )

x(t) y v(t), son coordenadas Lagrangianas que se mueven con la partı́cula.


N
X
f m (x, v, t) = δ [x − x (t)] δ [v − v (t)]
i=1

Esta es la función de distribución microscópica. Las unidades de la función de dis-


tribución es el recı́proco del volumen en el espacio-fase de seis dimensiones. El ı́ndice
m indica que la distribución es microscópica. d3 xd3 vf , es el número de partı́culas en un
elemento diferencial de volumen del espacio-fase de seis dimensiones.
La función de distribución está normalizada de tal forma que la integral sobre el es-
pacio de velocidades proporciona la densidad de las partı́culas.
Z N
m
X part
n (x, t) ≡ d3 vf m (x, v, t) = δ [x − x (t)] [ ].
m3
i=1

Al igual que la distribución f m , esta definición microscópica de la densidad forma un


pico en la posición instantánea de la partı́cula x = xi (t). La integración de la densidad
sobre el volumen V del plasma proporciona el número total de partı́culas de la especie
evaluada en el plasma. Z
d3 xn (x, t) = N
V

La trayectoria de las partı́culas xi (t), vi (t) para cada una de las partı́culas se obtiene
de las ecuaciones de movimiento en los campos magnéticos y eléctricos microscópicos,
Em y Bm .
dvi
m = q [Em (xi , t) + vi × Bm (xi , t)]
dt
7
dxi
= vi
dt

i = 1, 2, ..., N

El campo magnético y eléctrico microscópico se obtiene de las ecuaciones de Maxwell


en el espacio libre.
ρm
q
∇ · Em =
0

−∂B m
∇ × Em =
∂t

∇ · Bm = 0

∂Em
 
m m
∇×B = µ0 J + 0
∂t
Las fuentes microscópicas de carga y corriente se obtienen de los dos primeros mo-
mentos respecto de la función de distribución sobre el espacio de velocidades y realizando
la suma respecto de todas las especies.

X Z X N
X
ρm
q (x, t) ≡ qs d 3
vfsm (x, v, t) = qs δ [x − xi (t)]
s s i=1

X Z X N
X
m 3
J (x, t) ≡ qs d vvfsm (x, v, t) = qs vi (t)δ [x − xi (t)]
s s i=1

Las ecuaciones de la función de distribución, distribución de densidades, densidad


de carga y de corriente, junto con las condiciones iniciales y una relación termodinámica
para el calor producido, proporcionan una descripción microscópica completa del plasma.
Ellas describen el movimiento exacto de todas las partı́culas con carga en el plasma, y
la consecuente densidad de carga y corriente, los campos eléctricos y magnéticos que
generan, y el efecto de estos campos microscópicos en el movimiento de la partı́cula. To-
das estas ecuaciones deben calcularse simultáneamente y deben ser auto-consis tentes.
Integrar las N ecuaciones de movimiento sobre el tiempo, para obtener la descripción

8
evolutiva del plasma, sin embargo un plasma tı́pico contiene N > 1012 part, por lo que
esto serı́a algo poco práctico de realizar. Es necesario desarrollar un esquema donde
dichas propiedades microscópicas puedan ser promediadas en un conjunto de ecua-
ciones cuya solución pueda emplearse para obtener propiedades promedio, fı́sicamente
medibles como la densidad y la temperatura.
Para desarrollar el procedimiento bajo el cual se promedian las propiedades, es con-
veniente tener una sola ecuación evolutiva para la distribución microscópica f m en lugar
de tener que lidiar con un gran número de ecuaciones de movimiento. Dicha ecuación
puede obtenerse al calcular la derivada total del tiempo de:
N
X
f m (x, v, t) = δ [x − xi (t)] δ[v − vi (t)]
i=1

 N
df m

∂ dx ∂ dv ∂ X
≡ + · + · δ[x − xi (t)]δ [v − vi (t)]
dt ∂t dt ∂x dt ∂v
i=1
N  
X ∂ dxi ∂ dvi ∂
= + · + · δ[x − xi (t)]δ [v − vi (t)]
∂t dt ∂x dt ∂v
i=1
N  
X dxi ∂ dvi ∂ dxi ∂ dvi ∂
= − · − · + · + · δ [x − xi (t)] δ [v − vi (t)]
dt ∂x dt ∂v dt ∂x dt ∂v
i=1
=0

En donde:
xδ (x − xi ) = xi δ (x − xi )

vδ(v − vi ) = vi δ(v − vi )

∂ dxi ∂
δ [x − xi (t)] = − · δ [x − xi (t)]
∂t dt ∂x

∂ dvi ∂
δ [v − vi (t)] = − · δ [v − vi (t)]
∂t dt ∂v

9
Utilizando las funciones delta para cambiar la dependencia de las derivadas parciales
respecto de xi , vi a aquellas con respecto a x, v, y escribir df m /dt = 0, como:

df m ∂f m dx ∂f m dv ∂f m
≡ + · + ·
dt ∂t dt ∂x dt ∂v
∂f m ∂f m q ∂f m
= +v· + [Em (x, t) + v × Bm (x , t)] · =0
∂t ∂x m ∂v

Esta es la ecuación de Klimontovich. Matemáticamente incorpora las N ecuaciones


del movimiento de partı́culas en una sola ecuación ya que las caracterı́sticas matemáticas
de la ecuación diferencial parcial de primer orden en las siete variables continuas e in-
dependientes x, v, t, son:
dx
=v
dt

dv q
= [Em (x, t) + v × Bm (x , t)]
dt m
Esta ecuación se reduce a:

dvi
m = q [Em (xi , t) + vi × Bm (xi , t)]
dt

En el punto de posición de la partı́cula:

x → xi

v → vi

La ecuación diferencial parcial de primer orden avanza posiciones en el espacio-fase


de seis dimensiones x, v a lo largo de trayectorias descritas por la ecuación de movimiento
de la partı́cula única, independientemente de si existe una partı́cula en ese punto fase
(x, v).
Esta forma de Klimintovich de las ecuaciones de movimiento son las que serán pro-
mediadas para obtener la descripción cinética del plasma.
Para promediar la ecuación microscópica, definimos N6D en un elemento diferencial
de volumen del espacio-fase de seis dimensiones.

∆V ≡ ∆x∆y∆z

10
Y el volumen del espacio de velocidades:

∆Vv ≡ ∆vx ∆vy ∆z

Z Z
N6D ≡ d3 x d3 vf m
∆V ∆Vv
Es necesario considerar un elemento de volumen mayor al espacio promedio que ex-
iste entre las partı́culas en el plasma, con el fin de que exista una cantidad de partı́culas
en el elemento de volumen suficiente para evitar fluctuaciones estadı́sticas, sin em-
bargo debe cuidarse que no sean tantas que impliquen la variación significativa de las
propiedades macroscópicas como la densidad promedio: ∆x  n1/3 , en el espacio fı́sico,
vT
∆vx  1/3 , en el espacio de velocidades.
(nλ3D )
Para aplicaciones de plasmas el tamaño del elemento de volumen debe ser menor o
del orden de la longitud de Debye λD para la cual N6D ∼ (nλ3D )2  ∆.
El criterio dentro del cual deben encontrarse las dimensiones del elemento de volumen
está determinado por:
1
n− 3 < 4x < λD

La función de distribución promedio hf m i se define como el número de partı́culas


contenidas en el elemento de volumen del espacio-fase de seis dimensiones dividido entre
el volumen del elemento.

N6D
hf m (x, v, t)i ≡ lim
n−1/3 <∆x<λD ∆V ∆Vv
d3 x d3 vf m
R R
∆V ∆Vv
= lim R R
n−1/3 <∆x<λD ∆V d3 x ∆Vv d3 v

Las unidades de la función de distribución promedio son el número de partı́culas por


unidad de volumen del espacio-fase de seis dimensiones part. · s3 /m6 .
La desviación de la distribución microscópica completa f m de su promedio, la cual
por definición debe tener cero de promedio, estará expresada como δf m .

δf m ≡ f m − hf m i

hδf m i = 0
11
Esta es una distribución de partı́culas discreta. La función de distribución promedio hf m i
representa las propiedades suavizadas de las especies del plasma para una 4x ≥ λD , la
función de distribución microscópica δf m representa los efectos de “partı́cula discreta”
de las partı́culas con carga de forma individual para el intervalo:
1
n− 3 < 4x < λD

La distribución microscópica f m es “picuda” debido a que representa las partı́culas


puntuales por medio de funciones delta. Por otro lado la función de distribución promedio
hf m i indica el número de partı́culas sobre escalas de longitud que son grandes en com-
paración con el espacio entre partı́culas promedio. Finalmente la función δf m también
presenta un comportamiento de pulsos instantáneos, sin embargo tiene una base de -
hf m (x)i, por lo que su promedio se desvanece.
Además de las funciones de distribución, debemos separar el campo eléctrico y magnético,
y las densidades de carga y corriente en su parte de partı́cula discreta y la que tiene un
comportamiento suavizado.
E = hEm i + δEm

B = hBm i + δBm

ρm

m m
q = ρq + δρq

J = hJm i + δJm

Al sustituir estas expresiones en la ecuación de Klimontovich y promediando la ecuación


resultante mediante la definición previamente establecida, obtenemos la ecuación fun-
damental cinética del plasma:
∂ hf m i ∂ hf m i q ∂ hf m i
+v· + [hEm i + v × hBm i] ·
∂t  ∂x m ∂v
m

−q ∂δf
= [δEm + v × δBm ] ·
m ∂v
Los términos de la izquierda describen la evolución de la función de distribución
promedio suavizada en respuesta a los campos eléctricos y magnéticos, promediados y
12
suavizados presentes en el plasma. El término de la derecha representa las correlaciones
entre dos partı́culas discretas y con carga con una separación de la distancia de Debye
entre ellas, esto es los pequeños efectos de las colisiones de Coulomb que tienen sobre
la función hf m i. Al promediar las ecuaciones microscópicas de Maxwell y las fuentes de
densidad de carga y corriente, obtenemos ecuaciones promediadas y suavizadas que no
tienen términos extra de correlación, como el término derecho de la ecuación cinética.
La descripción por partı́cula es apropiada para plasmas de baja densidad contenidos
en campos magnéticos de gran intensidad y es de gran utilidad para comprender los
procesos fı́sicos involucrados.

1.2 Ecuaciones cinéticas del plasma

En esta sección vamos a escribir los parámetros microscópicos que se han promediado
y suavizado sin el ı́ndice m y sin corchetes de promedio. hf m i → f (x, v, t), hEm i → E,
hBm i → B, ρm → ρq , hJm i → J y donde C (f ), es el operador de colisiones de Coulomb


q

sobre la función f .
Tomando en cuenta estas especificaciones, podemos escribir la ecuación cinética
como:
∂ hf m i ∂ hf m i q ∂ hf m i
+v· + [hEm i + v × hBm i] ·
∂t  ∂x m ∂v
∂δf m

−q m m
= [δE + v × δB ] ·
m ∂v

df ∂f ∂f q ∂f
= +v· + [E + v × B] · = C (f )
dt ∂t ∂x m ∂v

f = f (x, v, t)

Además de la ecuación fundamental cinética del plasma necesitamos las ecuaciones


promedio de Maxwell y las densidades de carga y corriente.
ρq
∇·E=
0

∂B
∇×E=−
∂t
13
∇·B=0

 
∂E
∇ × B = µ0 J + 0
∂t

X Z
ρq (x, t) ≡ qs d3 vfs (x, v, t)
s

X Z
J (x, t) ≡ qs d3 vvfs (x, v, t)
s

Este es el conjunto de ecuaciones fundamentales que proveen una descripción cinética


completa del plasma. Todas las cantidades contenidas en ellas se encuentran linealizadas
y promediadas. Los efectos discretos de las partı́culas en el plasma están contenidos
en el operador de colisión de coulomb, mediante el procedimiento de promediado de
las funciones del plasma implı́citamente se asume que los efectos del comportamiento
discreto de las partı́culas no se extienden mas allá de la distancia de Debye λD .
Para plasmas de baja presión donde las corrientes del plasma son despreciables y
el campo magnético, si está presente, es constante en el tiempo es posible emplear una
aproximación electroestática para el campo eléctrico:

E = −∇φ

Con lo que las ecuaciones se reducen a:

∂f ∂f q ∂f
+v· + [−∇φ + v × B] · = C (f )
∂t ∂x m ∂v

ρq
−∇2 φ =
0

X Z
ρq = qs d3 vfs (x, v, t)
s

Este conjunto de ecuaciones proporciona una descripción, electrostática y cinética del


plasma completa.

14
En un plasma magnetizado con un radio de giro pequeño en comparación con las
escalas de longitud perpendiculares al gradiente ρ∇⊥  1, y procesos pequeños en com-
paración con la frecuencia de giro ∂/∂t  ωc , es conveniente cambiar las variables in-
dependientes del espacio-fase de x y v a las coordenadas de los centros guı́a: x, εg , µ,
la tercera variable del espacio-velocidad, es el ángulo ϕ del movimiento de giro que se
promedia para obtener las ecuaciones de movimiento de los centros guı́a.
Retomando la función de las ecuaciones del movimiento de partı́culas para obtener las
ecuaciones de Klimintovich, encontramos que en términos de las coordenadas de centros
guı́a, la ecuación cinética del plasma se transforma en:
   
∂f dxg dµ ∂f dεg ∂f
+ · ∇f + + = C (f )
∂t dt dt ∂µ dt ∂εg

El promedio de giro de la derivada temporal del momento magnético y ∂f /∂µ son


ambos pequeños en la pequeña expansión del radio de giro, por lo tanto su producto
puede ser despreciado en la ecuación cinética de primer orden del plasma.
La derivada temporal de la energı́a puede calcularse al orden más bajo, mientras se
desprecie la velocidad de deriva vD , utilizando la ecuación de centros guı́a, escribiendo
el campo eléctrico en su forma general:

∂A
E = −∇Φ −
∂t

d ∂ dxg ∂
= + ·
dt ∂t dt ∂x

∼ ∂
= + vk ∇k
∂t
Tenemos:
mvk2
!
dεg d dΦ dB
= +q +µ
dt dt 2 dt dt

∼ ∂Φ ∂B ∂A
=q +µ − qvk b̂ ·
∂t ∂t ∂t
Después de promediar la ecuación cinética del plasma sobre el ángulo ϕ, la ecuación
cinética del plasma para el promedio de giro, con la función de distribución de centros

15
guı́a fg puede escribirse en términos de las coordenadas momento y energı́a como:

∂fg dεg ∂fg


+ vk b̂ · ∇fg + vD⊥ · ∇fg + = hC (fg )iϕ
∂t dt ∂εg

fg = f (xg , εg , µ, t)

Esta es la ecuación cinética de deriva, en ella el operador de colisión se encuentra


promediado sobre la fase de giro y el gradiente espacial se toma con: (εg , µ, t) constantes.
Esto es:

∇≡ |ε ,µ,t
∂x g
Esta ecuación cinética de deriva a orden uno es suficiente para la mayorı́a de las apli-
caciones. Sin embargo al igual que las orbitas de centros guı́a en las que está basada,
es incorrecta en ecuaciones de segundo orden de la expansión del radio de giro, conse-
cuencia de ello es que no puede expresarse en forma conservativa. Existen ecuaciones
giro-cinéticas más generales y de mayor precisión.
Para muchos procesos del plasma estaremos interesados en escalas de tiempo cortas
durante las cuales los efectos de las colisiones de Coulomb son despreciables. Para estos
casos podemos considerar la ecuación cinética como:

df ∂f ∂f q ∂f
= +v· + [E + v × B] · =0
dt ∂t ∂x m ∂v

Esta es la ecuación de Vlasov, y también es conocida como la ecuación cinética sin col-
isiones. Dado que la ecuación de Vlasov no tiene efectos por correlaciones de partı́culas
discretas, es reversible en el tiempo y sus soluciones siguen las orbitas simples de partı́culas
sin colisión en el espacio-fase de seis dimensiones. Por lo tanto las funciones de dis-
tribución solución conservan la entropı́a, esto es no hay una relajación irreversible de
irregularidades en la función de distribución, y como las órbitas de la partı́cula, incom-
presible en el espacio-fase de seis dimensiones.

1.3 Ecuación de continuidad

La ecuación de continuidad describe el transporte de una cantidad fı́sica conservada.


Esta expresión matemática de la ley de conservación puede darse de forma integral en
16
términos de un flujo integral o diferencial en términos del operador de divergencia apli-
cado en un punto. Diversos fenómenos fı́sicos pueden describirse mediante el empleo
de ecuaciones de continuidad cuya estructura depende de la cantidad a describir. Las
ecuaciones de conservación son una expresión local de las leyes de conservación. Las
ecuaciones de continuidad incluyen términos fuentes y sumideros.
La ecuación de continuidad es aplicable en cantidades fı́sicas capaces de fluir o mo-
verse. La forma en que se desplaza dicha cantidad está descrita por su flujo. El flujo de
la cantidad descrita es un campo vectorial, denotado por j.
La forma integral de la ecuación de continuidad sostiene que:

• La magnitud de la cantidad fı́sica descrita se incrementa regionalmente cuando flujo


adicional de la cantidad se agrega a través de la frontera de la región, y disminuye
cuando el flujo es hacia fuera.

• La magnitud de la cantidad fı́sica descrita se incrementa regionalmente cuando la


cantidad se crea dentro de la región, y disminuye cuando se consume dentro de la
región.

• Más allá de estos dos procesos, no existe alguna otra forma para que la magnitud
de la cantidad se modifique dentro de la región.

La expresión matemática de la forma integral está definida como:


I
dq X
+ j · dS =
dt s
H
donde S, es una superficie cerrada imaginaria que encierra un volumen V , s dS, denota
la integral de superficie sobre la superficie cerrada S, q es el valor total de la cantidad
P
fı́sica de flujo, j es el flujo de q, t es el tiempo, y es la tasa neta de acumulación.
La expresión matemática de la forma diferencial está definida mediante el teorema de
la divergencia.
∂q
+∇·j=σ
∂t
donde ∇· denota la divergencia, q es la cantidad total de q por unidad de volumen, j es el
flujo de q, t es el tiempo, y σ es la acumulación de q por unidad de volumen y por unidad
de tiempo. Esta ecuación general debe utilizarse para deducir cualquier ecuación de
17
continuidad, desde la sencilla ecuación de continuidad de volumen hasta la de Navier-
Stokes.

1.4 Relajamiento del plasma

La integral de colisión de Landau posee una importante propiedad: se reduce a cero


cuando las partı́culas tienen una distribución de Maxwell en equilibrio estadı́stico. Para
funciones de distribución diferentes a la Maxweliana, la integral de colisión no desa-
parece. Ligada a esta propiedad del plasma, la solución de la ecuación cinética para un
plasma no magnetizado, es decir donde no existe un campo externo, se busca alrededor
de un comportamiento asintótico (t → ∞).
La solución de la ecuación cinética para un plasma no magnetizado, es decir donde no
existe un campo externo. La solución a este problema nos indica que independientemente
del valor inicial de la función de distribución F (p, 0) la solución a la ecuación cinética
que contiene la función de distribución de una sola partı́cula, para el caso en donde
E (e) = B (e) = 0 tiende asintóticamente conforme (t → ∞), hacia una distribución de
Maxwell; limt→∞ F (p, t) = FM (p).
   i ∂F (1) 
∂F (1) ∂F (1) 1h
+ v1 · +e E(e) (r1 ) + v1 × B(e) (r1 ) ·
∂t ∂r1 c ∂p1
 
∂F (1)
=
∂t c

Esta ecuación es válida mientras la función F (1) tenga una variación suficientemente
lenta en el tiempo (τF  τG ) y los campos externos sean suficientemente débiles. La
cantidad (∂F (1)/∂t)c denota la integral de colisión.
Si omitimos la integral de colisión en la ecuación cinética el comportamiento asintótico
no se cumple. Se puede decir que el proceso por el cual la función de distribución se
aproxima a una distribución de Maxwell, ocurre gracias a las colisiones binarias entre
las partı́culas del plasma. Este proceso es la relajación del plasma. Una postura más
general es que los efectos de las correlaciones de mayor orden son responsables por la
relajación en el plasma.
El tiempo en que la función F (p, t) se aproxima a una función de Maxwell FM (p) es el

18
tiempo de relajación.
3
l T2 1
τ∼ ∼ 4 m2
v̄ e nL
donde l, es la trayectoria libre promedio, v́ es la velocidad térmica promedio de las
partı́culas, L es el logaritmo de Coulomb.
Debido a que el tiempo de relajación depende de la masa de las partı́culas, este es
diferente para distintos tipos de partı́culas. Para las condiciones de un plasma termonu-
clear con temperaturas altas y densidad bajas, tanto la trayectoria libre promedio como
el tiempo de relajación pueden tomar grandes valores. l ∼ 106 cm, τe ∼ 10−4 s, τp ∼ 10−2 s.
El inverso de τ determina la frecuencia de colisión, definido como el promedio de
colisiones de una partı́cula por unidad de tiempo.
Debido a que la integral de colisión desaparece para una distribución de Maxwell,
para estimar la integral de colisión podemos asumir que:
 
∂F (p) 1
∼ − [F (p, t)] − [FM (p)]
∂t c τ

Con este resultado la ecuación cinética para un plasma uniforme y sin campo magnético
externo obtiene la forma:

∂F (p, t) 1
= − [F (p, t)] − [FM (p)] ,
∂t τ

cuya solución analı́tica es

t
F (p, t) = FM (p) + [F (p, 0) − FM (p)] e− τ .

Esta fórmula proporciona una descripción muy general del proceso de relajación. Esta
descripción emplea una descripción uniforme para las velocidades del plasma cuando en
realidad existen regiones de baja y alta velocidad en donde τ es mayor. El modelo más
simple de la integral de colisión no toma en cuenta esta variación, sin embargo el empleo
de una función más precisa conduce a una función integral-diferencial cuya solución
solo puede ser obtenida de forma numérica.
Debido a que el tiempo de relajación es proporcional a la raı́z cuadrada de la masa de
la partı́cula y la trayectoria media libre es independiente de la masa, las partı́culas lig-
eras como los electrones tienen un tiempo de relajación más corto que las partı́culas con

19
mayor masa. Este fenómeno propicia que el intercambio de energı́a entre partı́culas de
la misma especie se realice a una mayor velocidad que el que se requiere entre partı́culas
de especies diferentes. Por lo tanto el equilibrio entre iones y electrones caracterizados
por la misma temperatura será establecido solo después de que se establezcan distribu-
ciones de equilibrio termodinámico tanto para los iones y los electrones con diferentes
temperaturas. La distribución de los electrones se establecerá antes de que la de los
iones.

1.5 Descripción hidrodinámica de un plasma

Se puede describir adecuadamente el proceso oscilatorio de alta frecuencia en un plasma


quasi-estacionario empleando una ecuación con un campo auto-consistente sin una in-
tegral de colisión y en conjunto con las ecuaciones de Maxwell para el campo electro-
magnético en el vacı́o. Podemos emplear para el estudio del proceso oscilatorio en un
plasma la ecuación cinética con la función de distribución de una sola partı́cula, un
campo auto-consistente, y una integral de colisión junto con las ecuaciones de Maxwell.

ωτ  1

donde τ es el tiempo promedio entre colisiones binarias de partı́culas y ω es la frecuencia


de las oscilaciones.
Sin embargo dicho procedimiento es complicado. En cambio si consideramos ωτ  1,
como la condición para el tiempo de relajación significa que el lı́mite cuando τ → 0, condi-
ciona que l → 0, y l ∼ v̄τ . Este tipo de lı́mite implica que el plasma pueda ser considerado
como un ente continuo y que podemos despreciar su estructura molecular. Por ello
no es necesario el empleo de conceptos como: funciones de distribución de partı́culas,
funciones de correlación, y es suficiente utilizar una descripción hidrodinámica más
compacta en donde los conceptos básicos son la velocidad hidrodinámica del medio,
densidad, y presión. Todos estos parámetros macroscópicos, pueden variar en función
del tiempo y espacio, sin embargo estas variaciones son lentas. Si asignamos ω y L,
como la frecuencia y la longitud caracterı́stica para cambios apreciables en una escala
macroscópica, las siguientes condiciones se satisfacen:

ωτ  1
20
L1

Es factible utilizar una descripción hidrodinámica si el sistema es capaz de establecer


una distribución local de partı́culas en estado quasi-estacionario en un breve periodo
de tiempo. Los parámetros macroscópicos determinan estas distribuciones y varı́an de
acuerdo a las leyes de la hidrodinámica. Se pueden formular estas leyes a partir de la
descripción cinética.
Existen tres diferentes tiempos de relajación en el plasma: τee , τii , y τie . Cuyo equilibrio
se establece en el mismo orden. Estos tiempos de relajación están relacionados mediante:
r
mi
τii ∼ τee
me
mi
τei ∼ τee
me

Las distribuciones cuasi-estacionarias tanto de iones como de electrones están car-


acterizadas por diferentes valores de temperatura y diferentes valores promedio de ve-
locidad de partı́cula (hidrodinámica) por lo que podemos considerar un plasma multi-
componente. Cuando el equilibrio entre los iones y electrones del plasma es establecido
se tiene un valor de temperatura y de velocidad hidrodinámica de partı́cula común para
ambas especies por lo que se considera al plasma como un medio hidrodinámico de
una sola componente. En los sistemas multi-componentes e incluso en los de una sola
especie, los valores de velocidad y temperatura pueden cambiar a lo largo del espacio-
tiempo, sin embargo si se cumple:
ωτ  1

L1

Se puede considerar que estos cambios son suficientemente lentos.


Con las condiciones planteadas, es válido sugerir que si la frecuencia ω de las oscila-
ciones satisface: ωτee  1, ωτii  1 y ωτei  1, es posible considerar al plasma como un
sistema multi-componente. Sin embargo si:

ωτei  1

21
Podemos considerar al sistema como un medio hidrodinámico simple.
La descripción hidrodinámica está ligada con el brevı́simo tiempo de relajación esti-
mado para el plasma en el caso de un plasma con una tasa de colisión elevada. Sin em-
bargo para el caso de un plasma con una baja tasa de colisiones e incluso para un plasma
“frı́o” la velocidad de la partı́cula es también la velocidad hidrodinámica del medio. Para
la formulación de las leyes hidrodinámicas de un plasma, comenzaremos considerando
al plasma como un medio hidrodinámico simple.

1.6 Ecuaciones de la Magneto-hidrodinámica

Los parámetros básicos para la descripción hidrodinámica del plasma son: la velocidad
u = u (r, t) y la densidad ρm = ρ (r, t).
En términos de estas cantidades, podemos plantear la ecuación de continuidad como:

∂ρm
+ ∇ · ρm u = 0
∂t

Y la ecuación de Navier-Stokes para la conservación del momento lineal como:


 
∂u ∂u
ρm ≡ ρm + (u · ∇) u = −∇p + fe + fv ,
∂t ∂t

donde p es la presión, fe es la densidad de la fuerza pondero-motriz

1
fe = ρe E + [j × B] ,
c

con ρe la dendidad de carga, B la inducción magnética, j la densidad de la corriente de


conducción, fv la densidad de esfuerzo cortante definida como
 
2 1
fv = η∇ u + ζ + η ∇ (∇ · u)
3

donde η y ζ son los coeficientes de viscocisdad del plasma.


A estas ecuaciones hidrodinámicas es necesario añadir las ecuaciones de Maxwell:

1 ∂B
rotE = −
c ∂t

∇·B=0

22
1 ∂E 4π
rotH = − + j
c ∂t c

∇ · E = 4πρe

H = B − 4πM

donde H es la fuerza de campo magnético y M es el vector de magnetización. La ecuación


de estado del plasma:
p = p (ρm , θ)

Y la ecuación de transferencia de calor:

ds
ρm θ =q
dt

X ∂uα j2
q≡ παβ + (x∇θ) +
∂xβ σ
α,β

donde θ es la temperatura del plasma, παβ es el tensor de esfuerzos viscosos, y x es la


conductividad térmica del plasma.
La ecuación indica la cantidad de calor producida por unidad de tiempo y volumen,
debido a la viscosidad, la conductividad térmica y la conductividad eléctrica. En un
plasma ionizado el vector de magnetización M, se desvanece por lo que los vectores H, y
B coinciden.
En el caso de velocidades y frecuencias suficientemente bajas, cuando se cumplen las
siguientes desigualdades:
ω V V
 1,  1,  1,
σ σL c
donde V es la velocidad caracterı́stica, L es la longitud caracterı́stica y ω es la frecuen-
cia caracterı́stica, que caracterizan a las cantidades macroscópicas. Es posible hacer
algunas simplificaciones:
La corriente de desplazamiento es pequeña comparada con la corriente de conducción,
en la que se puede despreciar la corriente de convección ρe u.
Es posible despreciar la fuerza eléctrica ρe E en la fuerza ponderomotriz.
23
Tomando en cuenta que,
1 ∂E E
∼ ,
c ∂t
cT
donde T ω −1 es el tiempo caracterı́stico, tenemos:



1 ∂E 4π
 σ |E|
c ∂t
c

De las ecuaciones de los rotacionales obtenemos:

L
E∼ B
cT

E
ρe ∼
L
Entonces encontramos que:
ρ u V
e
∼ 1
σE σL

Utilizando:
E ∼ |[u × B] /c|

Obtenemos:
|ρe E| V
1 ∼
c |j × B| σL
Bajo las suposiciones anteriores, las ecuaciones básicas junto con las ecuacione s de
estado forman las ecuaciones de la magneto-hidrodinámica:

∂u 1
ρm = −∇p + [j × B] + fv
∂t c

∂ρm
+ ρm u = 0
∂t

∂s
ρm =q
∂t


rotB = j
c

 
1
j = σ E + [u × B]
c
24
divB = 0

1 ∂B
rotE = −
c ∂t
Expresando el campo eléctrico en términos del campo magnético.

1 1
E = − [u × B] + j
c σ

1 c
E = − [u × B] + rotB
c 4πσ
Sustituyendo esta ecuación en el rotacional del campo eléctrico se obtiene:

∂B
= rot [u × B] − rot (vm rotB)
∂t

c2
vm =
4πσ
El campo B, está determinado por esta ecuación y la ecuación de divergencia del
campo magnético.
La cantidad νm , es la viscosidad magnética. Si la conductividad σ, es independiente
de la posición la ecuación del cambio del campo con respecto al tiempo se reduce a:

∂B
= rot [u × B] + vm ∇2 B
∂t

Cuando la conductividad es suficientemente grande es posible despreciar el último


término de la ecuación. Esto se cumple cuando:

Rσ  1

Aquı́,
LV
Rσ = ,
vm
es el número magnético de Reynolds ó también se conoce como el número de Lundquist.
Teniendo en cuenta que:

rot [u × B] = (B · ∇) u − (u · ∇) B − B(∇ · u)
25
y,  
1 ∂ρm u
divu = − − · ∇ρm ,
ρm ∂t ρm
se puede demostrar que:
 
d B ∂ B
= + (u · ∇)
dt ρm ∂t ρm
 
B
= ·∇ u
ρm

Empleando un elemento infinitesimal de una lı́nea de fluido tenemos:

d
δl = (δl · ∇) u
dt

Cuando Rν  1
Rσ  1

Las lı́neas del campo magnético se mueven junto con las partı́culas del plasma situ-
adas en las lı́neas, e indica que la disipación de la energı́a causada por la conductividad
del medio es pequeña.
LV
Rχ =
χ
Es otro número adimensional donde se relaciona la conductividad térmica del plasma
mediante:
κ
χ=
cp
donde κ es la conductividad térmica y cp es el calor especı́fico.
Cuando Rχ  1, es suficientemente grande la energı́a disipada por la conductividad
térmica es una cantidad pequeña.
Finalmente si Rν  1, implica que la energı́a disipada por la viscosidad del medio es
pequeña. Si la desigualdad  1, de estos tres números adimensionales se cumple, es
posible despreciar cualquier disipación de energı́a. Y nos encontramos ante el caso de
un medio ideal, para el cual las ecuaciones magneto-hidrodinámicas se reducen a:

∂ρm
+ div (ρm u) = 0
∂t

ds
=0
dt
26
du 1
ρm = −∇p + [rotB × B]
dt 4π

∂B
= rot [u × B]
∂t

divB = 0

Las suposiciones de la nula disipación de energı́a, y que los procesos macroscópicos


se desarrollan lentamente son necesarias para aplicar la descripción hidrodinámica, en
donde:
ds
= 0,
dt
indica que se trata de un proceso adiabático. El valor de V en el caso de un flujo esta-
cionario es la velocidad de flujo, mientras que para el movimiento de las ondas es la
velocidad de fase. En el caso ideal y no relativista Rσ  1, implica Rχ  1 y Rν  1.

27
Capı́tulo 2

Aproximación de fluido

Para muchas aplicaciones del plasma, la descripción de la aproximación de fluido de


una especie con carga puede ser suficiente. Este es generalmente el caso cuando no hay
velocidades particulares ó regiones del espacio de velocidades donde hay partı́culas con
carga con un comportamiento diferente a las partı́culas térmicas de esa especie.
Para la deducción de las ecuaciones de la aproximación de fluido, es necesario definir
los momentos de velocidad de la función de distribución que vamos a requerir. Estos re-
sultan de integrar potencias de bajo orden de la velocidad v, multiplicadas por la función
de distribución f sobre el espacio de velocidades en el marco del laboratorio.
Z
d3 vv j f, j = 0, 1, 2 . . . (2.1)

Las integrales son todas finitas por que la función de distribución debe relajarse rápidamente
con la velocidad para que los momentos fı́sicos de bajo orden sean finitos. Esto quiere
decir que no es posible tener una cantidad grande de partı́culas con energı́as arbitraria-
mente grandes ya que la energı́a en las especies será infinita o divergente.
Los momentos de velocidad que nos interesan de la función de distribución f (χ, ν, t)
son:

• Densidad (part/m3 ): Z
n≡ d3 vf.

• Velocidad de flujo (m/s): Z


1
V≡ d3 vvf.
n

28
• Temperatura (J, eV ):
mvr2 mvT2
Z
1
T ≡ d3 v f=
n 3 2

• Flujo conductivo de calor (W/m2 ):

mvr2
Z
q≡ d3 vvr ( )f.
2

• Presión (N/m2 ):
mvr2
Z
p≡ d3 v f = nT.
3

• Tensor de presión (N/m2 ):


Z
P ≡ d3 vmvr vr f = pl + π.

• Tensor de esfuerzos (N/m2 ):

vr2
Z  
3
π≡ d vm vr vr − l f.
3

Donde
vr ≡ v − V (x, t) ≡ |vr |,

es la velocidad relativa.
Por definición tenemos:

Z
d3 vvr f = n(V − V ) = 0

Todas estas propiedades de la aproximación de fluido de alguna especie con carga en


particular contenida en un plasma son por lo general funciones de la posición espacial
x y el tiempo t. En este caso, la densidad n es el promedio suavizado de la densidad
microscópica, la velocidad de flujo V , es la velocidad de flujo macroscópico de las especies
de partı́culas, la temperatura T es la energı́a promedio de la especie de partı́culas, el
flujo de calor conductivo q, es el flujo de densidad de energı́a medido desde el resto del
marco inercial de la especie de partı́cula al igual que la temperatura. La presión p es
una función escalar que representa la parte isotrópica de los esfuerzos expansivos de las
partı́culas ya que el movimiento térmico causa que se expandan isotrópicamente lejos

29
de su posición inicial. El tensor de presión P representa la presión total de la especie,
con componentes isotrópicos y anisotrópicos. Finalmente el tensor de viscosidad π es un
tensor simétrico de seis componentes, que representa las componentes anisotrópicas del
tensor de presión.
Adicionalmente necesitaremos los momentos de velocidad de menor orden del oper-
ador de colisión de Coulomb C(f ). Se tiene conservación de la densidad en colisiones
Z
0 = d3 vC(f ).

La densidad de fuerza de fricción (N/m3 ) se define como


Z
R ≡ d3 vmvC(f ).

Densidad de energı́a de intercambio (W/m3 ):

(mvr2 )
Z
Q ≡ dv C(f ).
2
De la integral respecto de las velocidades para las colisiones (momento lineal) se indica
que las colisiones de Coulomb no crean ni destruyen partı́culas con carga, por lo que esta
integral se anula.
El momento del operador de colisión de Coulomb representa la ganancia o pérdida
por unidad de volumen de las especies de partı́culas cargadas que fluyen relativamente
a otras regiones.

Re ∼
= −me ne ve (Ve − Vi )
J
= ne e
σ

Ri = −Re

El momento cuadrático del operador de colisión representa la tasa de intercambio de


energı́a por unidad de volumen debido a las colisiones de Coulomb entre 2 partı́culas con
carga de diferente especie y diferente temperatura:

Qi = 3(me /mi )ve ne (Te − Ti )


30
J2
Qe ∼
= − Qi
σ
En un plasma magnetizado, la conductividad eléctrica a lo largo del campo magnético
es el valor de Spitzer:

σk = σSp

Pero perpendicular al campo magnético es la conductividad de referencia:

σ⊥ = σ0

Por lo tanto en un plasma magnetizado:

Jk J⊥
Re = −ne e(b̂ + )
σk σ⊥
y

Jk2 2
J⊥
Qe = + − Qi
σk σ⊥

2.1 Teorı́a general de momentos

Al igual que en la teorı́a cinética de los gases neutros, las ecuaciones de la aproximación
de fluido son deducidas de tomar el momento del espacio de velocidades de una ecuación
cinética relevante, para lo cual lo mas simple es usar la forma conservativa de la ecuación
cinética del plasma.

Z
∂f ∂ ∂ q
d3 vg(v)[ + · [vf ] + · [ + [E + v × B]f ] − C(f )] = 0
∂t ∂x ∂v m
En donde g(v) es la función de velocidad relevante para el momento de fluido que se
desea describir.
Evaluando esta ecuación para g = 1 obtenemos el momento de densidad. Dado que
las coordenadas Eulerianas espaciales de la velocidad v son estacionarias y por lo tanto

31
independientes del tiempo, la derivada temporal puede intercambiarse con la integral
sobre el espacio de velocidades. Por lo tanto la primera integral queda como:
Z
∂ ∂n
d3 vf =
∂t ∂t
R 3
De la misma manera dado que d v y el operador de la derivada espacial ∂/∂x se
conmutan pueden intercambiarse en el segundo término para convertirse en:
Z
∂ ∂
· d3 vvf = · nV ≡ ∇ · nV
∂x ∂x
Dado que el integrando en el tercer término se encuentra en forma de divergencia en
el espacio de velocidad, la integral puede ser transformada en una integral de superficie
utilizando el teorema de Gauss:

Z   Z  
∂ dv dv
d3 v · f= dSv · f =0
∂v dt dt
El cual se desvanece por que deben existir exponencialmente pocas partı́culas en la
superficie limitante del espacio de velocidad.

|v| −→ ∞

Con lo que todos los momentos algebraicos de la función de distribución son finitos y
por lo tanto existen.
El momento de densidad de la ecuación cinética del plasma da la ecuación de con-
tinuidad o ecuación de densidad:

∂n
+ ∇ · nV = 0,
∂t
es decir,
∂n
= −V · ∇n − n∇ · V
∂t
de donde
dn
= −n∇ · V
dt
Aquı́ para obtener la segunda forma de la ecuación se ha empleado un vector identi-
dad, y la última forma esta en términos de la derivada total temporal en un fluido que se
mueve con velocidad de flujo V:
32
d ∂
≡ |x +V · ∇,
dt ∂t
es la derivada total respecto al tiempo de un fluido en movimiento, a veces conocida como
la derivada sustancial.
Partiendo de:
∂n
= −V · ∇n − n∇ · V,
∂t
podemos observar como en una posición fija, el incremento en la densidad de una especie
en el plasma es causada por advección de las especies con velocidad de flujo V a través
del gradiente de densidad, de una región de alta densidad hasta un región local de baja
densidad (V · ∇n < 0), o por compresión del flujo (∇ · V < 0, convergencia). Recı́procamente
la densidad local disminuye si las especies del plasma fluyen de una región local de baja
densidad hacia una región local de mayor densidad, o si el flujo se expande (diverge). La
última forma:
dn
= −n∇ · V.
dt
Muestra que en un marco de referencia en movimiento con el flujo de velocidad V
(descripción de Lagrange), solo la compresión o expansión del flujo provoca la respectiva
variación en la densidad.
La ecuación de momentos para una especie del plasma se deduce de forma similar.
Utilizando g = mv, calculando los términos y utilizando vv = (vr + V) (vr + V) para
evaluar el segundo término, encontramos:

∂(nV)
m + ∇ · (pl + π + mnVV) − nq[E + V × B] − R = 0
∂t

Para la obtención del segundo hasta el último término utilizamos identidades vecto-
riales para escribir:
          
∂ dv ∂ dv dv ∂v
v · f = · v f − f· ,
∂v dt ∂v dt dt ∂v

es decir,     
∂ dv dv
· f − f.
∂v dt dt
Dado que:
∂v
≡l
∂t
33
y
dv dv
·l=
dt dt
El término que contiene la divergencia en el espacio de velocidades nuevamente de-
saparece por la conversión a la integral de superficie. Ahora se reescribe:

d ∂
≡ |x +V · ∇
dt ∂t

Utilizando:
∂n
+ ∇ · nV = 0,
∂t
obtenemos
∂n
= −V · ∇n − n∇ · V,
∂t
de donde
dn
= −n∇ · V.
dt
Para remover la contribución de:
∂n
∂t
y
∇ · mnVV = mV (∇ · nV) + mnV · ∇V

Para obtener:
dV
mn = nq [E + V × B] − ∇p − ∇ · π + R.
dt
En donde la derivada total respecto al tiempo del fluido en movimiento es aquella definida
como:
d ∂
≡ |x +V · ∇.
dt ∂t
Los siguientes dos términos representan la fuerza por unidad de volumen sobre las
especies que resulta del tensor de presión:

P = pl + π.

Tomando en cuenta las contribuciones tanto de la presión expansiva isotrópica p, y


del estrés anisotrópico π.

34
El término R representa la densidad de la fuerza de fricción en las especies que resul-
tan de la relajación por colisiones de Coulomb del flujo V hacia las velocidades de flujo
de otras especies de partı́culas con carga en el plasma.
Para la obtención de la ecuación de energı́a para una especie del plasma, tomando el
momento de energı́a de la ecuación cinética del plasma. Empleando g = mv 2 /2 en:

Z
∂f ∂ ∂ q
d3 vg(v)[ + · [vf ] + · [ + [E + v × B]f ] − C(f )] = 0,
∂t ∂x ∂v m
y procediendo de la misma forma que para la ecuación de momentos, obtenemos:

     
∂ 3 1 2 5 1 2
nT + mnV +∇· q+ nT + mnV V+V·π −
∂t 2 2 2 2
nqV · E − Q − V · R = 0

Utilizando el producto escalar de la ecuación de momentum con V para remover el


término ∂V 2 /∂t de la ecuación y utilizando la ecuación de densidad se puede simplificar
en:  
3 ∂p 5
= −∇ · q + pV + V · ∇p − π : ∇V + Q = 0,
2 ∂t 2|
3 ∂p 5
+ p∇ · V = −∇ · q − π : ∇V + Q.
2 ∂t 2
La primera forma de la ecuación de energı́a muestra que la tasa local (Euleriana) a la
que se incrementa la energı́a interna por unidad de volumen de las especies,

3 3
nT = p,
2 2
está dada por la suma neta (divergencia) de los flujos de energı́a dentro del volumen local
debido a la conducción de calor (q).
La convección de calor, (5/2)pV − (3/2)pV es la energı́a interna transportada a lo largo
con la velocidad de flujo V más pV debido al trabajo mecánico desarrollado sobre o por
las especies durante el movimiento.
Advección por presión, desde una zona de baja presión hacia una zona local de mayor
presión, es la disipación debido al estrés por flujo-gradiente-inducido en las especies
(−π : ∇V) y el intercambio de energı́a por colisiones (Q).

35
La ecuación de energı́a usualmente está escrita en la forma de una ecuación para la
derivada temporal de la temperatura. Esta forma se obtiene utilizando la ecuación de
densidad para eliminar el término ∂n/∂t implı́cito en ∂p/∂t en:

3 ∂p 5
+ p∇ · V = −∇ · q − π : ∇V + Q,
2 ∂t 2

para obtener
3 dT
n = −nT (∇ · V) − ∇ · q − π : ∇V + Q.
2 dt
Finalmente es útil escribir las ecuaciones de energı́a en términos de la entropı́a de
colisión en vez de la presión y la temperatura.
La entropı́a (adimensional) s para f ∼
= fM es:

3
Z !  
1 T2 3 p
s≡− d vf ln f ∼
3
= ln + C = ln 5 +C
n n 2 n3
donde C es una constante. La entropı́a representa el estado de desorden de un sistema.
Matemáticamente es el logaritmo del número de estados estadı́sticamente independi-
entes que una partı́cula puede tener en un volumen relevante del espacio-fase de seis
dimensiones.
La entropı́a en el sistema crece monótonamente en el tiempo mientras las colisiones
provocan que las partı́culas se dispersen en un mayor volumen del espacio-fase de lo que
inicialmente era un estado de mayor densidad.
El incremento de la entropı́a (ds/dt > 0) en un fluido en movimiento es causado por
el flujo de calor neto que entra en el volumen, la disipación debida al estrés por flujo-
gradiente-inducido en las especies del plasma, y el intercambio de energı́a por colisiones.
Las ecuaciones de momento del fluido para las especies con carga del plasma:

dn
= −n∇ · V
dt
dV
mn = nq [E + V × B] − ∇p − ∇ · π + R
dt
3 dT
n = −nT (∇ · V) − ∇ · q − π : ∇V + Q
2 dt
Son similares a las ecuaciones de momento obtenidas de las ecuaciones cinéticas de
un gas neutro. Las diferencias clave son:

36
• La densidad de fuerza promedio nF̄ en las especies del plasma esta dada por la
densidad de la fuerza de Lorentz n[E + V × B] en lugar de la fuerza gravitacional
−m∇VG .

• Los efectos de las colisiones de Coulomb entre las diferentes especies con carga en
el plasma deriva en la aparición de fuerza de fricción (R) e intercambio de energı́a
(Q) en las ecuaciones de momento y energı́a.

Para el caso de los plasmas existe la complejidad adicional en la cuestión de la densidad


y el flujo de varias especies debe sumarse de acuerdo a las ecuaciones:

X
ρq = n s qs
s

X
J= ns qs Vs
s

Para obtener las fuentes de densidad de carga y la corriente para las ecuaciones
de Maxwell que deben de resolverse con la finalidad de obtener los campos E y B en
el plasma, que determinan la densidad de la fuerza de Lorentz para cada especie de
partı́culas del plasma.
Es importante resaltar que mientras cada momento del fluido de la ecuación cinética
es una ecuación exacta, las ecuaciones de momento del fluido representan una jerarquı́a
de ecuaciones que, sin alguna otra especificación, no forman un conjunto completo de
ecuaciones.
La ecuación de densidad, que es la ecuación de momento con el menor orden. En
principio podrı́a resolverse para la evolución de la densidad n en el tiempo, si la velocidad
de flujo de las especies V estuviese especificada. En cambio la velocidad de flujo es
determinada de la ecuación del siguiente orden, la ecuación de momento:

dV
mn = nq [E + V × B] − ∇p − ∇ · π + R
dt

Sin embargo para resolver esta ecuación para V necesitamos conocer la presión (p =
nT ) de la especie en cuestión y por lo tanto la temperatura T y el tensor de estrés π. La
temperatura es obtenida de la versión isotrópica de la ecuación de momento de siguiente

37
orden, que es la ecuación de energı́a. Pero esta ecuación depende del flux de calor q. Por
lo tanto las ecuaciones de momento, densidad y energı́a no son ecuaciones completas
pues aún no se especifica el momento de “cierre” de mayor orden en las ecuaciones, que
corresponden al flux de calor y el tensor de estrés q y π.
Para determinarlos es posible tomar momentos de orden mayor de la ecuación cinética,
algunos de los cuales no tienen una interpretación fı́sica simple y no tienen una medición
sencilla. Fı́sicamente aún los momentos del orden mas elevado dependen de detalles cada
vez mas finos de la función de distribución f , por lo tanto es posible considerarlas como
despreciables, particularmente si consideramos los efectos de las colisiones de Coulomb
al momento de suavizar caracterı́sticas de escalas finas de la función de distribución del
espacio de velocidades. Las ecuaciones de fluido deducidas hasta este momento pro-
porcionan ecuaciones evolutivas para las propiedades fı́sicas medibles mas importantes
(n, V, T ) de las especies del plasma, por lo que es deseable cerrar la jerarquı́a de ecua-
ciones en este nivel.
El procedimiento general para cerrar la jerarquı́a de las ecuaciones de fluido es la
obtención de los momentos de cierre requeridos, en ocasiones denominados relaciones
constitutivas, a partir de integrales de la función de distribución cinética f para π y q:
Z  2
3 mr
q ≡ d vvr f
2

vr2
Z  
3
π ≡ d vm vr vr − l f
3
La función de distribución debe resolverse desde una ecuación cinética que tome en
cuenta la evolución de los momentos de fluido de menor orden:

• n(x, t)

• V (x, t)

• T (x, t)

Estos son los parámetros de la función de distribución de equilibrio dinámico de Maxwell


fM de menor orden. La ecuación cinética resultante y el procedimiento para determinar
la función de distribución y los momentos de cierre son conocidos como la aproximación
de Chapman-Enskog. En esta aproximación las distorsiones cinéticas de la funcién de
38
distribución son conducidas por las fuerzas termodinámicas ∇T y ∇V , ambos gradi-
entes de los parámetros de la distribución de Maxwell de menor orden, la temperatura
se utiliza para q, y la velocidad del flujo para π. Para situaciones en donde los efec-
tos de las colisiones sean dominantes, la ecuación cinética resultante puede resolverse
asintóticamente por medio de un esquema perturbativo y los momentos de cierre q y π
representan la difusión de calor y momento inducida por las colisiones en el medio.
En un plasma dominando por las colisiones de Coulomb el flux de calor q inducido por
el gradiente de temperatura ∇T será determinado por el proceso de difusión por colisión
de caminata aleatoria.
q ≡ −κm ∇T = −nχ∇T,

y la conductividad térmica
κm (4x)2
χ≡ ∼ .
n 24t
Esto constituye el flujo de calor de Fourier, en donde ∆x es el tamaño del paso aleatorio
en el espacio que toman las partı́culas en un tiempo ∆t.
Para los procesos de colisión de Coulomb en un plasma sin magnetizar tenemos que
4x ∼ λ es la longitud de colisión y 4t ∼ 1/v es el tiempo de colisión. Por lo tanto, el
escalamiento para la conductividad térmica es:
5
2v2 T2
χ ∼ νλ = T ∝ .
v n
En un plasma magnetizado estos procesos de colisión siguen ocurriendo libremente a
lo largo del campo magnético, siempre y cuando λ∇k  1. Perpendicularmente al campo
magnético el giro de la partı́cula delimita el tamaño de paso perpendicular ∆x en relación
al radio de giro ρ. Por lo tanto, en un plasma magnetizado con colisiones tenemos que la
conducción de calor paralela se define como

qk = −nχk b̂∇k T, χk ∼ vλ2 ,

y la conducción de calor perpendicular

q⊥ = −nχ⊥ ∇⊥ T, χ⊥ ∼ vρ2 ,

donde
∇k ≡ b̂ · ∇,
39
 
∇⊥ = ∇ − b̂∇k = −b̂ × b̂ × ∇

con
B
b̂ ≡
B
La razón de difusión de calor, en la dirección perpendicular a paralela es:
 2
χ⊥  ρ  2 v
∼ ∼ 1
χk λ ωc

La cual por definición es muy pequeña para un plasma magnetizado. Por lo tanto, en
un plasma magnetizado la difusión del calor por colisiones es mucho menor a través de
las lı́neas de campo que a lo largo de ellas, tanto para iones como para electrones.
Para electrones como para iones con, aproximadamente, la misma temperatura la
frecuencia de colisión de los electrones es mayor:
 1
νe mi 2
& 43  1
νi me

Las longitudes de colisión generalmente son de mismo orden

λe ∼ λi ,

y el radio de giro de los iones es mucho mayor:


 1
ρi mi 2
& 43  1
ρe me

Por lo tanto, las colisiones a lo largo de las lı́neas de campo magnético son las que
causan que en los electrones la difusión de calor sea a una velocidad mayor que la de los
iones. En el caso perpendicular a las lı́neas de campo, la difusión de calor por parte de
los iones es el proceso dominante.
De forma similar, el tensor de estrés viscoso p causado por la difusión por colisiones
(por caminata aleatoria), en un plasma sin magnetizar y en presencia de un gradiente en
la velocidad de flujo V es
π ≡ −2µm W

donde
µm (4x)2

nm 24t
40
es el tensor de viscocidad, donde W es la forma simetrizada del gradiente de velocidad
de flujo de las especies  
1 T 1
W≡ ∇V + (∇V ) − l (∇ · V)
2 3
Al igual que para el flujo de calor, el coeficiente de difusión para un plasma con colisiones
no magnetizado, se escala como:
µm
∼ vλ2 .
nm
Y para un plasma magnetizado se tiene:

πk = −2µm
k Wk ,

µm
k
vλ2 ,
nm
π⊥ = −2µm
⊥ W⊥ ,
µm

∼ vρ2 .
nm
 
Dado que la termodinámica considera que Wk ≡ b̂ b̂ · W · b̂ b̂ y W⊥ son cantidades
tensoriales, el álgebra se vuelve algo complicada principalmente en presencia de cam-
pos magnéticos no homogéneos. Al igual que para la difusión de calor, la difusión por
colisión del momentum a lo largo de las lı́neas de campo magnético se da a mayor ve-
locidad que a través de ellas. Debido al factor de la masa en el coeficiente de viscosidad
µm , los efectos de la viscosidad de los iones son dominantes tanto paralelamente como
perpendicularmente al campo B:

µm  1
ke me 2 1
∼ .  1,
µm
ki mi 43

3
µm

me 2
⊥e
∼ . 1.3 × 10−5  1.
µm
⊥i mi
En muchos plasmas el radio de giro ρ es mucho menor a la escala longitudinal per-
pendicular empleada para los gradientes de temperatura y flujo, por lo tanto los términos
proporcionales al radio de giro son usualmente despreciados en comparación con el resto
de los demás términos. Esto es particularmente cierto para los electrones, dado que su
radio de giro es mucho menor en comparación con el de los iones. En la aproximación
del radio de giro pequeño los flujos son generalmente pequeños en comparación con su
41
velocidad térmica, por lo tanto los términos de flujo pueden despreciarse a excepción
quizá de los términos de iones.
Cuando no existe una producción significativa de entropı́a en la escala de tiempo de
interés, la entropı́a es una constante del movimiento del fluido. Entonces obtenemos la
ecuación adiabática de estado para las especies:

ds 1 d p
≡ ln ' 0 ↔ p ∝ nΓ , T ∝ nΓ−1 ,
dt Γ − 1 dt nΓ
(N +2)
donde Γ = N con N grados de libertad. Para el caso de tres dimensiones N = 3 y
Γ = 53 .
Otras ecuaciones de estado utilizadas en la fı́sica de plasmas son: la ecuación isotérmica
de estado. (Γ = 1) p ∝ n, T = cte., la ecuación de estado de especies frı́as p ' 0, T ' 0, y la
ecuación de flujo incompresible ∇ · V = 0, n = cte. (Γ → ∞).
Para la la ecuación de flujo incompresible, asumiendo ds/dt = 0 y utilizando la ecuación
de densidad, encontramos:

1 dp 1 dn 1 1 dp
= −Γ = Γ∇ · V ↔ ∇ · V =
p dt n dt Γ p dt

De la ecuación de flujo incompresible obtenemos que (Γ → ∞), independientemente


de la evolución en la presión de cada especie, la ecuación de densidad se transforma
como
dn ∂n
= + V · ∇n = −n (∇ · V) = 0
dt ∂t
Por lo tanto, la densidad es constante en el elemento de fluido en movimiento para
un flujo incompresible (caso Lagrangiano), sin embargo la densidad cambia en el tiempo
bajo un enfoque Euleriano debido a la advección del fluido hacia regiones del espacio
con diferentes densidades. Dado que la presión no está determinada por la ecuación de
estado para un flujo incompresible, aún se necesita obtener por separado en este modelo.
Cuando se utiliza alguna de las ecuaciones de estado, esta proporciona un relación
de cierre relacionando la presión p o la temperatura T a la densidad n; por lo tanto, esta
ecuación reemplaza la ecuación de energı́a o entropı́a para cada especie.
Cuando se utiliza la ecuación de estado para flujo incompresible, esta solo actúa como
una condición de restricción en el flujo. Para este caso, debe resolverse adicionalmente

42
una ecuación de energı́a o entropı́a relevante para obtener la evolución de la presión o la
temperatura de la especie en términos de la densidad n y otras variables.

43
Capı́tulo 3

Transporte neoclásico y el código


ASTRA

Los códigos de transporte han sido una herramienta esencial en modelar el compor-
tamiento de fluidos de distinta naturaleza. En la investigación de la fusión nuclear con-
trolada, la implementación y dominio de estos códigos es de gram relevancia dado la com-
plejidad del comportamiento del plasma, y el elevado costo de este tipo de experimentos.
Las primeras simulaciones de transporte y códigos de transporte para Tokamaks comen-
zaron a desarrollarse a finales de los años 60’s. El desarrollo posterior de estos códigos se
ha visto influenciado por los avances en la teorı́a de plasmas de fusión, los experimentos
en Tokamaks, y en el desarrollo de nueva tecnologı́a y técnicas de cómputo.
A pesar de los grandes avances logrados en este campo, los procesos de transporte
en este tipo de sistemas aún se encuentran lejos una comprensión total del fenómeno.
Pues aunque se han desarrollado constantemente modelos teóricos y la capacidad de
predicción mejora con cada uno, hoy en dı́a la mayor parte de los modelos empleados
están basados en correlaciones deducidas del escalamiento de resultados experimen-
tales ası́ como de aproximaciones por prueba y error. En base a estas cuestiones es
razonable suponer que se debe realizar una gran cantidad de cálculos antes de lograr
una concordancia aceptable entre los resultados teóricos y los experimentales.
El código ASTRA resuelve ecuaciones de transporte acopladas, unidimensionales, y
dependientes del tiempo, para la descripción de partı́culas, calor, y corriente. También
resuelve sistemas bidimensionales para el equilibrio MHD auto-consistente en la ge-
ometrı́a real del reactor. ASTRA es un código capaz de crear modelos numéricos para

44
la interpretación y predicción de modelos de transporte, análisis de estabilidad, y proce-
samiento de datos experimentales. El código ASTRA es personalizable a las necesidades
del modelo, lo que lo convierte en un código muy flexible. Se pueden adaptar subruti-
nas, módulos, mÃľtodos de resolución, etc. para representar diferentes procesos fı́sicos,
simular equipo auxiliar, y procesamiento de datos. ASTRA tiene la capacidad de generar
modelos interactivos para observar la evolución temporal de los parámetros del plasma
durante la ejecución, la presentación de los datos de salida, y el control de parámetros
para modificar el curso de la simulación en tiempo real.

3.1 Sistema de coordenadas, promedios de superficie y flujos

Para propósitos descriptivos de la geometrı́a en el codigo ASTRA, se considera el sistema


magnético con simetrı́a toroidal. En la vida real, pequeñas violaciones a esta simetrı́a
siempre están presentes y deben ser consideradas para las propiedades de confinamiento
del plasma, sin embargo son despreciadas en el equilibrio.
Para la descripción de la geometrı́a de un Tokamak se emplean dos sistemas de coor-
denadas:

• Un sistema cilı́ndrico, con el eje polar coincidiendo con el eje mayor del toroide.
{r, ϕ, z}

• El segundo queda descrito por las variables {a, θ, ζ}y es un sistema relacionado a
la geometrı́a magnética del Tokamak, donde: a es una variable “radial” y es una
etiqueta arbitraria para una superficie de flux magnético. θes el ángulo poloidal y
ζ es el ángulo toroidal y esta relacionado al primer sistema coordenado mediante
ζ = −ϕ

La superficie de flujo está definida mediante a = Sa , introduciendo una integral de


superficie sobre el interior V de esta superficie de flux:
Z a Z Z a Z 2π Z 2π

Z
dSa
V = dV = da = da dζ gdθ
V 0 Sa |∇a| 0 0 0
Z a Z 2π

= 2π da gdθ, (3.1)
0 0

45
donde g es el determinante del tensor métrico
 
D (x, y, z)
g= = (∇a∇θ∇ζ)−2
D (a, θ, ζ)

y V es el volumen de un tubo de plasma


La superficie de flujo promedio de una función arbitraria f (r) está definida como:
Z a Z 2π Z 2π

Z
∂ ∂
hf i = f dV = da dζ gdθ
∂V V ∂V 0 0 0
∂a 2π √
Z
= 2π gf dθ (3.2)
∂V 0

Las derivadas parciales en la ecuación indican que todas las funciones de superficie,
es decir las funciones de la coordenada espacial a, también son funciones dependientes
del tiempo t. La expresión 3.2 se puede representar de forma equivalente como
Z Z I
∂ dSa ∂ψ dlθ
hf i = f dV = f = f
∂V V Sa |∇V | ∂V Bpol
I I
dlθ dlθ
= f / (3.3)
Bpol Bpol

Para cualquier función axisimétrica f = f (a, θ)se tiene:


Z

hB·∇f i = hdiv (f B)i = div (f B) dV
∂V V
I

= f B · ds = 0 (3.4)
∂V Sa

Empleando la ecuación (3.3), obtenemos el área de la superficie magnética:


I Z
∂ ∂V
Sa = dSa = |∇V | dV = h|∇V |i = h|∇a|i (3.5)
Sa ∂V V ∂a

Algunas propiedades importantes del procedimiento para promediar son:


Z 2π
∂V √
hF (a) f (a, θ)i = F (a) hf i , = 2π gdθ,
∂a 0
∂ ∂Γ
hdivgi = hg · ∇V i = (3.6)
∂V ∂V

3.2 Ecuaciones de transporte

El código ASTRA como muchos otros códigos empleados para la simulación de reactores
Tokamak, se compone de un sistema de ecuaciones difusivas en una dimensión para la
46
densidad y la temperatura de los diferentes componentes del plasma, una ecuación de
equilibrio en dos dimensiones, y una variedad de módulos adicionales para la descripción
del calentamiento, la inducción de corriente y para otros procesos no difusivos en el
plasma.
En esta sección se presenta el conjunto de ecuaciones de transporte implementadas
en el código ASTRA, que permiten la variación del campo toroidal B0 y en consecuencia,
la compresión adiabática en un menor radio.
El conjunto básico de las ecuaciones de transporte del código ASTRA incluye ecua-
ciones para la densidad electrónica ne , la temperatura electrónica Te , la temperatura
iónica Ti , y el flujo poloidal ψ.

!
1 ∂ Ḃ0 ∂ 1 ∂
ρ V 0 ne + 0 Γe = Se


V0 ∂t 2B0 ∂ρ V ∂ρ
!  
3 − 5 ∂ Ḃ0 ∂ h 5 i 1 ∂ 5
V0 3 − ρ V 0 3 ne Te + 0 qe + Te Γe = Pe ,
2 ∂t 2B0 ∂ρ V ∂ρ 2
!  
3 0
 − 5 ∂ Ḃ0 ∂ h 5
0 3
i 1 ∂ 5
V 3
− ρ V ni Ti + 0 qi + Ti Γi = Pi , (3.7)
2 ∂t 2B0 ∂ρ V ∂ρ 2
!
J 2 R0 ∂ G2 ∂ψ V0
 
∂ψ Ḃ0 ∂ψ
σk − = − (jBS + jCD ) .
∂t 2B0 ∂ρ µ0 ρ ∂ρ J ∂ρ 2πρ

Si la densidad principal de los iones ni y su flujo Γi no se encuentran explı́citamente


definidos, se asume que ni = ne /Zi y Γi = Γe /Zi , este tipo de aproximación no toma en
cuenta las impurezas. Todos los flujos definidos en el conjunto de ecs.(3.7), estos son; el
flujo electrónico Γe , el flujo de calor electrónico qe , y el fluo de calor de los iones qi , son
considerados como flujos totales a través de la superficie de flujo ρ = cte. Con la intención
de cerrar el sistema de ecuaciones, los flujos deben expresarse en términos de las fuerzas
termodinámicas tomadas como derivadas con respecto de ρ. En una simulación, los
flujos ası́ como la conductividad σk y la densidad de corriente bootstrap, puede tomarse
de alguna teorı́a de transporte o de las observaciones experimentales en un Tokamak.
Independientemente del origen de la función utilizada para la representación de los flujos,
se asume en el código ASTRA que la matriz de transporte tiene la siguiente forma:

47
 Γe
    1 ∂ne 
ne Dn De Di DE ne ∂ρ
qe 1 ∂Te
   χen χe χe χe   Te ∂ρ

ne Te
 = −V 0 G1  i E ·
. (3.8)
  
qi i i i   1 ∂Ti
χ χ χ χ

i
ni Ti
 n e Ti ∂ρ
  E  
µ0 V 0 G1 jBBS Cn Ce Ci 0 Ek
p Bp

Todos los coeficientes del menor superior izquierdo de tercer orden en la matriz de
transporte tiene dimensiones m2 /s , mientras la columna derecha y la fila del fondo son
 

adimensionales. El campo magnético poloidal Bp está dado por:

def 1 ∂ψ
Bp =
2πR0 ∂ρ
def def ∂ψ 1 ∂ψ Bp R0
i = µ = = = , (3.9)
∂Φ 2πB0 ρ ∂ρ B0 ρ
def 1 ∂Φ
q = =
µ ∂ψ

En la práctica rara vez se utilizan todos los términos del lado derecho del sistema
de ecuaciones de transporte 3.7, y la matriz de transporte completa 3.8 se utiliza en el
modelado de transporte. Más aún el empleo del conjunto de ecuaciones (3.7 es excesivo.
De hecho el cómputo de la densidad electrónica ne puede reemplazarse por mediciones
experimentales. El código ASTRA proporciona un camino fácil para retener los términos
y ecuaciones de los conjuntos 3.7 y (3.8) que son necesarios para el problema especı́fico
que se considera. Estos son codificados como un conjunto especial de instrucciones.
Por otro lado, el conjunto de ecs.(3.7) es insuficiente para algunos problemas de interés,
como el calentamiento menor de los iones ó la remoción de “ceniza” de hélio. En tales
casos, el código puede complementarse con ecuaciones de difusión adicionales.
Adicional al conjunto principal de ecuaciones de transporte, ecs. (3.7), el código AS-
TRA proporciona un variedad de modulos para la descripción de diversos procesos en el
plasma de un Tokamak y para el cómputo de coeficientes de transporte en las ecs.(3.7)-
(3.8). Los procesos modelados por estos módulos incluyen; calentamiento auxiliar y con-
ducción de corriente, el comportamiento de especies presentes en menor proporción, el
borde del plasma, plasma SOL, estabilidad MHD, etc.

48
3.3 Fuentes y Sumideros

Las fuentes de electrones Se en la primera de las ecuaciones del conjunto de ecs.(3.7), es:
(e) (i)
Se = sion N ne + sion N ni − srec ni ne (3.10)

donde N es la densidad de los átomos neutros del gas de trabajo, srec es la tasa de
(e) (i)
recombinación, y sion y sion son las tasas de ionización por colisión por medio de iones y
electrones. Las fuentes de energı́a por iones Pi y electrones Pe están definidas como:

Pe = POH − PΓ − Pei − PeRAD − PeN + PeH + P F U S (3.11)

Pi = PΓ + Pei + PiN + PiH

estas definiciones incluyen el calentamiento óhmico POH , la perdida de potencia por ra-
diación PeRAD , la potencia de calentamiento auxiliar para iones PiH y electrones PeH , in-
tercambio de energı́a entre iones y electrones

D E Γ ∂ (n T )
e i i
PΓ = (∇ρ)2 (3.12)
ne ∂ρ
3me ne
Pei = (Te − Ti )
mi τe
y términos de pérdidas debidas a procesos atómicos
 
(e) 3
PeN = sion N εion + srad N ε1 + srec ni Te ne (3.13)
2
3
PiN = (sion N TN + scx N (TN − Ti ) − srec ne Ti ) ni
2
donde TN es la temperatura de los átomos neutros, scx y srad son las tasas de intercambio
de carga y radiación, con; εion = 13.6eV , ε1 = 10.2eV .
Aunque la corriente de amarre jBS esta expresada en e las ecs. (3.8) en términos
de las fuerzas termodinámicas, la corriente de amarre junto con la corriente de con-
ducción jCD pueden tomarse como fuentes del flujo poloidal. La densidad de la corriente
de conducción jCD ası́ como N , TN , PeH y PiH deben ser provistas por módulos separados
los cuales están incluidos en el código. Oscilaciones de diente de sierra y el compor-
tamiento del plasma en la región del diversor son ejemplos de procesos que requieren de
un tratamiento por separado y los módulos se encuentran incluidos en el código.
49
Para cerrar el sistema de ecuaciones de transporte (3.7) necesitamos, adicionalmente
a las fuentes y sumideros de las ecuaciones (3.10)-(3.13), la matriz de transporte (3.8),
la conductividad eléctrica σk y la corriente de conducción jCD , también se debe especi-
ficar las funciones de superficie V 0 , I, G1 , G2 , G3 . Después con las condiciones iniciales
apropiadas y de frontera, podemos determinar el tiempo de evolución de las distribu-
ciones radiales en la densidad electrónica ne , las temperaturas Te,i y la densidad de cor-
riente.

3.4 Condiciones Iniciales

Usualmente las condiciones iniciales no juegan un papel primordial en el modelado del


transporte en Tokamaks, porque las condiciones normales de descarga en estado esta-
cionario siguen bajo estudio. Sin embargo, para el análisis de procesos transitorios,
como la corriente de elevación del plasma, las condiciones iniciales tienen un papel
fundamental. Cuando están disponibles, las cantidades evaluadas experimentalmente
pueden utilizarse para proporcionar las distribuciones, pero estas no siempre lo están,
especialmente para el flujo poloidal ψ y la correspondiente densidad de corriente. La
simulación de ciertos procesos puede ser sensible a la elección inicial del perfil de la
densidad de corriente, y esta elección puede requerir de cierto cuidado. Usualmente es
suficiente establecer las condiciones iniciales auto-suficientes, para proporcionar una
suave evolución del plasma durante los instantes iniciales de la corrida.
Para las primeras 3 ecuaciones de (3.7), las condiciones iniciales naturales son:

ne (ρ, t) |t=0 = ne0 (ρ)

Te (ρ, t) |t=0 = Te0 (ρ) (3.14)

Ti (ρ, t) |t=0 = Ti0 (ρ)

Una condición inicial similar para el flux poloidal se leerı́a:

ψ (ρ, t) |t=0 = ψ0 (ρ)

Sin embargo esta no tiene un uso práctico porque usualmente es conveniente expresar
la densidad de corriente jk o la transformación rotacional µ, en vez del flux poloidal ψ tal
50
cual. Por lo tanto, las condiciones iniciales apropiadas son alternativas,

jk (ρ, t) |t=0 = j0 (ρ) ó µ (ρ, t) |t=0 = µ̄ (ρ) (3.15)

Estas condiciones involucran ya se la segunda o la primera derivada de la función incógnita


ψ (ρ), respectivamente. Con cualquier condición inicial, ψ (ρ) debe ser determinada ya sea
empleando

def hj · Bi
jk =
B0

 
2πR0 2 ∂ −1 ∂ψ
= J G2 J
µ0 V 0 ∂ρ ∂ρ
o las ecs.(3.9). Una cantidad mejor definida es la corriente total del plasma Ipl , que siem-
pre se mide de experimentos, mientras la distribución de corriente durante la primera
fase de la descarga en el Tokamak usualmente es desconocida. Por lo tanto es razonable
requerir que Ipl tome un valor dado, y las normalizaciones de la densidad de corriente
o los perfiles de transformación rotacional son ajustados para ser auto-consistentes uti-
lizando

jk (ρ, t) |t=0 =

j0 (ρ) ,
R ρB −2 0 (3.16)
0 J j0 V dρ = 2πR0 Ipl

µ̄ (ρ) µ0 Ipl
µ (ρ, t) |t=0 = (3.17)
µ̄ (ρB ) 2πB0 ρB G2 (ρB )
respectivamente. La condición de frontera J (ρB ) = 1 es utilizada en la ec.(3.16).
También es conveniente y útil la condición inicial con la corriente del plasma dis-
tribuida de acuerdo a la condición de estado estacionario

∂2ψ ∂ψ
=0⇒ = Upl (ρ) = cte
∂t∂ρ ∂t
Utilizando la ley paralela de Ohm

jk = σk Ek + jBS + jCD

51
y suponiendo Ḃ0 (t = 0) = 0, re-escribimos la condición como:

 2πρ ρ
jk − jBS − jCD |t=0 = 0 σk (ρ) Upl = Cσ 0 σk (ρ) (3.18)
V V

Al igual que las ecs.(3.17) y (3.16) el factor Cσ en la ec.(3.18) debe obtenerse de la


condición que plantea la corriente total como Ipl . Aunque el tiempo de relajación de la
corriente en un TOKAMAK es largo, la condición estacionaria inicial ec.(3.18) es ampli-
amente utilizada cuando las mediciones directas de la corriente, no están disponibles.
La condición de la ec.(3.18) tiene relevancia fı́sica cuando procesos no difusivos causan
una rápida relajación de la corriente, o cuando una larga fase de relajación previa no es
de interes.

3.5 Condiciones de frontera para la densidad y la temperatura

Las condiciones de frontera para las ecs.(3.7) en el eje magnético ρ = 0 están impuestas
por la geometrı́a del problema y el requerimiento que todos los flujos se desvanecen en
ρ=0

Γe |ρ=0 = qe |ρ=0 = qi |ρ=0 = 0 (3.19)

Las condiciones de frontera en ρ = ρB para las primeras 3 ecuaciones del sistema de


ecs.(3.7) usualmente toma una de dos formas posibles

ne (ρB ) = neB (t) ó Γe (ρB ) = ΓeB (t)

Te (ρB ) = TeB (t) ó qe (ρB ) = qeB (t) (3.20)

Ti (ρB ) = TiB (t) ó qi (ρB ) = qiB (t) (3.21)

Todas las funciones de la derecha en las ecs.(3.20), además de depender explı́citamente


del tiempo t, pueden también depender de todos los demás parámetros del plasma. Las
condiciones de frontera de las ecs.(3.20) pueden proponerse en diferentes combinaciones.
Vale la pena notar que, en general la cuestón de las condiciones de frontera no es tan
directa. De hecho, en la configuración del espacio,

52
se supone que la frontera del plasma debe coincidir con alguna superficie de flujo SB .
Esta superficie fronteriza sirve como entrada para un solver de la ecuación de equilibrio
de Grad-Shafranov. El solver devuelve una posición de esta superficie en coordenadas
de flujo ρB . Aún en el caso más simple, cuando la frontera del plasma SB está fija en
el espacio y no tiene movimiento, el correspondiente valor ρB puede variar en el tiempo
dependiendo de la presión del plasma y las distribuciones de densidad de corriente. Por
lo tanto, la frontera del plasma ρB (t) se mueve en el marco de referencia Lagrangiano de
coordenadas de flujo ρ. Por lo tanto el problema bajo consideración es siempre el de una
frontera móvil.
También es posible que las condiciones de frontera requeridas no puedan ser expre-
sadas explı́citamente en alguna de las formas dadas en las ecs.(3.20). Este puede ser el
caso cuando condiciones de frontera más elaboradas son impuestas al considerar mod-
elos tipo SOL. La correcta descripción de estos casos pueden requerir el uso de teorı́a
cinética ó MHD. En el código ASTRA, estas situaciones son implementadas mediante el
llamado de subrutinas especiales.

3.6 Condiciones de frontera para el flujo poloidal

La condición de frontera del lado izquierdo para la cuarta ecuación del conjunto de
∂ψ
ecs.(3.7) es directamente, ∂ρ |ρ=0 = 0 . La condición de frontera en ρ = ρB es más
complicada. Generalmente, la condición debe de obtenerse resolviendo un problema de
magnéto-estática en la región exterior (respecto al plasma) tomando en cuenta todas las
corrientes poloidales con sus inductancias mutuas, magnetización del núcluo de hierro,
corrientes inducidas en la cámara de vacio, etc. Es mucho más práctico usar una de las
tres posibilidades descritas a continuación.

3.6.1 Corriente preescrita del plasma

Esta es la condición de frontera utilizada con mayor frecuencia para el modelado de


transporte. Haciendo uso de

53
Z Z
def 1
Ipl = j · dSζ = (j · ∇ζ) d3 x
Sζ 2π V
Z ρ Z ρ 0
1 J V
= V 0 jtor dρ = j dρ
2πR0 0 2πR0 0 J 2 k
G2 ∂ψ 2πB0 2πR0
= = ρG2 µ = G2 Bp
µ0 ∂ρ µ0 µ0

escribimos la condición como

∂ψ µ0
|ρ=ρB = Ipl (t)
∂ρ G2 (ρB )

donde Ipl (t) es la corriente total del plasma con una dependencia del tiempo preescrita.

3.6.2 Búcle de voltaje preescrito

Esta usualmente no es una buena elección como condición de frontera. De hechom en


plasmas ohmı́cos, esta condición de frontera da lugar al surgimiento de una inestabilidad
térmica. Por otro lado, esta inestabilidad es muy lenta y fácilmente puede ser estabi-
lizada. Por lo tanto, esta condición también se encuentra incluida en el código como una
posible opción. En el caso más simple donde la forntera del plasma se encuentra fija se
lee

∂ψ
|ρ=ρB = UplB (t)
∂ρ
con UplB siendo el búcle de voltaje a la frontera preescrito.

3.6.3 Ecuación del circuito externo

Esta proporciona una condición de frontera que incluye las 2 condiciones previas como
casos particulares. Utiliza una descripción simplificada de un circuito externo con la
única ecuación
d
Uext = (Lext Ipl ) + UplB
dt
donde Lext es la inductancia externa del plasma, y UplB = Upl (ρ = ρB ) es el voltaje de búcle
calculado al borde

∂ρB ∂ψ ∂ρB
UplB = Utor |ρB = Upl (ρB ) + 2πB0 ρB µ (ρB ) = |ρB +2πB0 ρB µ (ρB )
∂t ∂t ∂t
54
Uext es una función de tiempo dada asociada con el voltaje producido por la bobina pri-
maria.

3.7 Cierre de las ecuaciones de equilibrio y transporte

Las cuatro ecuaciones del sistema (3.7) no conforman un sistena cerrado, aún cuando
todos los flujos y lados derechos se encuentren determinados. De hecho, las funciones
V (ρ, t), J (ρ, t), G1,2 (ρ, t) incluidas en el sistema (3.7) aún permanecen desconocidas. Más
aún, la variable independiente ρ también es una función desconocida del tiempo y co-
ordenadas espaciales: ρ = ρ (r, z, t). Es la solución de la ecuación de Grad-Shafranov:
 
4∗ ψ = r2 div ∇ψ
r 2 = −4π 2 µ r 2 ∂p + I ∂I , la que da las relaciones instantaneas de todas
0 ∂ψ ∂ψ

las coordenadas de flujo, tales como ρ, V, J,etc. entre ellas y con respecto al sistema
coordenado de laboratorio {r, z}. Estas relaciones pueden variar en el tiempo debido a
la evolución temporal de parámetros del plasma como lo describen las ecs.(3.7). Por lo
tanto para estudiar la evolución del plasma en un Tokamak debemos resolver el conjunto
de ecuaciones de transporte de una dimensión (3.7) y (3.8) de manera simultanea con la
ecuación de equilibrio de dos dimensiones.
!
2πµ R J 2πµ R 2 J 2 r 2 ∂p
0 0 0 0
4∗ ψ =
2 2 jk +

2 2
− 2
B /B0 B0 ρµ B /B0 R0 ∂ρ
" #
J2 r2
 
R0 J ∂p ∂p
= 2πµ0 R0
2 2 jk + − (3.22)
B /B0 B0 ρµ ∂ρ B0 R0 ρµ ∂ρ

La ecuación de equilibrio (3.22) depende de t paramétricamente pero no explı́citamente.


Fı́sicamente esto significa asumir que el plasma esta siempre en equilibrio, y no se con-
sideran el pocesos de relajamiento. La justificación es que de la relajación al equilibrio, es
un proceso varios órdenes de magnitud más rápido que cualquier proceso de transporte

3.8 Parámetros del Tokamak TCABR

El laboratorio de fı́sica de plasmas del Instituto de Fı́sica de Sao Paulo en la universidad de


Sao Paulo es uno de los laboratorios participantes en el “Proyecto de Coordianción de In-
vestigación” en la “Investigación conjunta empleando pequeños tokamaks” de la Agencia
Internacional de Energı́a Atómica. Los experimentos conjuntados (JE) son organizados

55
por los laboratorios participantes en cooperaciï"Čï·şn con la IAEA. El objetivo mayor de
las actividades es el de estimular la colaboración entre grupos de investigación y labora-
torios, la promociï"Čï·şn de experimentos, y la creación de sinergias para contribuir de
manera efectiva al desarrollo en la ciencia de la fusiï"Čï·şn nuclear.
Los parámetros del Tokamak TCBAR son [7]:

Radio mayor 0.61 m


Radio menor 0.18 m
Campo magnético toroidal BT = 1.07 T
Corriente máxima al plasma Ipmax = 100 kA
Densidad electrónica máxima promedio ne = 4.5 ∗ 1019 part
m3
Temperatura electrónica máxima Te = 500 eV
Temperatura iónica máxima Te = 200 eV
Duración del pulso 100 ms

Tabla 3.1: Parámetros del Tokamak TCABR

El programa de investigación cientı́fica que se desarrolla en el TCABR se compone de


las siguientes lı́neas de investigación:

• Interacción de ondas electromagnéticas RF en el rango Alfvén con plasmas de Toka-


mak.

• Fı́sica del borde de plasma en regı́menes óhmicos y de confinamiento mejorado del


borde del plasma .

• Rotación del plasma toroidal y poloidal.

• Desarrollo de diagnósticos del plasma.

• Control remoto y adquisición de datos.

3.8.1 Análisis del transporte del Tokamak TCABR con el código ASTRA

Conociendo los perfiles de la densidad del plasma, temperatura electrónica e iónica, se


puede establecer el modelo para analizar los diferentes coeficientes de transporte y los
parámetros mas importantes para el análisis del transporte en el Tokamak TCABR.

56
Datos experimentales y modelo de ASTRA

Los datos experimentales son proporcionados en el formatao para ASTRA en la forma


TCABR - Brazilian tokamak
Name Time Value Error
AB .18
RTOR 0.615
BTOR 1.1
ELONM 1.
ELONG 1.0
TRIAN 0.0
IPL 0.12
AMJ 2.
ZMJ 1.
NNCL 1. .001 cold 1019 m−3
NNWM .0001 warm 1019 m−3
QICR 0.1
Number of points in arrays (radial positions are assumed to be distributed
equidistantly in the radial coordinate - equivalent radius of flux surfase)
POINTS 16
Name Time Values Aerr Rerr
NE 0.20 3.43 3.40 3.32 3.20 3.04 2.87 2.70 2.52 2.36 2.20 2.04 1.87 1.67 1.42 1.13 0.80
TE 0.09 0.60 0.57 0.56 0.54 0.50 0.43 0.34 0.25 0.16 0.08 0.06 0.04 0.03 0.02 0.01 0.008
TI 1.60 .40 0.38 .35 .32 .28 .25 .21 .17 .14 .11 .08 .05 .02 .009 .008 .006
ZEF 1.60 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80 1.80
END
El modelo que se introduce en ASTRA se escribe en la forma
NEQUIL=20
! MESHEQ=21;
NEUT:.02:::N;
! Normalized to 1 MW Gaussian function: CAR1=P 0 ∗ exp{−[(r − cf 3)/(a ∗ cf 4)]2 }

57
FGAUSS(CF3*0.5,CF4*0.03,CAR1):;

NE:EQ; NE=CF1*3.0*FPR; NEB=.3*CF1;


PEX=QICR*FRAMP(1.,2.)*CAR1;
! DN=0.3*CHE1*HAALC+HNGSE;
DN=0.3*CHE1/SQRT(1+CF11*PEX*CF12*SQRT(AMJ*TI)/(ZMJ*BTOR));
SN=SNEBM;
SNN=SNNEU;
CF5=0.3;
TE*:EQ; TE=CF2*0.5*FPR; TEB=CF2*0.05;
PRAD=CF6*(.02+.03*FPR);
PE=PJOUL-PRAD+PEX;
! HE=CHE3*HAALC+HNGSE;
HE=15.*CHE3/SQRT(1+CF11*PEX*CF12*SQRT(AMJ*TI)/(ZMJ*BTOR));

TI*:EQ; TI=CF5*FPR; TIB=CF5*0.1;


! PIX=QICR*FRAMP(1.,2.)*CAR1
! PI=PINEU+PIX; XI=1.8*CHE4;

CU:EQ; CC=CNHR; CU=CC;


! HC=HCSA; DC=DCSA; XC=XCSA;

Radial output
Te\TE\−2; j\CU; jbs\CUBS; tria\TRIA;
Ti\TI; ne\NE
NEX; mu\MU\1; shif\SHIF;
Tex\TEX\-2; Upl\UPL; He\HE\-1; beta\BETAJ;
Tix\TIX; Qe\QETOT; q\1./MU\5; elon\ELON;

PiN\-PINEU\-3; sigm\CC; He\HE\-1; Pei\PEICL\-3;


PeN\PENEU+PENLI\-3; N\NN\1; Hinc\CHE4*HNGSI\-1;PEGN\PEGN;

58
0.12 4.5

0.1

3.5

0.08 3

2.5

0.06

0.04 1.5

0.02

0.5

0 0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

(a) (b)

Figura 3.1: (a) Perfil de β = nTe /B 2 (2µ0 ) para dos tiempos diferentes. (b) Perfil de la
densidad de corriente del plasma para dos tiempos diferentes. En este caso vemos que
ambas cantidades se saturan rápidamente con el tiempo.

PGN\PEGN+PIGN\-3; TN\TN; Hi\XI\-1; PeiG\PEIGN\-3;


PEC\PEX\-3; Zeff\ZEF; HAlc\CHE1*HAALC\-1;PIGN\PIGN;

Time output

Pinp QEXB+QJOULB; hnei NEAVB; ne NECHC; Te0 TEC -3;


V(a) UPLB; ; Ti0 TIC; hTei TEAVB -3;

Pech QEXB; taux WTOTB/(QJOULB) .1; Wblk WTOTB -1; taui TAUEIB;
We WEB -1; tauE TAUEB .1; Wi WIB -1; IT89 TITER .1;

Poh QJOULB; li LINTB; betj BETAJB 1; q(0) 1./MUC;


Prad QRADB; Ibs IBSB; Xs CMHD4;
Poh\PJOUL\-3;Pec\PEX\-3;

59
2 3

1.8

2.5
1.6

1.4

1.2

1
1.5

0.8

0.6 1

0.4

0.5
0.2

0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18
0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

(a) (b)

Figura 3.2: (a) Perfil de la densidad del plasma para dos tiempos diferentes. (b) Perfil
del factor de seguridad para dos tiempos diferentes. En este caso vemos que ambas
cantidades se saturan rápidamente con el tiempo.

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

Figura 3.3: (a) Perfil del flujo de calor de los electrones para dos tiempos diferentes. En
este caso vemos que ambas cantidades se saturan rápidamente con el tiempo.

60
0.5 0.16

0.45
0.14

0.4

0.12

0.35

0.1
0.3

0.25 0.08

0.2
0.06

0.15

0.04

0.1

0.02
0.05

0 0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

(a) (b)

Figura 3.4: (a) Perfil de la temperatura electrónica para dos tiempos diferentes. (b) Perfil
de la temperatura iónica para dos tiempos diferentes. En este caso vemos que ambas
cantidades se saturan rápidamente con el tiempo.

Parámetros obtenidos de ASTRA

Como ejemplo se muestra el comportamiento de los diferentes parámetros obtenidos a


través del código ASTRA.
Las gráficas elaboradas con los datos calculados mediante ASTRA para el TCABR
muestran que los procesos de estabilización del plasma ocurren en cortos periodos de
tiempo, pues para ambas simulaciones los perfiles obtenidos son prácticamente iguales.
Los valores desarrollados para el factor β cerca del borde están dentro de los valores más
elevados desarrollados por los Tokamaks, sin embargo los parámetros operacionales y
los valores de plasma desarrollados se ubican muy lejos de una tasa de fusión relevante.
La densidad del plasma disminuye su valor del centro al borde hasta en un orden de
magnitud, la distribución de partı́culas a lo largo del radio disminuye demostrando una
mayor concentración del plasma al centro del toroide y una difusión radial de materia,
el cambio de las dimensiones en la geometrı́a anidada con repecto al radio contribuye a
la disminución del valor de densidad, sin embargo este descenso no ocurre de forma lin-
eal, comportamiento que hace evidente la acción de mecanismos adicionales de difusión
como colisiones energéticas entre las partı́culas, y el efecto de la trampa mganética para
confinar la materia.
La temperatura de los electrones disminuye de forma radial hasta prácticamente un

61
orden de magnitud. La temperatura de los iones presenta un comportamiento similar,
manteniendose por debajo de la temperatura de los electrones e indicando transferencia
de energı́a de estas partı́culas hacı́a los iones. El plasma pierde energı́a radialmente de
forma no lineal. El flujo de calor debido a los electrones incrementa con el radio, ya que
refleja el comportamiento de perdida de energı́a en el sistema. La densidad de corriente
disminuye radialmente, con un comportamiento casi lineal hasta cero en el borde del
plasma.

Conclusiones

La teorı́a neoclásica toma en cuenta el efecto de partı́culas atrapadas, producto de la ge-


ometrı́a del campo magnético. En un plasma con baja tasa de colisiones, el confinamiento
en el régimen de la trayectoria libre promedio es dominado por la tasa de pérdida que
poseen las partı́culas atrapadas en los espejos magnéticos. Los coeficientes de transporte
neoclásico representan el transporte en una geometrı́a de toroide, debida únicamente a
las colisiones de Coulomb entre los centros guı́a.
Las colisiones entre centrso guı́a las desvı́an de sus trayectorias inı́ciales, en una
especie de caminata aleatoria, originando lo que se conoce como difusión radial. Las
partı́culas “brincan” de una superficie magnética a otra debido a estas colisiones. Los
datos predichos mediante la teorı́a neoclásica tienen una concordancia, aunque irregular,
respecto de datos experimentales con plasmas que no son portadores de corrientes in-
tensas, y solo para el transporte de calor en iones. Sin embargo los resultados obtenidos
para la conducción de calor en electrones y el transporte de partı́culas están subestima-
dos. Todos estos resultados se obtienen en máquinas Tokamak. Los modelos neoclásicos
predicen el transporte con cierto éxito bajo circunstancias especı́ficas.
Cuando la deducción de las ecuaciones de transporte se efectua en una geometrı́a
toroidal, no es posible emplear la simetrı́a poloidal para establecer los gradientes poloidales
de temperatura y densidad. Los términos de fricción se generalizan para incluir coefi-
cientes numéricos ajustados y el término de fricción térmica, que surgen de la depen-
dencia a la frecuencia de colisión que tiene la temperatura. Estas consideraciones dan
origen a una expresión para el flujo radial de partı́culas sobre el ángulo poloidal.
Los efectos de la gometrı́a toroidal, que están asociados con el gradiente de campo

62
magnético y rezagos en la curvatura, dan como resultado un incremento de al menos un
orden de magnitud en los flujos de transporte radial, sobre los valores clásicos.
La naturaleza helicoidal del campo, acarrea efectos adicionales en el transporte. La
formación de un pozo magnético en el campo en la parte exterior del toro y un máximo
en el campo en el borde interior del toro. Las partı́culas atrapadas en el pozo magnético,
producen un flujo de partı́culas radial cuando tienen el tiempo suficiente para ejecutar
una órbita de partı́culas atrapadas, antes de ser desviadas 900 .
Los valores de operación del reactor TCABR propician que el modelo neoclásico sea
apropiado para la descripción del sistema. En el dispositivo no se observan fenómenos
de turbulencia, por lo que los datos predichos por el modelo mantienen relación con las
mediciones experimentales. La núla variación entre los perfiles de las simulaciones con
distintos tiempos de ejecucı́ón muestran la tendencia al estado estacionario del modelo,
pues no se manifiestan perturbaciones en los perfiles debido a la evolución temporal del
campo electro-magnético.

63
Bibliografı́a

[1] Jardin, Stephen. Computational Methods in plasma Physics, Florida: Chapman &
Hall/CRC. 2010. [Cap. 1, 5, 6]

[2] Stacey, Weston M. Fusion: An Introduction to the physics and technology of magnetic
confinement fusion, Atlanta,GA: Wiley-VCH. 2010. [Cap. 3]

[3] Akhiezer A.I., Akhiezer I.A., Plovin R.V., Sitenko A.G., and Stepanov K.N., Plasma
Electrodynamics, Oxford: Pergamon Press, 1975, [Cap. 1].

[4] Pereverzev, G.V. Yushmanov, P.N. ASTRA: Automated System for TRansport Analysis,
Munchen: Max-Planck-Institut fur Plasmaphysik. 2002. [Cap. 2, 3, 4]

[5] D’haeseleer, W.D. Hitchon, W.N.G. Callen, J.D. Shohet, J.L. Flux Coordinates and
Magnetic Field Structure: A Guide to a Fundamental Tool of Plasma Theory, Berlin:
Springer-Verlag. 1990. [Cap. 3, 6]

[6] Bird, Robert. Fenómenos de Transporte, 2da ed. México: Limusa Wiley. 2010. [Cap.
3. Apend. D]

[7] TCABRJE - Joint Experiment on TCABR - Host Laboratory Experiment.


http://tcabrje.if.usp.br

[8] Callen, James D. Fundamentals of Plasma Physics, Wisconsin: University of


Wisconsin-Madison.2006. [Cap. 3, 5]

64

También podría gustarte