TOPICSによるシミュレーションの現状
林伸彦、白井浩

原研那珂研、炉心プラズマ解析室

平成１５年３月１７-１８日、核燃焼プラズマ統合コード作業会、九大応力研

TOPICS(TOkamak Prediction and Interpretation Code System)
コンパイル : Fortran 90
実行環境 : IBM SP3, Compaq(Dec) alfa WS, PC Linux, SGI Origin
解析手法 : 定常（実験解析）
非定常解析（実験解析、シミュレーション）
実行モジュールは同じで、入力ファイルで切替え
図形出力 : Graphic library : GDlib + Eagle（原研オリジナル）
定常解析 : 計算終了後自動的に図形出力
非定常解析 : ファイル出力し、別のソフトToppltで図形出力
移植性 : 実験解析用ファイル入力と図形出力ルーチン以外は基本的に修正なしで可能
マニュアル : まとまったものはない（計画あるが実行できず）
Web上に資料をまとめ中各コードとのインターフェース :基本的に独立させようとしているが完全ではないTOPICS構成図
輸送係数異常輸送
Empirical, hi model

Sawtooth, Tearing,Ballooning

新古典輸送
C&H, H&H, H&S models,
NCLASS

Bootstrap電流
MI method, NCLASS

MHD現象

プラズマ輸送NBI/NBCD

1D transport

Stix, 1or2DFP

電流拡散

ECH/ECCD

1D diffusion

Ray tracing,
Relativistic FP

新古典抵抗

MHD平衡

H&H model, NCLASS

2DGrad-Shafranov

ペレット
Fehlberg, Parks models

放射損失

SOL/Divertor
1D transport or simple model

中性粒子
Monte-Carlo

不純物
1D transport

TOPICSの構成
プラズマ輸送(等・不等間隔mesh: 50-400)

トロイダルフラックスで定義された規格化小半径：r = (F(r)/F(1))0.5 (0 < r < 1)

∂n i ˆ∂n i∂ Ê2VrD=—¢Á˜ +Si∂r ¯∂tV ¢∂r Ë（多成分イオン可能）∂T j ˆ∂ Ê3∂ Ê2
ˆÁ n T ˜ =Á V ¢ —r n j c j
˜ + Pj∂t Ë 2 j j ¯ V ¢∂r Ë∂r ¯( j = e, i)

新古典拡散 :
(1) Chang & Hinton model(2) Hinton & Hazeltine model(3) Matrix Inversion法（Hirshman & Sigmar model）(4) NCLASS (Houlberg)
異常拡散 :(1) 経験則（大河、Kadomtsev）(2) Romaneli ηi model

MHD平衡 : MEUDAS (mesh: 128x128, 256x256, 512x512)
2D Grad-Shafranov eq.
Fixed or Free boundary, 固有値 or FCT解法
電流拡散式¸∂ Ï ∂ Ê ∂Y ˆ∂ Ê ∂Y ˆ=rDESÁ˜Á˜Ì˝∂t Ë ∂F ¯ ∂r Ó ∂r Ë ∂F ¯˛

新古典抵抗 : Hirshman & Hawryluk model (or NCLASS)
Bootstrap電流 :Matrix Inversion法（係数はHirshman & Sigmar model）(or NCLASS)NBI : Stix定常解 or 1D or 2D Fokker-Planck eq.
NBCD : Mikkelsen model (Stix & 1DFP)
Polevoi (ASTRA) model (2DFP)
ECH/ECCD : Ray tracong and Relativistic Fokker-Planck eq.

中性粒子 : Monte-Carlo法
不純物輸送 : Full charged impurity ionor IMPACT (fluid eqs. + NCLASS)
放射損失 : Given profile or Calculated from carbon Lz(Te)or IMPACT
ペレット : Fehlberg model or Parks model (including ExB drift)
MHD現象 : Sawtooth : Kadomtsev model
Tearing mode : Cylinder model for stability index D'
Modified Rutherford eq. for island width
Bollooning mode : Stability analysis
SOL / Divertor : IMPACT or five-point model (予定)Modified Rutherford equation
Lq 2 Ê 1 ˆ—rWm 0 dW22 12—kebr= k1 ( D¢ - aW ) —r + k2m 0 Lq jBS1Á˜B p W 2 + W d 2 3 s p rsL p Ë q 2 ¯Wh dtClassicalÊ r pi Lq ˆ-k e b p Á˜Ë Lp ¯
1.54 sBootstrap
2—r2GGJLq —r
1Wkmj
h5 0rs B p EC EC W 2
W 4 + W p4PolarizationEC

