スライド 1 - JT

第18回若手科学者によるプラズマ研究会
2015/03/04-06
SOL-Divertor Plasma Simulations by
Introducing Anisotropic Ion Temperatures and
Virtual Divertor Model
非等方イオン温度と仮想ダイバータモデルを導入
したSOL-ダイバータプラズマシミュレーション
Satoshi Togo, Tomonori Takizukaa, Makoto Nakamurab,
Kazuo Hoshinob, Kenzo Ibanoa, Tee Long Lang, Yuichi Ogawa
Graduate School of Frontier Sciences, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-8568, Japan
aGraduate school of Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japan
bJapan Atomic Energy Agency, 2-166 Omotedate, Obuchi-aza, O-aza, Rokkasho 039-3212, Japan
2
Motivation
The SOL-divertor plasma code packages (SOLPS, SONIC, etc.)
・used to estimate the performance of the divertors of future devices
・some physics models are used in the plasma fluid model (e. g. viscosity)
・physics models are valid in the collisional regime
Result from PARASOL code
parallel momentum transport equation (1D)


mi nV 

V 

mi nV 2  nTi  nTe     i , //
  Mm
t
s
s 
s 
The parallel ion viscous flux -ηi,//(∂V/∂s)
・the approximated form of the stress tensor p
(π = 2n(Ti,// - Ti,⊥)/3)
・derived under the assumption that π << nTi.
The kinetic simulations showed a remarkable
anisotropy in the ion temperature even for the
medium collisionality.
collisional
collisionless
A. Froese et al., Plasma Fusion Res. 5 (2010) 026.
The boundary condition Mt = 1 has been used in the conventional codes.
However, the Bohm condition only imposes the lower limit as Mt ≥ 1.
Momentum Eq. & Virtual Divertor Model
3
Introduction of the anisotropic ion temperatures, Ti,// and Ti,⊥, to the fluid model
・changes the momentum transport equation into the first-order
・makes the explicit boundary condition at the divertor plate unnecessary
Parallel-to-B component of the Boltzmann equation

f
f eE f
f
 v//


t
s m v// t

mi nV 

mi nV 2  nTi , //  nTe  M m
t
s
(isotropic Te is assumed)

p
2




mi nV

mi nV  nTi  nTe 
 Mm
t
s
s
nTi  nTi , //  2Ti ,  3
p  2nTi , //  Ti ,  3
mi nV 

V 

mi nV 2  nTi  nTe     i , //
  Mm
t
s
s 
s 
coll
nTi, //  nTi  p
nTi,   nTi  p 2
conventional codes
(effective isotropic Ti)
Instead of the boundary condition Mt = 1, we modeled the effects of the divertor
plate and the accompanying sheath by using a virtual divertor (VD) model.
Flow velocity is not determined by downstream
‘waterfall’ but by upstream condition.
P. C. Stangeby, The Plasma Boundary of Magnetic Fusion Devices.
4
Plasma Fluid Eqs. & Artificial Sinks in VD
Artificial sinks in the virtual divertor (VD) region
Eq. of continuity
n nV

S
t
s
S VD  
Eq. of momentum transport
mi nV 

mi nV 2  nTi , //  nTe  M m
t
s


n
M mVD  
 VD
QiVD
, //  
mi nV
 VD

1 1
1

2
m
nV

g
nTi , // 

i
i , //
VD
2
 2

nTi ,  Ti , // 
 rlx

me nTe  Ti , // 
nTe
V
mi
e
s
t

nTi ,V
s
 1  c 
qieff
,
s
 Qi , 
nTi ,  Ti , // 
 rlx

2me nTe  Ti , 
mi
e
Eq. of electron energy transport
qeeff
3m nTi  Te 
nTe
 3
  5

 Qe  e
V
 nTe    nTe V 
t  2
s
mi
e
s
 s  2

S. Togo et al., J. Nucl. Mater. (2015) in Press.
1  3

g
nT


e
e
 VD  2

Periodic boundary condition
Eq. of perpendicular (⊥) ion energy transport
nTi ,
QeVD  
※according to the
image of a waterfall
Eq. of parallel (//) ion energy transport
qieff
 1
1
3
  1

, //
2
2
 mi nV  nTi , //    mi nV  nTi , // V  c
t  2
2
2
s
 s  2

 Qi , // 
g i , nTi ,

VD V 


 mi nDm
 QiVD
,
s 
s 
 VD
Plasma Fluid Eqs. & Artificial Sinks in VD
Eq. of continuity
n nV

S
t
s
Eq. of momentum transport
mi nV 

mi nV 2  nTi , //  nTe  M m
t
s

Because the parallel internal energy convection is
3 times as large as the parallel internal energy, Ti,//
becomes lower than Ti,⊥.

Ion pressure relaxation time rlx = 2.5i.
(E. Zawaideh et al., Phys. Fluids 29 (1986) 463.)
Eq. of parallel (//) ion energy transport
qieff
 1
1
3
  1

, //
2
2
 mi nV  nTi , //    mi nV  nTi , // V  c
t  2
2
2
s
 s  2

 Qi , // 
nTi ,  Ti , // 
 rlx

me nTe  Ti , // 
nTe
V
mi
e
s
c = 0.5 is used.
Eq. of perpendicular (⊥) ion energy transport
nTi ,
t

nTi ,V
s
 1  c 
qieff
,
s
 Qi , 
Heat flux limiting factors;
ai,// = ai,⊥ = 0.5, ae = 0.2
are used.
qseff = (1/qsSH + 1/asqsFS)-1
nTi ,  Ti , // 
 rlx

2me nTe  Ti , 
mi
e
Eq. of electron energy transport
qeeff
3m nTi  Te 
nTe
 3
  5

 Qe  e
V
 nTe    nTe V 
t  2
s
mi
e
s
 s  2

Neutral is not considered
at first.
5
Results (Anisotropy of Ti and the Profiles)
Anisotropy of Ti vs normalized MFP of i-i collision
VD
n (1019 /m3)
collisional case
M
weakly collisional case
Ti,// (eV)
Ti,⊥ (eV)
Te (eV)
Gsep and Psep are changed.
6
Results (Bohm condition
M   V cs
cs 
Te   a Ti , //
M*
= 1)
Time evolutions of M* for various VD
mi
a = 3 for adiabatic, collisionless sound speed
PARASOL
T. Takizuka and M. Hosokawa, Contrib. Plasma Phys. 40 (2000) 3-4, 471.
M* saturates at ~ 1 independently of the value of VD.
Characteristic time Bohm depends on VD and 2 orders shorter than the quasi-stationary time
~ 10-3 s.
For the simulations of transient phenomena, such as ELMs, Bohm has to be smaller than
their characteristic times.
7
Results (Sheath heat transmission factors)
Boundary conditions for the heat flux at the divertor plate in the conventional codes;
5
T
1

qi   mi nV 2  nTi V   i , // i   i nTiV
2
x
2

From the sheath theory,
i  3
 me  Ti 
1  
 e  2.5  0.5 ln  2p
mi  Te 

qe 
T
5
nTeV   e, // e   e nTeV
2
x
For Ti = Te,
e ≈ 5 for H+ plasma
e ≈ 5.3 for D+ plasma
Relation between sheath heat transmission factors and gs
i scarcely depends on gs because convective heat flux dominates conductive one.
e can be adjusted to the values based on the sheath theory by VD model.
8
Results (Dependence of the profiles on
Decay length in VD region: Ld ~ VtVD.
As long as Ds < Ld < LVD, the profiles in
the plasma region do not change.
If Ld > LVD (when VD = 5×10-5 s in
figure), the profiles become invalid.
If Ld < Ds, numerical calculation diverges.
The profiles just in front of the divertor
plates are affected by the artificial sinks in
VD region due to numerical viscosity.
This problem will be solved by
introducing a high-accuracy difference
scheme and an inhomogeneous grid.
VD)
9
Results (Supersonic flow due to cooling)
C = (RG/Rp)(Tt/TX)1/2, T = Ti,// + Te(,//)
M = V/cs
cs 
Te (, //)  Ti , //
mi
Isothermal sound speed
PARASOL
(at the plate)
(at the X-point)
Cooling term Qe = -nTe/rad is set in
the divertor region.
RG = Rp = 1
Mt well agrees with the theory.
The reason for MX > 1 is under
investigation.
T. Takizuka et al., J. Nucl. Mater. 290-293 (2001) 753.
10
Results (Ion Viscous Flux vs Stress Tensor)
Two kinds of ion viscous flux,
p
BR

i , //
  i , //
V
x
 0.96nTi i 
p lim
 1
1 

  BR 
nTi 
p
1
 = 0.7 :viscosity limiting factor
are compared to the stress tensor,
pdef = 2n(Ti,//-Ti,⊥)/3, in the particle-sourceless-region (s = 20.005 m) and particlesource-region (s = 17.005 m).
pBR becomes 2~3 orders larger than pdef as
lmfp/L becomes large.
In the particle-source-region, the correlation between plim and pdef becomes worse
especially in the collisional region.
In the particle-source-less region, plim with  = 0.7 comparably agrees with pdef. However, 
depends on the anisotropy of ion pressure which might change with the neutral effects.
Therefore, it is necessary to distinguish between Ti,// and Ti,⊥.
11
12
Self-Consistent Neutral Model (in VD)
present
conventional
Recycling neutral
Diffusion neutral
The coordinate x
: poloidal direction x = (Bp/B)s.
boundary condition
Recycling neutral (inner plate)
nn,inrecy
t


 nn,inrecyVn,inrecy
x
n nV

S
t
s
Diffusion neutral
 
S VD  
nn,diff
t
Recycling neutral (outer plate)
n
0
periodic boundary condition
 VD
n
 VD
 2
nn,diff
nn,outrecy
 n,VDdiff
t


 nn,outrecyVn,out
recy
x
 
n
0
 VD
(Eq. of continuity for plasma)
nn ,diff
 

   Dn
x 
x
n

   n,VDdiff
 n,diff

n,diffVD : input
 2
nn,diff
 n,VDdiff
13
Self-Consistent Neutral Model (in Plasma)
Recycling neutral (inner plate)
nn,inrecy
t


 nn,inrecyVnin,recy
x
n
in
n, recy
 n,recy
in
in
 S iz,
recy  S cx, recy
where VFC = (2εFC/mi)1/2 with
Franck-Condon energy εFC = 3.5 eV.
Recycling neutral (outer plate)
nn,outrecy
t


 nn,outrecyVn,out
recy
x
n
Diffusion neutral
nn,diff
n
 
   Dn n ,diff
t
x 
x
Dn 
out
n, recy
 n,recy
in
Vn,out
recy  Vn, recy  VFC 2
out
out
 S iz,

S
recy
cx, recy
n

   n,diff  S iz,diff  S rc  S cx,recy
 n,diff

 n,recy 
1
 n,recy
 n,diff 
1
 n,diff
  m 
 a L  FC i 


d


 T m
 aL  i i
 d

aL : input




Ti
mi  cx   iz   n,diff 
T. Takizuka et al., 12th BPSI Meeting, Kasuga, Fukuoka 2014 (2015).
The coordinate x: poloidal direction x = (Bp/B)s.
14
Atomic Processes
Recycling neutral (outer plate)
Recycling neutral (inner plate)
nn,inrecy
t


 nn,inrecyVnin,recy
x
n
in
n, recy
 n,recy
nn,outrecy
in
in
 S iz,
recy  S cx, recy
t


 nn,outrecyVn,out
recy
n
out
n, recy
 n,recy
x
out
out
 S iz,
recy  S cx, recy
Diffusion neutral
nn,diff
t

n
 
  Dn n ,diff
x 
x
n

   n,diff  S iz,diff  S rc  S cx,recy
 n,diff

Source terms for plasma:
S  S core  S iz  S rc


out
in
out
in
M m  mi 2 VFC sin  S iz,
recy  S iz, recy  S cx, recy  S cx, recy  miV S cx  S rc 

  T S
 
  T S
(θ = Bp/B)

2
2
Qi , //  Qi,core
//  Ti 2 S iz,2  S cx,2   mi 6 VFC S iz, recy  S cx, recy  mi 2 V  Ti, // 2 S cx  S rc 

2
Qi,  Qi,core
  mi 3VFC S iz, recy  S cx, recy
Qe  Qecore   iz S iz  3 2Te S rc
i
iz,diff
(εiz = 30 eV)
 S cx,diff
i, 
cx
 S rc 
(Ti = (Ti,// + 2Ti,⊥)/3)
15
Result (Low recycling condition)
X-point
Near the plate
aL = 1
Recycling rate ~ 0.17
Ti,///Ti,⊥ ~ 0.6
Recycling neutral dominant
Ti,⊥
Ti,//
Te
nn,recy
nn,diff
16
Result (High recycling condition)
X-point
Near the plate
aL = 0.1
Recycling rate ~ 0.92
Ti,///Ti,⊥ ~ 1
Diffusion neutral dominant
Ti,// Ti,⊥
Te
nn,diff
nn,recy
Conclusions
17
1D SOL-divertor plasma model with anisotropic ion temperatures has been
developed. In order to express the effects of the divertor plate and the accompanying
sheath, we use a virtual divertor (VD) model which sets artificial sinks for particle,
momentum and energy in the additional region beyond the divertor plate. In addition,
VD makes the periodic boundary condition available and reduces the numerical
difficulty.
For simplicity, the symmetric inner/outer SOL-divertor plasmas with the
homogeneous magnetic fields are assumed. In order to simulate more general
asymmetric plasmas with the inhomogeneous magnetic fields, the effects of the
plasma current and the mirror force have to be considered. In addition, it is
necessary to introduce a high-accuracy difference scheme and an inhomogeneous
grid in order to avoid the numerical errors at the divertor plate. These are our future
works.
DmVD & gi,// in VD region
DmVD and gi,// in VD region have Gaussian shapes.
The length of V-connection-region L0 = 1.6 m.
DmVD
 cD
g i , // 
g i,c//
L20
 VD
 4s 2 
exp  2 
 L0 
 4s 2 
exp  2   g i,e//
 L0 
7
Results (Bohm condition)
cs

Te , //    a Ti , //
Mach profiles for various Ds
mi
a = 3 for adiabatic, collisionless sound speed
M* ≈ 1 with no cooling effects.
Plasma
T. Takizuka and M. Hosokawa, Contrib. Plasma Phys. 40 (2000) 3-4, 471.
VD
The effect of artificial sinks in VD region
numerically diffuses in the plasma region.
Appendix ~ Collisionless Adiabatic Flow ~
1D equations in the collisionless limit;
d
nV   0
ds


d
mi nV 2  nTi , //  nTe  0
ds
dnTe
d 1
3

3
 mi nV  nTi , //V   V
ds  2
2
ds

 2 3Ti , //  Te  dn 3n dTe
V 


m
ds
mi ds
i


Refer to Sec 10.8 of Stangeby’s text
20
The effect of gs on t (heat transmission factor)
The boundary condition for the heat flux at the divertor plate in the usual codes;
5
T
1

qi   mi nV 2  nTi V   i , // i   i nTiV
2
x
2

qe 
T
5
nTeV   e, // e   e nTeV
2
x
The heat transmission factors, i and e, are input parameters.
The VD model, however, does not use this boundary condition but the periodic boundary
condition with the cooling index gs (s ∈ i//, i⊥, e). Therefore i and e are back calculated
using these relations.
The conduction heat fluxes are limited by the free-streaming heat fluxes with limiting
coefficients as as qseff = (1/qsSH + 1/asqsFS)-1.
qsSH   sSH
Ts
x
qiFS
, //  nTi, //
Ti , //
mi
qiFS
,   nTi, 
Ti , //
mi
qeFS  nTe
Te
me
Thus the effective conduction heat fluxes are smaller than the free-streaming heat fluxes
times limiting coefficients asqsFS so that t has the maximum.
 9f 6 1
M 2  3 f ani

i 
 f eq   ani

2  f ani  2
2
f

4
M
ani

5 a
e   e
2 M
mi
me
 f ani  2 f eq
3 f ani   f ani  2  f eq

3 f ani
3 f ani
3 
 ca i , //

 1  c a i ,
f ani  2
f ani  2  3 f ani   f ani  2  f eq

f ani 
Ti , //
Ti ,
f eq 
Te
Ti
Calculation condition
Calculation condition
H plasma and ni = ne = n
Symmetric inner/outer SOL
Length of the plasma L
44 m
SOL width d
2 cm
Separatrix area
40 m2
Particle flux from core Γsep
1~5×1022 /s
Heat flux from the core Psep
1~4 MW
Cooling index for i,//
1
Cooling index for i,⊥
1.2
Cooling index for e
2.5
VD
5×10-6 s
Heat flux limiter for ion
0.5
Heat flux limiter for electron
0.2
Discrepancy
Edge transport code packages, such as SOLPS and SONIC, are
widely used to predict performance of the scrape-off layer (SOL) and
divertor of ITER and DEMO. Simulation results, however, have not
satisfactorily agreed with experimental ones.
Comparison of results (EXP vs SIM)
M. Wischmeier et al., J. Nucl. Mater. 390-391, 250 (2009).
Why does Ti,// become lower than Ti,⊥?
Integration over x from the
stagnation to x yields,
Reduced eq. of parallel (//) ion energy transport
d 3
1

 nTi , //V   Qi , //  Qi
dx  2
3

Ti , // 
2Qi x
9nV
Ti , //
1
3
Reduced eq. of perpendicular (⊥) ion energy transport
dnTi ,V
dx
2
 Qi ,  Qi
3
Ti ,
By considering the kinetic energy term and force term, Ti,///Ti.⊥ ~ 0.2.
Eq. of parallel (//) ion energy transport
qi , //
 1
1
3
  1

2
2
 mi nV  nTi , //    mi nV  nTi , // V  c
t  2
2
2
x
 x  2

nTi ,  Ti , //  me nTe  Ti , // 
nTe
 Qi , // 

V
 rlx
mi
e
x
eff
~
Ti , 
2Qi x
3nV
Qualitative derivation of the viscous flux
Assumption of p << nTi
Simplified system equations
d
nV   0
ds
d
mi nV  nTi  p   0
ds


(A)

(E)’
(B)
d
3p
mi nV 3  3nTi , //V  
ds
 rlx
(C)
d
nTi,V   3p
ds
2 rlx
(D)
d
mi nV  nTi   0
ds


Then
(E)
(B)’
By (A) and (B)’, LHS of (E)’ becomes
d
dV
mi nV 3  2nTiV  2nTi
ds
ds
From (C) – (D)
d 
7  
9p

 mi nV 3   2nTi  p V   
ds 
2  
2 rlx


d
9p
mi nV 3  2nTiV  
ds
2 rlx
4
9
p   nTi rlx
dV
dV
  i , //
ds
ds
Necessity of artificial viscosity term
M VD  
mi nV
 VD
 
V 
 VD

x 
x 
artificial viscosity term

conservation of ion particles
conservation of parallel plasma momentum
d V 


dx



V
d
2
V  2nT  
dx

V
2
 c s2

2
dcs2
dV c s

V
dx

dx
When V is positive, RHS becomes positive. If V becomes supersonic,
dV/dx becomes positive and V cannot connects.
Discretization
discretization scheme
general conservation equation
  
 
 
 V     G
S
t
x
x 
x 
full implicit
upwind
staggered mesh (uniform dx)
central
Calculation method
matrix equation
Gx  s (ex. N = 5)
 ae2
0
0
 a p1

 ae2
0
  a w2 a p 2
 0
 a w3 a p3
 a e3

 0
0
 a w4 a p 4

0
0
 a w5
  a e5
  a w1 


 0 
u 0 


 0 



a
 w5 
1
 
0
v  0
 
0
 
1
 a w1  1   s1 
   
0   2   s 2 
0   3    s 3 
   
 a e 4   4   s 4 

a p5   5   s 5 
 a p1  a w1

  a w2
0
A  

0

0

Matrix G becomes cyclic
tridiagonal due to the
periodic boundary condition.
This
matrix
can
be
decomposed by defining
two vectors u and v so that
G  A  uv where A is
tridiagonal.
 a e2
0
0
a p2
 a e2
0
 a w3
a p3
 a e3
0
 a w4
a p4
0
0
 a w5


0


0

 ae4 

a p5  a e5 
0
Calculation method
Sherman-Morrison formula
A  uv
1
1
1

1
 A  A u 1 v  A u
zv 

x  G s  I 
y
 1 v  z 
1

1
vA1
where Ay  s and Az  u.
y and z can be solved by using tridiagonal matrix
algorithm (TDMA).
calculation flow
No
Ion // energy
Ion ⊥ energy
Elec. energy
Momentum
Particle
Yes
The number of equations can be changed easily.
Continuity of Mach number
conservation of ion particles
conservation of parallel plasma momentum
d nV 
S
dx

ncs 1  M 2


d
mi nV 2  2nT  0
dx





2
nc
M
1

M
dM
dT
 1 M 2 S  s
dx
2T
dx
O. Marchuk and M. Z. Tokar, J. Comput. Phys. 227, 1597 (2007).
Due to the continuity of Mach number, RHS has to be zero at the sonic
transition point (M = 1).
RHS > 0
RHS ≦ 0
Sonic transition has to occur at the X-point when T = const.
Result (particle flux & Mach vs nsep)
Γt ∝ nsep
→ accords with conventional
simulations
Supersonic flow (Mt > 1)
→ observed when nsep is low
Subsonic flow (Mt < 1)
→ observed when nsep is high
→ numerical problem?
Larger nsep (like detached plasmas)
is future work.
Result (Mach number near the plate)
Near the plate
Plasma
Δs = 5cm
Δs = 2cm
M > 1 is satisfied in the near-plate VD region.
Smaller Δs results in a better result. → Numerical problem?
VD
Result (Mt vs nsep)
The recycling neutrals are not
ionized or do not experience the
charge exchange near the plate (red
line).

Download Report