24P038-P
Analysis of classical impurity transport in fusion plasmas by Monte-carlo
Binary Collision Model
モンテカルロ 2 体衝突モデルによる核融合プラズマ中での不純物輸送の検討
Yuu Sawada , Mitsunori Toma , Akiyoshi Hatayama
澤田悠 , 藤間光徳 , 畑山明聖
Faculty of Science and Technology, Keio University
3-14-1, Hiyoshi, Kohoku-ku, Yokohama, Kanagawa 223-8522, Japan慶應義塾大学理工学部〒223-8522 神奈川県横浜市港北区日吉3-14-1
A numerical model for classical/neo-classical cross field diffusion of impurity ions in the magnetic
fusion devices are being developed based on the Binary Collision Monte-Carlo Model. As a first step,
the model validation has been done for a simple case with a uniform magnetic field. Reasonableagreement between the numerical result and the theoretical has been obtained. In addition to theconventional self-diffusion, the cross field transport towards the direction of the background density
gradient predicted by the theory is also observed in the simulation.

1. Introduction
Understanding and control of impurity transport isone of the important issues to reduce the impurity
radiation in the fusion core plasma region. Especially,
impurity accumulation due to the neo-classicaltransport in the magnetic confinement devices, suchas tokamaks, i.e., cross field transport by Coulomb
collisions of impurity ions with background fuel ions
(hydrogen) is one of possible causes. Impurity core
penetration of tungsten impurity ions has beenreported in the experiment during ICRF heating [1].
The reason, however, has not been clearly
understood yet up to now.
The purpose of the present study is to examine
whether the neo-classical transport process of
impurity ions is correctly incorporated in theimpurity transport code -”IMPGYRO” [2]. As a firststep, we have validated impurity cross field diffusion
in a uniform magnetic field in the IMPGYRO-codeby comparing with a classical theory of cross field
diffusion [3].
2. Cross-Field Diffusion of Charged Particles [3]We first briefly summarize the classical theory ofcross field diffusion based on Ref. 3. Consider now
two types of particle, which we denote by subscript
1 and 2. The general expression for the flux F1 ( X )of guiding centers of type 1 due to stochastic
process, i.e., Coulomb collisions of particles 1 and2, is given by,< ∆X > 1 ∂ [ N ( X ) < (∆X ) >],(1)
F = N (X )
−
2

1

11

∆t

12 ∂X

1∆t

where N ( X ) and ∆X are the density of the guiding
center of the particle species 1, and the step of theguiding center.
The step ∆X is related to the change in they-component of the velocity as
111(2)Under the uniform magnetic field Bz and theassumption that the both particle species aremagnetizedwiththeirLarmor
radii
rL1 = v1 / ωc1 , rL 2 = v2 / ωc 2 . Here, v i and ω ic = ei B z / mi( ei : electric charge, mi : mass) are the particle speed
and the cyclotron frequency of the i-th species,
respectively. For a guiding center of particle 1 at X ,the associated particle 1 is at x = X − v1 y / ω1c . The
density of particles 2 at x is the same as density ofguiding centers of particles 2 at x + v2 y / ω2 c
= X − v1y / ω1c + v2 y / ω2c . Hence, the probability per unit
time that a particle 1 with guiding center at X will
be involved in a collision, with scattering into solid
angle dΩ is given by∆X 1 = ∆v1 y / ω1c .

P(v1 , v2 , Ω)d 3 v1 d 3 v2 d Ω =(3)d 3 v1 d 3 v2 g1 (v1 ) g 2 (v) N 2 ( X + δ )vσ (Ω)d Ω,
where g1 , g 2 σ (Ω) , and v are the normalized

velocity distributions, the differential scattering crosssection of Coulomb collision, and the relative
velocity v = v1 − v 2 , respectively. In addition, δ isgiven by δ = v2 y / ω2 c − v1 y / ω1c from Eq. 2.
Using Eq. 3, we can calculate the first and thesecond moment of ∆X in Eq. 1 as follows:
12< ∆X >= αν r (1 +112L1e1e2)