Lq = q / ( dq/dr) , rs : Rational surface position
W : Magnetic island width in r coordinate, D'(W) : Cylindrical model, a = (m/rs)2
Wd : Finite c^/ c // effect (R.Fitzpatrick, Phys. Plasmas 2(1995)825.)
Wp : Validity limit of the model, Wp = es0 .5rpi assumed (H.Reimerdes, LRP 700/01(2001).)
jEC : Averaged EC current, jEC = IEC/pa2
hEC : Efficiency calculated on assumed island structure(C.C.Hegna and J.D.Callen, Phys. Plasma 4(1997)2940.)

hEC =dacos( ma ) jEC
2pdardÚ Ú 2p jECÚ drÚECH/ECCD code (K.Hmamatsu, et al., Plasma Phys.Control. Fusion 42(2000)1309.)
Poloidal plane

EC ray trajectory : Standard ray tracing method

2Profiles of driven current and heating :Relativistic Fokker-Planck eq.
Example :Fundamental O-mode wave of 110 GHz
No modulation in phase with island rotation
Poloidal location

Z (m )

1-1

100.450.6
rBeam divergencejEC (MA/m2)

0.5

Peakvalue
0°jEC0WEC
0.45

r

2°
4°
0.6

jEC0 (MA/m2)

jEC (MA/m2)

0.4

0000.03

rEC ~ rs0 =0.515

4 MW
IEC ~130 kA
93

3260

0.05

WEC

JT-60SC-2

1234

R (m )

Toroidal plane

20°R0.07Ip

IMPACT (IMPurity Analysis Code for Tokamak plasma)
Fluid transport eqs. + 新古典拡散係数(NCLASS)
解析例： 定常解析、Ar0.
40.
0 51*1-0.
00.00 0.
1*2

Ele.
Ion

