Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
SIMULACIÓN DE LA INYECCIÓN DE UN HAZ NEUTRO EN TOKAMAKS CON
TRIANGULARIDAD POSITIVA Y NEGATIVA
SIMULATION OF NEUTRAL BEAM INJECTION IN POSITIVE AND NEGATIVE
TRIANGULARITY TOKAMAKS
F. Sheffield*1 y R. Farengo1,2
1Instituto Balseiro – Universidad Nacional de Cuyo
San Carlos de Bariloche – Río Negro – Argentina
2Centro Atómico Bariloche - Comisión Nacional de Energía Atómica
San Carlos de Bariloche – Río Negro – Argentina
Recibido: 09/02/2023; Aceptado: 14/09/2023
En este trabajo se estudia la dinámica de partículas cargadas de alta energía en tokamaks con triangularidad positiva
y negativa mediante la simulación de la inyección de un haz neutro en equilibrios MHD analíticos usando el código
de órbita exacta FOCUS. Se observa que la configuración de triangularidad negativa presenta menores pérdidas de
partículas no termalizadas y una deposición de energía sobre el plasma más localizada que su contraparte.
Palabras Clave: plasma, tokamak, triangularidad, FOCUS.
In this paper the dynamics of high energy charged particles in negative and positive triangularity tokamaks is studied
through the simulation of a neutral beam inyection on analytical MHD equilibria using the exact orbit code FOCUS. It
is observed that the negative triangularity configuration presents fewer losses of non-thermal particles, and an energy
deposition more localized than its counterpart.
Keywords: plasma, tokamak, triangularity, FOCUS.
https://doi.org/10.31527/analesafa.2023.34.4.97 ISSN 1850-1168 (online)
* facundo.sheffield@ib.edu.ar
©2023 Anales AFA
I. INTRODUCCIÓN
Los tokamaks son dispositivos de confinamiento magnético de plasma con geometría toroidal utilizados en investi-
gaciones sobre fusión nuclear controlada. Los tokamaks más avanzados en operación tienen triangularidad positiva, es
decir, poseen una sección transversal en forma de ’D’. Esto se debe a que se ha observado que las configuraciones de
triangularidad positiva resultan ser más estables frente a perturbaciones magnetohidrodinámicas macroscópicas que una
configuración de sección transversal circular [1].
Sin embargo, en los últimos años se han realizado experimentos con equilibrios con triangularidad negativa (sección
transversal en forma de ’D’ invertida) en los tokamaks TCV [2] y DIII-D [3] que indican que se pueden obtener parámetros
asociados al confinamiento y a la estabilidad similares a los obtenidos en un equilibrio con triangularidad positiva. Más
aún, se sabe que una geometría con triangularidad negativa tendría una serie de ventajas ingenieriles que reducirían
problemas de sobrecalentamiento en el tokamak [4,5].
Es por ello que resulta de interés estudiar en más detalle las configuraciones con triangularidad negativa. En este
trabajo se decidió estudiar la inyección de un haz neutro, que es un método comúnmente usado para calentar el plasma, en
configuraciones de triangularidad positiva y negativa. Para ello se simuló la dinámica de partículas inyectadas mediante
un haz neutro de deuterio en equilibrios analíticos con triangularidad positiva y negativa. En la sección II se discute
cómo se obtuvieron los equilibrios analíticos con triangularidad positiva y negativa, en la sección III se muestran los
resultados obtenidos en las simulaciones de la inyección del haz sobre ambas geometrías, y en la sección IV se presentan
las conclusiones.
II. EQUILIBRIOS MHD
Los equilibrios analíticos utilizados fueron obtenidos en el marco de la teoría magneto-hidrodinámica (MHD) ideal
siguiendo un procedimiento análogo al desarrollado por R. Farengo [6]. Para calcular equilibrios MHD toroidales con
rotación se puede utilizar la siguiente ecuación [7]:
∂2ψ
∂z2+∂2ψ
∂r2−1
r
∂ ψ
∂r=−µ0r2exp(mΩ2r2
2T)×
dG
dψ+Gmr2d
dψΩ2
2T−FdF
dψ
(1)
donde se asume que la temperatura total del plasma Tes constante sobre las superficies magnéticas del equilibrio. En la Ec.
(1), ryzson las componentes radial y vertical en coordenadas cilíndricas, mes la suma de las masas iónica y electrónica
del plasma, Ωla velocidad de rotación toroidal del plasma, ψel flujo poloidal dividido por 2π,Ges función de la presión
del plasma y Fes función de la corriente poloidal. Las cantidades ψyFse relacionan con el campo magnético de la
siguiente manera:
B=F
r
ˆ
θ+1
r∇ψ׈
θ,(2)
mientras que G=pexp(−mΩ2r2
2T), siendo pla presión del plasma. Asumiendo que la cantidad Ω2/(2T)es constante, es
necesario especificar la dependencia de GyFcon el flujo poloidal para poder resolver la Ec. (1). Para ello, se considera
el Ansatz de Solov’ev [8]:
dG
dψ=−S1,FdF
dψ=S2.(3)
De manera que la ecuación a resolver sea lineal. Bajo estas suposiciones se puede adimensionalizar la Ec. (1). En este
trabajo se adimensionalizan las longitudes con respecto a R0, el radio mayor del toroide; el campo magnético respecto a
B0, el valor del campo toroidal de vacío producido por las bobinas en R0; el flujo con B0R2
0; G con B2
0µ0y F con B0R0.
Hecho esto, se introduce el número de Mach (M), definido como la relación entre la velocidad de rotación toroidal en R0
y la velocidad térmica del plasma: M2=mΩ2R2
0/(2T). Realizado este procedimiento se llega a una ecuación de la forma:
∂2ψ
∂z2+∂2ψ
∂r2−1
r
∂ ψ
∂r=r2exp(M2r2)S1−S2,(4)
donde todas las cantidades son adimensionales. Al ser una ecuación lineal su solución puede descomponerse en una
particular y una homogénea. Para la solución particular se utiliza la obtenida por Maschke y Perrin [7]:
ψp(r,z) = −S1
4M4h1+M2r2−eM2r2i−S2
2z2.(5)
Para la ecuación homogénea, Zheng et al. [9] mostraron que la solución puede escribirse como
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
ψh(r,z) = ∑
n=0,2,...
n/2
∑
k=0
H(n,k,r)zn−2k,(6)
donde la forma explícita de H(n,k,r)se encuentra en el mismo trabajo escrita como G(n,k,r). La cantidad de términos
de la solución homogénea depende de la cantidad de condiciones de contorno que se quiera imponer, y estos términos
pueden variar según se quiera un equilibrio simétrico en zo no. Para este trabajo, se toman términos simétricos en zhasta
n=6, por lo que la solución homogénea resulta:
ψh(r,z) =
8
∑
i=1
ciψi(r,z),(7)
con ψilas primeras soluciones de la Ec. (6) con potencias pares en z.
Para finalizar el cálculo de los equilibrios es necesario determinar los valores de los coeficientes c1,...,c8, y de S1yS2.
Esto se logra introduciendo ocho condiciones de borde y dos condiciones sobre cantidades globales. Estas condiciones
son
ψ(Rin,0) = 0,(8)
ψ(Rout ,0) = 0,(9)
ψ(Rx,Zx) = 0,(10)
ψ(Ra,Za) = 0,(11)
∂ ψ
∂r(Rx,Zx) = 0,(12)
∂ ψ
∂z(Rx,Zx) = 0,(13)
∂2ψ
∂z2(Rin,0)−αin
∂ ψ
∂r(Rin,0) = 0,(14)
∂2ψ
∂z2(Rout ,0)−αout
∂ ψ
∂r(Rout ,0) = 0,(15)
ZS(ψ)−S1reM2r2+S2
rdrdz =It,(16)
−2S1
VZS(ψ)
ψeM2r2
B22πrdrdz =β.(17)
Las Ecs. (8), (9), (10), (11) corresponden a pedir que los puntos correspondientes sean parte de la última superficie
cerrada (separatriz), donde Rin yRout son los radios interior y exterior del toroide, respectivamente. Los valores de (Rx,Zx)
y(Ra,Za)se obtienen fijando valores geométricos de los equilibrios como la elongación (κ) y la triangularidad (δ). Su
relación viene dada por las ecuaciones:
R0=Rin +Rout
2,a=Rout −Rin
2,
δ=R0−Rx
a,κ=Zx
a,
(18)
con ael radio menor del toroide.
Las Ecs. (12)y(13) implican que el campo magnético sea nulo en el punto (Rx,Zx), conocido como punto X.
En las Ecs. (14) y (15)αin yαout son los valores de curvatura de la separatriz en Rin yRout , respectivamente, y permiten
fijar la curvatura de la sección transversal en esos puntos [10].
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
FIG. 1: Sección transversal del plasma en un tokamak con perfil en forma de D. Se observan las longitudes características de la
sección.
(a) Triangularidad positiva (δ>0) (b) Triangularidad negativa (δ<0)
FIG. 2: Equilibrios de triangularidad positiva y negativa utilizados. Se observan las superficies magnéticas en cada caso.
La visualización de estos parámetros geométricos y su definición resultan más claras al observar la Fig. 1, donde se
muestran las distancias características de la sección transversal del equilibrio.
Las Ecs. (16)y(17) son las condiciones sobre las cantidades globales. La Ec. (16) fija el valor de la corriente toroidal
del plasma, mientras que la Ec. (17) fija el parámetro βdel plasma, que relaciona la presión del plasma con la presión
ejercida por el campo magnético. En ambas ecuaciones se integra dentro de la sección transversal, donde ψ>0.
Las ocho condiciones de borde permiten dejar expresados los coeficientes c1,..., c8en función de S1yS2de forma ana-
lítica. Una vez hecho esto, se pueden obtener los valores de S1yS2resolviendo las condiciones globales numéricamente.
En el caso de este trabajo, se implementó el método de Newton-Rhapson de dos variables para resolverlas.
Vale la pena mencionar que debido al uso del modelo MHD ideal la temperatura y densidad del plasma no quedan
unívocamente determinadas, sino que es necesario considerar un perfil de densidad o temperatura, y luego obtener el otro
mediante la ecuación de estado del plasma, que se supone como la de un gas ideal.
Para este trabajo se emplearon equilibrios estáticos (M=0) de triangularidad positiva y negativa con valores iguales
de radio mayor, menor, B0y corriente a los utilizados por Austin et al. en el tokamak DIII-D [3]. En la Fig. 2se muestran
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
FIG. 3: Perfiles de temperatura ecuatoriales para ambos equilibrios. Se observa que la temperatura máxima es de ∼3keV para δ>0
y de ∼2.6keV para δ<0.
ambos equilibrios estáticos utilizados, que poseen una triangularidad de δ=±0.61, una elongación de κ=1.43 y un
β=0.015, donde este último se define como el promedio de la presión del plasma y la presión magnética sobre todo el
volumen.
Los equilibrios fueron elegidos de manera que su sección transversal sea la misma y por ende que la densidad de
corriente media sea igual. En ambos casos mostrados se utilizó un perfil de densidad uniforme de 0.5×1020 m−3, de
manera que el perfil de temperatura resulta diferente entre los equilibrios, como se muestra en la Fig. 3.
III. INYECCIÓN DE UN HAZ NEUTRO
Para la simulación de la inyección del haz neutro se utilizó el código FOCUS [11], que permite calcular las órbitas
exactas de partículas en un tokamak en forma paralela. El código en sí tiene incorporados módulos que permiten tener en
cuenta procesos de colisiones elásticas e inelásticas, por lo que es posible simular la ionización de un átomo neutro de alta
energía y su deposición sobre el plasma.
En la Fig. 4se muestran las posiciones de ionización de las partículas inyectadas por el haz, junto con una compara-
ción con el código NUBEAM [12]. Dado que este código es una rutina extensivamente testeada, el buen acuerdo en la
ionización sirve para verificar el correcto funcionamiento de FOCUS.
FIG. 4: Posiciones de ionización de partículas simulando la inyección de un haz neutro en el plano X-Y. Se detallan en rojo la simulación
de NUBEAM, en azul los resultados de FOCUS, y en línea punteada el eje magnético de la configuración.
En este trabajo se simuló la inyección de 2×105partículas de deuterio con una energía de 80 keV durante un tiempo de
81 ms sobre los equilibrios mostrados en la Fig. 2. El tiempo fue elegido de manera tal que una proporción significativa
de las partículas se haya termalizado.
Además de estas dos simulaciones se realizó otra sobre un equilibrio de δ<0 análogo pero con un mayor parámetro
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
β=0.017 para que los perfiles de densidad y temperatura tengan igual magnitud. El objetivo de esta simulación fue
corroborar que las diferencias obtenidas se deban al cambio de geometría y no a una diferencia entre los perfiles como la
observada en la Fig. 3.
Partículas escapadas
En la Fig. 5se observa la proporción de partículas que se escapan de los equilibrios a lo largo del tiempo en las
simulaciones.
(a)
(b)
FIG. 5: Partículas perdidas a lo largo del tiempo en configuraciones de triangularidad positiva y negativa. a) Simulación con igual β
en ambos equilibrios. b) Simulación con distinto βy perfiles semejantes de n y T .
Puede verse como durante gran parte de la simulación la fracción de partículas escapadas resulta menor para los casos
de δ<0, igualándose con su contraparte solo para tiempos cercanos a 80 ms, cuando las partículas ya se encuentran
termalizadas.
Dado que el comportamiento es el mismo en ambas simulaciones con δ<0, resulta claro que este fenómeno depende
de la geometría del equilibrio usado. Más aún, el hecho de tener menores pérdidas a tiempos cortos para δ<0 implica
una mayor deposición de energía sobre este equilibrio, ya que las partículas tendrán energías menores antes de escapar de
la configuración. Esto último sería un beneficio adicional para los equilibrios de triangularidad negativa por incrementar
la eficiencia del calentamiento por inyección de haces neutros.
Deposición de energía
Más allá del análisis sobre las pérdidas también se observaron diferencias en la deposición de energía. En la Fig. 6
se observa la región donde se deposita la energía en ambos equilibrios en lo que sería el estado estacionario asociado a
la inyección del haz. En este contexto se obtuvo un estado estacionario de inyección sumando las contribuciones de la
dinámica de las partículas a diferentes tiempos.
Resulta claro que el calentamiento se da en una región más elongada para δ<0. Comparando la Fig. 6con la Fig. 2
parece razonable suponer que esta diferencia se debe a la geometría de las superficies magnéticas, siendo estas más
elongadas para triangularidad negativa. Por otro lado, se observa que la deposición de energía resulta más localizada para
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
(a) δ>0
(b) δ<0
FIG. 6: Potencia por unidad de área depositada por el haz neutro en ambas geometrías, asumiendo una potencia incidente típica de
2.5MW.
la configuración con δ<0, siendo la potencia máxima depositada por unidad de área un 10% mayor. Esto último muestra
una clara ventaja de una configuración con δ<0, ya que una deposición más localizada implica un calentamiento más
localizado y un mayor control sobre el perfil de deposición de energía. Es importante recalcar que en este caso no hubieron
diferencias apreciables entre las simulaciones con δ<0, de manera que se puede afirmar que lo obtenido es un producto
directo del cambio en la geometría de los equilibrios.
IV. CONCLUSIONES
Se estudió la dinámica de partículas asociada a la inyección de un haz neutro en equilibrios con triangularidad positiva
y negativa. Para ello se utilizó un método para obtener equilibrios MHD analíticos sobre los cuales se utilizó el código
FOCUS para simular la inyección por haz neutro de un ensemble de partículas. Se observó que se pierden menos partículas
no-termalizadas en triangularidad negativa, lo que implica una mayor eficiencia en la energía depositada por el haz. Por
otro lado, se observó que la deposición de energía sobre el equilibrio se da en una región más elongada para triangularidad
negativa, y que dicha deposición resulta más localizada en este caso, permitiendo un mayor control sobre la deposición
de energía en un equilibrio de triangularidad negativa.
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
Más estudios podrían plantearse con los métodos usados en este trabajo. En particular, podría compararse el efecto de
inyectar múltiples haces o utilizar equilibrios rotantes (M>0) para observar la dinámica. Aunque en este último caso
debería también incorporarse un campo eléctrico radial que justifique la rotación del plasma.
REFERENCIAS
[1] J. P. Freidberg. Ideal MHD (Cambridge University Press, Nueva York, 2014).
[2] Y. Camenen, A. Pochelon, R. Behn, A. Bottino, A. Bortolon, S. Coda, A. Karpushov, O. Sauter, G. Zhuang y the TCV team.
Impact of plasma triangularity and collisionality on electron heat transport in TCV L-mode plasmas. Nucl. Flusion 47, 510
(2007).
[3] M. E. Austin, A. Marinoni, M. L. Walker, M. W. Brookman, J. S. deGrassie, A. W. Hyatt, G. R. McKee, C. C. Petty, T. L. Rhodes,
S. P. Smith, C. Sung, K. E. Thome y A. D. Turnbull. Achievement of Reactor-Relevant Performance in Negative Triangularity
Shape in the DIII-D Tokamak. Phys. Rev. Lett. 122, 115001 (2019).
[4] M. Kikuchi, T. Takizuka, S. Medvedev, T. Ando, D. Chen, J. Li, M. Austin, O. Sauter, L. Villard, A. Merle, M. Fontana, Y.
Kishimoto y K. Imadera. L-mode-edge negative triangularity tokamak reactor. Nucl. Flusion 59, 056017 (2019).
[5] M. Kikuchi, A. Fasoli, T. Takizuka, P. Diamond, S. Medvedev, X. Duan, H. Zushi, M. Furukawa, Y. Kishimoto, Y. Wu, O.
Sauter, L. Villard, S. Brunner, G. Merlo, G. Zheng, K. Mishra, M. Honda, H. Urano, M. Pueschel, D. Told, A. Fujisawa, K.
Nagasaki y F. Sano. Negative Triangularity Tokamak as Fusion Energy System en Proceedings of 1st International e-Conference
on Energies (MDPI, 2014).
[6] R. Farengo. Extended Solov’ev type equilibria for rotating plasmas with positive and negative triangularity. Phys. Plasmas 27,
122502 (2020).
[7] E. K. Maschke y H. Perrin. Exact solutions of the stationary MHD equations for a rotating toroidal plasma. Plasma Phys. 22,
579-594 (1980).
[8] L. S. Solov’ev. The theory of hydromagnetic stability of toroidal plasma configurations. Sov. Phys. JETP 26, 400-407 (1967).
[9] S. B. Zheng, A. J. Wootton y E. R. Solano. Analytical tokamak equilibrium for shaped plasmas. Phys. Plasmas 3, 1176-1178
(1996).
[10] A. J. Cerfon y J. P. Freidberg. “One size fits all” analytic solutions to the Grad–Shafranov equation. Phys. Plasmas 17, 032502
(2010).
[11] C. F. Clauser, R. Farengo y H. E. Ferrari. FOCUS: A full-orbit CUDA solver for particle simulations in magnetized plasmas.
Comput. Phys. Commun. 234, 126-136 (2019).ISSN: 0010-4655.
[12] A. Pankin, D. McCune, R. Andre, G. Bateman y A. Kritz. The tokamak Monte Carlo fast ion module NUBEAM in the National
Transport Code Collaboration library. Comput. Phys. Commun. 159, 157-184 (2004).ISSN: 0010-4655.
Sheffield et al. / Anales AFA Vol. 34 Nro. 4 (Diciembre 2023 - Marzo 2024) 97-102