1 ∂N

2N ∂X∆t ,(4)2(5)
Here, α is a numerical factor depending on themass m1 and m2 , and ν 12 is the collision frequency.
It should be noted that Eq. 5 gives the self-diffusionof the particle species 1, while Eq. 4 gives thetransport of the particle species 1 towards the density
gradient of the species 2.
22< (∆X ) >= 2αν r ∆t.
112 L 1Fig. 3 Impurity transport < ∆X > / ∆t due to backgrounddensity gradient as a function of impurity charge state.

Fig. 1 Time evolution of the average position < X >1

1

Fig. 2 Time evolution of the variance < ( X ) >2

1

Fig. 4 Self-diffusion coefficients < ( ∆X ) 2 > / ∆t incomparison with the theoretical value
1

3. Numerical Model

The trajectory of each test particle is followed by
Boris-Buneman algorithm [4]. In order to calculate
the velocity change ∆v due to Coulomb collision,
the Binary Collision Monte-Carlo (BCM) model has
been used, i.e., the scattering angle Θ and Φ aregiven in the following manner; Φ = 2π U ,Θ = 2arctanδ , where U is a uniform random numberand δ is sampled from the normal distribution withthe mean value < δ > and the variance < δ > [5] ,222< δ >= 0, < δ >=N 2 (e1 e2 ) ln Λ
8πε 02 M 2 v 3∆t c .

(6)

Here ln Λ is Coulomb logarithm, M is the reduced
mass of the species 1 and 2, ε 0 is the permittivity of
vacuum and ∆tc is the time-duration between two
successive binary collision. The particle density N
in Eq. 6 is evaluated at the actual position of eachtest particle.
24. Results and Discussion

Test simulations have been done under following
conditions as a reference case: 1) particle species 1:
tungsten (W) with charge state Z=1 as particle
species 1, 2) particle species 2: background
Deuterium (D), with the density scale length
N 2 /(dN 2 / dx) = −0.1 . 3) The initial velocity of thespecies 1and the velocity of background species 2
are sampled from the Maxwell velocity distribution
with a same temperature, and 4) uniform magnetic
field ( Bz = 5T ).
Figure 1 and 2 show the simulation results of< Xɶ > and < ( Xɶ ) > . The results are obtained by taking
an ensemble average over all the test particles( 105 particles) and normalized by Larmor radius.
2

1

1The theoretical values < Xɶ > / ν t in Fig. 1 and2
< ( Xɶ ) > / ν t in Fig. 2 are obtained from Eq. 4 andfrom Eq. 5, respectively.
Figure 3 shows the simulation results < ∆X > / ∆tfor the different charge state Z of the species 1.
Figure 4 compares the simulation results2< ( ∆X ) > / ∆t with the theoretical values for thedifferent values of the magnetic field ( B = 2T～5T ).
The results are plotted as a function of themagnetized parameter ω1cτ (= ω1c /ν 12 ) . The
numerical values of < ∆X > / ∆t and < ( ∆X ) 2 > / ∆tin Fig. 3 and 4 are obtained by a linear fitting of thetime evolution of < X > and < ( X ) > .
111

1

11

21

1

5. Conclusion and future study

As seen from Figs. 1-4, reasonable agreement has
been obtained between the numerical and theoreticalresults. The model improved in the present study will
be incorporated in the IMPGYRO code.
Neo-classical cross-field diffusions in realistic
tokamak magnetic configuration will be studied inthe future.
Reference[1] R.Dux et al., Plasma Phys. Control. Fusion 45,(2003) 1815.
[2] A. Fukano, M. Noritake, K. Hoshino, M. Toma etal.: J. Nucl. Mater. 363-365 (2007) 211-215.
[3] C.L.Longmire et al.: Phys. Rev. 103, 3, (1957) 508.
[4] C.K.Birdsall and A.B.Langdon, Plasma Phys.
Computer Simulation (McGraw-Hill, New York,1985).
[5] T.Takizuka et al.: J.Comp.Phys.25, (1977) 205.