Ele.
Ion
fast0.
4)
3* 0* .
m 3/02** 001 .
( 2NED0.
12-0.
3)
3**m/0W .
M 1(WOP)
Vek0
(.
2MET
0.
1

n0.2

Prad

T
0.4

0.6
RHO

0.8

0.
00.0

1.0

0.2

0.4

0.6
RHO

0.8

0.
01Q
J//
BSC
RFCD
NBCD
ECCD

0.
8

jX0A.
M6
J/Jq0,.
4Q0.
20.
00.2

0.4

0.6
RHO

0.8

1.00.
00.0

1.0

5-0.0NBIENBIIJOULRADECHRFERFI0.2

0.4

0.6
RHO0.8

1.04-01*0) .
3 2**m/02**01( 0y .
1tisneD0.
001*chrg 1chrg 2chrg 3chrg 4chrg 5chrg 6chrg 7chrg 8chrg 90.
)23**m/02**010(.
y1tisneD

nimp
9.00.
0
0.0

10.0
RHO

*10 -1chrg10chrg11chrg12
chrg13chrg14chrg15chrg16chrg17chrg18nimp

0.2

0.4

0.6
RHO

0.81.0ダイバータプラズマ５点モデル( Ref. : N.Hayashi, et al., J.Nucl.Mater. 266-269 (1999)526. )Stagnation point

Stagnation point

ThroatDivertor A
Inside plate
Heat flux

SOL
current

ThroatDivertor B
Outside plate(G0=0)SOL

Divertor region B
Throat B
Plate BCore

Divertor region A
Throat A

jPlate A

SOL

SSOL WSOL

Sdiv Wdiv
Ldiv

lstag

LSOL-lstag

Sdiv Wdiv
Ldivz(along B)
Particle &energy
source/sink∑ S ystem of model is divided into four regions, i.e., two SOL regions andtwo divertor regions.
∑ Fluid equations are integrated along B in each region and reduced to a setof nonlinear algebraic equations with physical variables at the five positions.
∑ All the quantities are those in a magnetic flux tube the nearest to the separatrix.

５点モデル結果例

Asymmetric equilibrium with current
31023
Ts=10eV1022

2040801022

1021

1Qsep (MW)Dynamic response to ELM particle/heat pulsew/
momentum loss

1A
01020

5030

R

40

50

T0Ts

040n/nnor

0.1

25
R 010 20Ts/DE

R = 1 2 5 10 201023Fsep (s-1)

Fsep (s-1)

1024

w/o momentum loss

2T/DE

High/Low recycling (symmetric) equilibrium

Bnsn0010-4

10-2t /tnor

160シミュレーション例（Tearing stability, sawtooth）
Parameters :I p = 1.5 MA, BT = 3.5 T, R = 3.2 m, a = 0.8 m, V = 60 m3 , k = 1.5, d = 0.1
n e = 1.0 ¥ 1019 m-3 , Te ª 1.8 keV, Ti ª 1.3 keV, b p ª 0.1, W s ª 0.4 MJ, POH ª 0.7 MW(Ohmic plasma)

Sawtooth : Kadomtsev modelÚ ( 1 q - 1) r dr = 0rc

1

30fdV = fAfterrc

Ú dV

r 2/1

q axis0Density and pressure profiles within 0<r<rc
are redistributed. ( f Before → f After )Sawtooth period : 200 ms

1

Tearing stability index D' : Cylinder model2
3 MA/m2

r 3/2

0.5

r

0Beforeq axisÚrc

Sawtooth
040

05jq

42/1

3D'20

t=0s

8s
1

0
0.0

6s
4s2s2

t=0s

1

0.2

0
0.0

0.4

0.6

r

0.8

1.02s4s
6s
8s

0.2

0.4

r

3/2

0.6

0.8

1.0

-40024time (s)68シミュレーション例（Tearing stability）
NBI before sawtoothm/n=3/2

0.6

1.5Sawtooth

D'020

NBI 15MW

0-40
0510j (MA/m2)

rsrs

NBI after sawtooth
0.62rs

SawtoothNBI

D'

050-150510time (s)15j (MA/m2)

10 MW
20 MW

t = 3.5 s(NBI start)

0015time (s)

5.6 s

jjBS
r

5 MW

1t = 8.7 s

j1 30 MW

rsjBS
00r1

シミュレーション例（新古典テアリングモード）
Typical JT-60SC Parameters :Ip=3.0 MA, BT=3.8 T, R=2.9 m, a=0.85 m, PNBI=30 MW, ne=4.2×1019 m-3,
Te=6.4 keV, Ti= 7.5 keV
Diffusivities : Empirical model
0.2
Di = c e =ci= 0.06(1 + 2 r 2 )(1 + PNBI 0.5 ) m2 / s2W

NBI source & Radiation : Fixed profile
Flattening effect : Amplified diffusivities inside islandNBI
30 MW

BS current modifies profile. NTM is triggered by D'>0.

02

1.5
1

t=3s

40

7s0

rs0

2

468 10 12 14

time (s)

00

D'(0)

D'(W)
rsTotalBootStrap

0.5

0

1D'bN2r

1-800rj (MA/m2)

2.50
48time (s)12シミュレーション例（新古典テアリングモード）

6 12s
4200

16121214st=2sr

840100.6r0rr

1r143210145qj (MA/m2)
080061jBS (MA/m2)

2.5Ti (keV)

8s

816Te (keV)

ne (1019m-3)

100r

100

シミュレーション例（ECCDによるNTM安定化）
NTM is fully stabilized by EC current on the island center.
Off-center EC current can decrease island width, but not fully stabilize NTM.
EC current moves rational surface through background current profile modification.
→ Decrease stabilizing efficiency of EC

j (MA/m2)

12

EC

EC 2 MWIEC ~ 60 kA

30.2
Bootstrap
at 10 s
00.4

0.5 r

0.6

0.54

D'>0D'>0EC
2MW

0
1.9

12

2

3
0.5

w/o ECbN

rsrEC 3

0.16

W

0.4

1EC 2 MW
10time (s)w/o EC
1.4

118

10

12time (s)14TOPICSシミュレーションの課題
・異常輸送モデル導入
・ダイバータモデル導入（５点モデル）
・マニュアル・各コードの独立性とインターフェース（ユニット化）
・実行環境に依存しない汎用性（主にGraphic lib.の問題）
・ユーザー不足
・プログラム改良（goto削除、implicit none等）