Physics">
Deducción de Ecuaciones de Transporte Noeclásico en Un Tokamak Axisimétrico
Deducción de Ecuaciones de Transporte Noeclásico en Un Tokamak Axisimétrico
Deducción de Ecuaciones de Transporte Noeclásico en Un Tokamak Axisimétrico
Diciembre, 2014.
Contenido
0.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Aproximación de fluido 28
2.1 Teorı́a general de momentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1
Resumen
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 θ
3
Capı́tulo 1
• Descripción de fluido
• Descripción cinética
• Descripción Giro-cinética
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
ρ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.
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)]
(x, y, z, vx , vy , vz )
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
−∂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
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
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
x → xi
v → vi
∆V ≡ ∆x∆y∆z
10
Y el volumen del espacio de velocidades:
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
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
δ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
B = hBm i + δBm
ρm
m m
q = ρq + δρq
J = hJm i + δJm
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)
∂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
E = −∇φ
∂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
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
∂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 = f (xg , εg , µ, t)
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.
• 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.
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 τ
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
ωτ 1
20
L1
L1
ωτ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.
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
1
fe = ρe E + [j × B] ,
c
1 ∂B
rotE = −
c ∂t
∇·B=0
22
1 ∂E 4π
rotH = − + j
c ∂t c
∇ · E = 4πρe
H = B − 4πM
ds
ρm θ =q
dt
X ∂uα j2
q≡ παβ + (x∇θ) +
∂xβ σ
α,β
1 ∂E 4π
σ |E|
c ∂t
c
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
4π
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
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
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
27
Capı́tulo 2
Aproximación de fluido
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.
28
• Temperatura (J, eV ):
mvr2 mvT2
Z
1
T ≡ d3 v f=
n 3 2
mvr2
Z
q≡ d3 vvr ( )f.
2
• Presión (N/m2 ):
mvr2
Z
p≡ d3 v f = nT.
3
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
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 ).
(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
σk = σSp
σ⊥ = σ0
Jk J⊥
Re = −ne e(b̂ + )
σk σ⊥
y
Jk2 2
J⊥
Qe = + − Qi
σk σ⊥
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 + π.
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
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.
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)
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
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
λe ∼ λi ,
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
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
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.
• 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
ζ = −ϕ
45
donde g es el determinante del tensor métrico
D (x, y, z)
g= = (∇a∇θ∇ζ)−2
D (a, θ, ζ)
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
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πρ
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
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:
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.
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,
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
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
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.
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.
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
∂ψ µ0
|ρ=ρB = Ipl (t)
∂ρ G2 (ρB )
donde Ipl (t) es la corriente total del plasma con una dependencia del tiempo preescrita.
∂ψ
|ρ=ρB = UplB (t)
∂ρ
con UplB siendo el búcle de voltaje a la frontera preescrito.
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.
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 ρµ ∂ρ
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]:
3.8.1 Análisis del transporte del Tokamak TCABR con el código ASTRA
56
Datos experimentales y modelo de ASTRA
57
FGAUSS(CF3*0.5,CF4*0.03,CAR1):;
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;
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.
Time output
Pech QEXB; taux WTOTB/(QJOULB) .1; Wblk WTOTB -1; taui TAUEIB;
We WEB -1; tauE TAUEB .1; Wi WIB -1; IT89 TITER .1;
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.
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
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]
64