移送磁気プラズマ中性化過程のハイブリッド・シミュレーション

平成２４年度 修 士 論 文

移送磁化プラズマ中性化過程のハイブリッド・シミュレーション

指導教員 髙橋 俊樹 准教授

群馬大学大学院工学研究科

電気電子工学専攻

渡辺 隆之

目次

第 1 章 序論･･････････････････････････････････････････････････････････1

1.1. 研究背景･･･････････････････････････････････････････････････････1

1.2. 核融合発電･････････････････････････････････････････････････････3

1.3. Compact Torus プラズマ･･･････････････････････････････････････････5

1.4. 超高速フロー入射実験･･･････････････････････････････････････････7

1.5. 研究目的･･･････････････････････････････････････････････････････8

第 2 章 シミュレーション・モデル･･････････････････････････････････････9

2.1. 平衡計算･･･････････････････････････････････････････････････････9

2.1.1. CT パラメータ･･･････････････････････････････････････････････9

2.1.2. Grad-Shafranov 方程式････････････････････････････････････････11

2.1.3. プラズマ・パラメータ････････････････････････････････････････13

2.2. ハイブリッド・シミュレーション･････････････････････････････････15

2.2.1. イオンの初期分布･･･････････････････････････････････････････15

2.2.2. 粒子軌道計算･･･････････････････････････････････････････････17

2.2.3. 磁場と電場の時間発展･･･････････････････････････････････････18

2.2.4. クーロン衝突･･･････････････････････････････････････････････19

2.2.4.1. Slowing-Down Collision･･･････････････････････････････････19

2.2.4.2. ピッチ角散乱･･･････････････････････････････････････････20

2.2.5. PIC 法･･････････････････････････････････････････････････････23

2.2.5.1. 勾配補正フィルタ･･･････････････････････････････････････23

2.2.6. 平滑化･････････････････････････････････････････････････････27

2.3. 中性化････････････････････････････････････････････････････････28

2.3.1. 軸方向 NBI モデル･･･････････････････････････････････････････28

2.3.2. 中性流体･･･････････････････････････････････････････････････29

2.3.3. 中性化判定･････････････････････････････････････････････････30

2.3.4. 反応速度係数･･･････････････････････････････････････････････31

2.3.5. 低速イオン･････････････････････････････････････････････････33

2.4. 高速並列演算環境･･････････････････････････････････････････････34

2.4.1. ハードウェア･･･････････････････････････････････････････････35

2.4.2. ソフトウェア･･･････････････････････････････････････････････35

2.4.3. 演算能力（姫野ベンチマークテスト）･･････････････････････････36

2.4.4. 同期通信と並列化･･･････････････････････････････････････････37

2.5. FRC プラズマの場合･････････････････････････････････････････････39

2.5.1. 平衡解･････････････････････････････････････････････････････39

2.5.2. プラズマ・パラメータ････････････････････････････････････････40

2.5.3. 電場･･･････････････････････････････････････････････････････41

第 3 章 シミュレーション結果･････････････････････････････････････････42

3.1. 演算結果（Spheromak の場合）････････････････････････････････････42

3.1.1. イオン密度の推移･･･････････････････････････････････････････42

3.1.2. CT 中性化率･････････････････････････････････････････････････46

3.1.3. 磁場・電場の推移････････････････････････････････････････････47

3.1.4. 電流減衰･･･････････････････････････････････････････････････51

3.2. 演算結果（FRC の場合）･･････････････････････････････････････････55

3.2.1. イオン密度の推移･･･････････････････････････････････････････55

3.2.2. CT 中性化率･････････････････････････････････････････････････59

3.2.3. 磁場・電場の推移････････････････････････････････････････････60

3.2.4. 電流減衰･･･････････････････････････････････････････････････69

3.3. Spheromak と FRC での中性化挙動差･･･････････････････････････････72

3.3.1. 電場モデルの差･････････････････････････････････････････････72

3.3.2.

J



B





p

の比較･･･････････････････････････････････････････73

第 4 章 結論･････････････････････････････････････････････････････････79

4.1. まとめ････････････････････････････････････････････････････････79

4.1.1. Spheromak の場合････････････････････････････････････････････79

4.1.2. FRC の場合･････････････････････････････････････････････････79

4.2. 今後の課題････････････････････････････････････････････････････80

謝辞････････････････････････････････････････････････････････････････81

参考文献････････････････････････････････････････････････････････････82

1 / 83

第 1 章 序論

1.1. 研究背景

さる，2011 年（平成 23 年）3 月 11 日に発生した東北地方太平洋沖地震及びその後の余震 により引き起こされた大規模地震災害はまだ記録に新しい．特に，発電所の緊急停止に伴 う停電，津波による冷却系統遮断によって連鎖的に発生した福島第一原子力発電所事故な どは，節電要請や立ち入ることのできない警戒区域・避難指示区域という形で 2013 年現在 でも影響が継続している． 上記の福島第一原子力発電所事故により，日本国のエネルギー政策は大きな軌道修正を 迎えつつある．1955 年 12 月 19 日に成立した原子力基本法をもとに，我が国は原子力発電 の比率を徐々に増やし，大震災直前では京都議定書[1]_{の「二酸化炭素排出量を 2012 年まで} に 1990 年比で平均マイナス 6%」を目標に原子力発電のさらなる増強を目指していた．し かし，大震災後の世論変化及び原子力発電のコスト見直し等により，新規原子力発電所の 建設は困難になり，既存発電所の再稼働もストレス・テストや破砕帯などの再調査の必要 性によって見通しが立っていない[2]_{のが現状である．} Fig. 1.1 日本の電力発電構成の将来展望図．（2012 年度エネルギー白書 [3] より．）

2 / 83 Figure 1.1 にあるように，2010 年度実績における原子力発電の割合は約三割に迫る勢いで あり，京都議定書の二酸化炭素排出量を 2012 年までに 1990 年比で平均マイナス 6%にする という目標達成には，火力発電の割合を減らし，原子力発電を 45%まで増強する必要があ るという算定であった．（2012 年度におけるエネルギー基本計画[3] より．） しかし，原子力発電の見直しを求める声が無視できない中，仮に原子力発電の割合を段 階的に 0%，15%，20～25%にそれぞれ減らす，脱原発シナリオが政府によって公開された． この中で，Fig. 1.1 のゼロ，15，20～25 の各シナリオに共通する点として火力発電が半分程 度のウェイトを占めていることに注目するべきである． 本来，京都議定書の「二酸化炭素排出量を 2012 年までに 1990 年比で平均マイナス 6%」 という目標を達成するには，火力発電の割合を 35%程度に減少させる必要がある．しかし ながら，脱原発によるエネルギー不足を再生可能エネルギーで補おうとしても，既存火力 発電無くしてはたちまち電力供給が滞ってしまう．また，再生可能エネルギーの普及が進 まない場合は Fig. 1.2 のようにほぼ火力に偏った発電構成になってしまう．これでは，明ら かに京都議定書目標の達成が困難である．また，2012 年以降に制定される予定のポスト京 都議定書における削減目標を達成するのも困難が予想される． これらの課題に対し，現在は再生可能エネルギーを筆頭とする新エネルギーの開発が期 待されている．それらは風力，太陽光，波力などの代替エネルギーを活用した発電方法で あり，核融合発電もその中のひとつに当たる．そのための研究として，大型核融合実験炉 ITER の建設が始まっている． Fig. 1.2 日本の 2011 年度電力発電実績割合．（2012 年度エネルギー白書 [3] より．）

3 / 83

1.2. 核融合発電

核融合発電とは，核融合炉と呼ばれる装置内で核融合反応をおこし，そのときに発生す る膨大な熱エネルギーを利用し，既存火力発電と同様にタービンを回して発電する方法（Fig 1.3 参照）である．現在，最も実現しやすい核融合反応のひとつとして D-T 反応（Fig. 1.4 参照）があるが，この方法では重水素や三重水素が消費され，燃料粒子不足となって長時 間の核融合反応を維持できない．そのため，長時間運転するには外部から粒子を供給する 必要がある． Fig. 1.3 核融合発電の模式図． （自然科学研究機構，核融合科学研究所のパンフレット[4]_より．） Fig. 1.4 D-T 反応の模式図．

4 / 83 従来の燃料供給法として，ガスパフ法と呼ばれる，プラズマに直接中性粒子ガスを吹き つける方法や，高速ペレット入射という氷のかけらを高速入射し，プラズマの持つ温度に より溶解，電離を起こし粒子を供給する方法，CT (Compact Torus) 入射[5] と呼ばれる大きな プラズマに小さなプラズマを，プラズマガンによって 100 [km/s]を超える高速で入射し粒子 供給する方法などがある．ITER ではプラズマ温度が高くなるので，ガスパフ法，高速ペレ ット入射ではプラズマ周辺での電離が多くなり，プラズマ中心での粒子供給が望めない（高 速ペレット入射はある程度の速度を持つことで中心での粒子供給は可能となる）． また，Fig 1.5 のように ITER は大型コイルを使用した磁気閉じ込め方式を採用している． 磁気閉じ込め方式とは真空容器の中に高温，高密度のプラズマを生成し，磁力線によるカ ゴで閉じ込めることにより核融合反応を起こす方式である．CT 入射では対象の大型プラズ マの持つ漏洩磁場による減速や，閉じ込め磁場により CT が影響を受け中心まで届かない[6] など，中心部へ供給するための解決すべき問題が残されている． Fig. 1.5 国際熱核融合実験炉（ITER）模式図 [7] ．

5 / 83

1.3. Compact Torus プラズマ

CT 入射プラズマの CT とは Compact Torus の略であり，磁気閉じ込め方式の一つである． CT に分類されるものとして FRC，Spheromak があるが大きな特徴は容器の中心にコイルが 存在しないことから移送ができるということである．Fig. 1.6, 1.7 に FRC と Spheromak の磁 場構造を示す． なお，本研究では両プラズマについて取り扱う． Fig. 1.6 FRC の磁場構造 [8] ． Fig. 1.7 Spheromak の磁場構造 [9] ． Spheromak はトカマク型に類似するプラズマであるが，トカマクに比べて，中心軸にコイ ル，ブランケット等を必要としない利点があり，より簡素な装置で生成することが可能で ある．また，中心軸に何もないということは，移送できるという特徴がある．すなわち， 同軸ガンなどの装置から射出することも可能である．

6 / 83 なお，トカマクと同様にトロイダル磁場を持つが，中心軸にコイルを有しない理由は， このトロイダル磁場を自らのプラズマ電流だけで作り，それをまた閉じ込めに使うからで ある．この磁場配位を長時間維持出来れば，外部トロイダル磁場に浪費されるエネルギー が不要になり，高効率な核融合装置となる．一例として，効率の目安となるベータ値は， トカマクの工学的ベータ値は 5 [%] 程度であるが，Spheromak は 40~50 [%] である．

しかし，Spheromak は生成する装置によっては，Fig. 1.8 に示すような tilting instability[9]

が 発生することが実験及び計算で確認されている．これは，プラズマ平衡達成後，崩壊過程 にてトロイダル磁場配位が傾いていく現象である．この現象はプラズマのアスペクト比（円 筒座標における R,Z 比）を小さく（1:1 に近く）することによって抑制できる．

→ →

7 / 83

1.4. 超高速フロー入射実験

本章の 1.2 で説明した CT 入射は磁場に影響されることにより，ITER の炉内プラズマ中 心部への粒子供給は見込めない．その問題点を改善するため，兵庫県立大学の永田等は CT の中性化入射[10] を提案した．これは，2006 年ロシアの Globus-M[11] における CT 入射実験に て，入射した CT が移送中に一部のイオンが再結合して中性粒子となって入射されているこ とが発見されたことに注目した結果である． CT の中性化入射では，CT 射出装置以降の移送管に水素ガスを詰めた中性化セルを設置 し，粒子の衝突による荷電交換反応を利用して CT を中性粒子フローに変換する．Figure 1.9, 1.10 に実験装置の模式図を示す．この方法は，ガスパフ法の磁場に依存しない点と CT 入射 の 100 [km/s]を超える入射速度を併せ持つことで ITER への粒子供給が期待できる． Fig. 1.9 CT 中性化入射（超高速フロー入射）の模式図． Fig. 1.10 CT 中性化セルの実験装置模式図 [12] ． （兵庫県立大 FACT 装置．）

8 / 83

1.5. 研究目的

本章の 1.4 で説明した CT 中性化において，優秀な点が多い一方，課題も山積みである． 特に，中性化過程による速度低下や挙動変化，また，中性化した後のガスフローの状態が 問題となる．中性化前の CT 内部の粒子はイオンであるため，磁場によって閉じ込められ形 を維持しているが，中性化後は利点である磁場の影響を受けない中性ガスフローであるこ とが仇となり，膨張していく可能性がある．また，膨張した場合，粒子がバラバラになり ながら管壁に衝突し，減速する可能性もあり，管壁への影響も予想される． 本研究では，これらの懸念される事項の内，中性化過程による速度低下や挙動変化を検 討・評価するため，磁場閉じ込め装置である Flux Conserver 内（Fig. 1.11 参照）での Spheromak プラズマを初期状態とし，射出後の中性化セル通過時のプラズマ時間推移を得る，高速並 列演算器を活用したハイブリッド・シミュレーションを構築した．また，もう一つの CT プ ラズマである FRC を用いた場合のハイブリッド・シミュレーションも構築した．

9 / 83

第 2 章 シミュレーション・モデル

2.1. 平衡計算

プラズマの平衡状態とはプラズマの巨視量が時間と共に変化しない状態のことである． 本研究では，プラズマの初期状態は平衡状態にあるという仮定からスタートする．ここで は，平衡状態にある CT プラズマを扱う際に用いられる方程式並びに CT パラメータを説明 する.

2.1.1. CT パラメータ

シミュレーションにおける CT プラズマ（Spheromak；FRC については後述する）の装置 及びプラズマ・パラメータは以下の Table 2.1 のように設定した[8]_．

Table 2.1 CT パラメータ．

電子温度

10 [eV]

CT 半径

5.0 [cm]

CT 長

7.5 [cm]

偏長度

1.5 [ - ]

CT 磁場

0.2 [ T ]

CT 入射速度

200 [km/s]

ここでの偏長度は，CT 長と CT 半径の比である．Figure 2.1 は CT パラメータをわかりや すく図で示したものである．Figure 2.1 はプラズマの R, Z 平面における二次元断面図（ポロ イダル面）であり，引き伸ばしたドーナツ（もしくはバームクーヘン）の断面のようなも のである． プラズマのポロイダル方向，トロイダル方向は Fig. 2.2 に示す．

10 / 83 Fig. 2.1 CT パラメータ図説 [8] ． Fig. 2.2 トロイダル方向とポロイダル方向 [14] ． なお，Fig. 2.1 の電子密度に関しては次節で説明する Grad-Shafranov 方程式によって求め る平衡解のパラメータに影響されるため，この段階では決定することができない．（ここで の CT パラメータは，予め想定される実験装置などによって決定されるものであり，他のパ ラメータに影響されない．）

11 / 83

2.1.2. Grad-Shafranov 方程式

プラズマの平衡は一般的に Grad-Shafranov 方程式を解くことによって求められる． Grad-Shafranov 方程式の導出過程については卒業論文[15]にある． 円筒座標系における Grad-Shafranov 方程式は， 2 2 0 2

1

z

r

r

r

r

d

dF

F

d

dp

r













































(2.1) であり，圧力分布とトロイダル磁場関数（ポロイダル電流関数）を与えることによって磁 束関数を求めることができる．そして，これを有限差分で解くことにより Fig. 2.3 のような 配位を作ることができる． Fig. 2.3 R, Z 平面における磁束分布図 [8] ． なお，圧力関数とポロイダル電流関数については以下のように設定した． ・圧力関数 0 lim 0 max

)

(















p

p

(2.2)

01

.

0

0 lim







(2.3)

12 / 83 ・ポロイダル電流関数















































0 0 0 0 max

0

0

















e

F

F

　　　

　　　

(2.4) なお，圧力関数とポロイダル電流関数内の



₀や



_eは Fig. 2.4 で示されている磁束面での 磁束値であり，



₀は最大磁束，



_eはセパラトリクスでの磁束を表す． Fig. 2.4 Spheromak の R, Z 平面（磁束，磁場）の図 [16] ． ここでの

p

や

F

は規格化（無次元化）された物理量である．プラズマを扱う式は，質量 m や電荷 e などを含んでいる．これらは 10-27 [kg] や 10-19 [C]のオーダであり数値計算する 際にはこれらを入力すると数値誤差を生じたり，無視小として勝手に取り扱われたりされ る恐れがある．また，次元を含むより無次元化したほうが都合良いことが多い． 以下に規格化用の定数を示す． max max max max wall

,

,

,

,

F

F

F

p

p

p

z

z

z

r

r

r

















物理量の上に添えた「－」は規格化したことを示す．

13 / 83 規格化及び無次元化を行った Grad-Shafranov 方程式は， 2 2 2 2

~

1

1

~

z

E

r

r

r

r

d

F

d

F

q

d

p

d

r

B













































(2.5) となり，それぞれの変数は， wall max 2 0 2 max 2 2 0 max 4 0

,

~

,

~

r

z

E

F

a

q

p

a















(2.6) で表される． (2.6)式の

B

~

はベータ値（1.3 節参照）の次元を持った変数であり，任意値を与えることに よってプラズマの形状を変化されることができる．また，

q

~

は安全係数と呼ばれ，トロイ ダル磁場の強さ等に影響する．E は偏長度（エロンゲーション）である．

2.1.3. プラズマ・パラメータ

2.1.1 節で与えた CT パラメータ及び前節(2.6)式より，平衡解算出後の各プラズマ・パラメ ータを求めることができる． ・最大磁束 プラズマの最大磁束



_max は CT 磁場の最大値 0.2 [T]より，平衡計算で得られた Z 方向磁 場の最大値を

B

_Zmaxとすると， 2 wall max max

2

.

0

r

B

_Z





(2.7) となり，本計算では

5

.

0



10

5

[Wb]

となった．

r

_wallは CT 半径（装置半径）である．

14 / 83 ・最大圧力 最大磁束と

B

~

より最大圧力

p

_maxは， 4 wall 0 2 max max

~

r

p









(2.8) と表すことができ，本計算では

31

.

82

[

Pa]

となった． ・最大電子密度 最大電子密度

n

_max は最大圧力より，

)

(

_{i e} max max max

T

T

k

p

kT

p

n







(2.9) と表すことができ，本計算では

9

.

93



10

18

[

m

-3

]

となった． （Ti及び Teはそれぞれ，イオン温度と電子温度 [eV]であり，同値である．） ・規格化時間 本シミュレーションでは，イオンの 1 ラーモア回転に必要な時間をもとに規格化を行っ ている．その周期は， max 2 wall





r

e

m

_i



(2.10) と表すことができ，本計算では

1

.

04



10

6

[s]

となった．miはイオン質量である．

15 / 83

2.2. ハイブリッド・シミュレーション

本研究では，水素イオンの軌道計算を行ない，PIC (Particle In Cell) 法を用いて超粒子の 情報を電子流体の情報に変換し，電子流体情報を場の計算に反映させ，その場の中でまた 粒子を動かすハイブリッド・シミュレーションを行っている．ここでは，その際に用いる 計算モデルやハイブリッド・シミュレーション法について説明する． 本節の本論は卒業論文[15] に記載してあるため，ここでは変更点について詳しく記述する．

2.2.1. イオンの初期分布

プラズマ内イオンを粒子的に扱い，マクロなプラズマを解析するためには多数の粒子を 同時に計算させる必要がある．正確な解析のためには全プラズマ粒子の軌道計算を行う方 法が最も良いが，現実の実験室のプラズマでは 1 [cm3 ]あたり 1010～1014_{個の粒子が存在する} ため，それらをすべてシミュレーションすることは不可能である．そこで，粒子シミュレ ーションでは，取り扱う粒子は実際の粒子の電荷質量比を一定に保ったまま，多数の粒子 の電荷と質量をまとめた超粒子として扱う．本研究では，水素イオンを扱っており，約 7 億個の粒子を配置し，軌道計算を行っている．（卒業研究時は約 1000 万個であったが，計 算機環境の改善より大幅に増強した．） Spheromak，FRC 等の CT プラズマは軸対称であるため，θ 方向を配慮しない R, Z の 2 次 元平面を計算することで，解析できる．そこで，Fig. 2.5 に示すように粒子を計算領域内に 配置し，軌道計算を行う．R 方向を 1024 分割，Z 方向を 2048 分割し，それぞれのグリッド 点間に様々な速度ベクトルを持つ粒子を配置する．速度については Fig. 2.6 のような Maxwell 分布を仮定し，それを 7 点分割して代表的な速度を取り出し，乱数で幅を持たせて 与えている．また，速度は R, θ, Z の 3 方向に関して与える必要がある．つまり，Maxwell 分布に沿った分布を作るために，1 グリット点ごとに 3

7

通りの速度を持った粒子を配置す ることになる．これらに分布関数，

















_

_

















i 2 2 2 2 3 i i i

2

)

v

v

(v

exp

2

T

m

T

m

n

f

r  z



(2.11) を用いて求めた重みを持たせる．ここで，n_i,m_i,T_iはそれぞれイオンの密度，質量，温度（エ ネルギー [eV]）を示す．

16 / 83 密度と分布関数の関係は，





f

d

v

n

(2.12) であるから，分布関数 f に微小速度dvをかけたものは，密度の次元をもっている． したがって，超粒子の重みを



とすると，

V

f



d



d



v



(2.13) となる． Fig. 2.5 格子点間粒子配置の模式図． ※卒業論文時のモデルでは，計算環境の制約上，R 方向を 129 分割，Z 方向を 257 分割し， なおかつグリッド間ではなく，グリッド点上に粒子を配置して計算を行った．（計算精度上 は粒子が多いほど良い．）また，速度分布に乱数的操作を導入していなかったため，初期速 度が同じ粒子が多数存在していた．

17 / 83 Fig. 2.6 Maxwell 分布図 [17] ．

2.2.2. 粒子軌道計算

粒子の軌道計算には運動方程式，











_









        

v

B

E



v

u

v

m

q

t

m

d

d

(2.14) を用いている．ここで，



はテスト粒子で，



はバックグランド粒子である．



_はテス ト粒子とバックグランド粒子との衝突周波数[18]_であり，

 



_





















        

_





_{2 2 3} 0 2 2

v

1

4

x

m

m

m

q

q

n

(2.15)

 

   





T

m

x

t

t

e

x

x t

2

v

,

d

2

2 0







(2.16) で与えられる．本研究のシミュレーションでは，テスト粒子は水素イオンを用い，バック グランド粒子は電子を用いている． また，(2.14)式では軌道を計算するのに磁場 B ，電場 E ，バックグランド粒子の流速uが 必要である．これらをどのように求めるかは，次節で説明する．

18 / 83

2.2.3. 磁場と電場の時間発展

ここで，Grad-Shafranov 方程式(2.5)より，プラズマの初期磁束が計算されているとする． 初期磁場は各方向成分別に，

r

r

B

r

F

B

z

r

B

_{r z}





















1



,

_

,

1



(2.17) と，求めることができる． 電場に関しては，電子の運動方程式より， e ei e e

1

en

p

en

R

E









(2.18) を使用する．ここで，

R

_eiは次節で説明するクーロン衝突による衝突項である． 衝突項の計算に必要な電子流速

u

_eは，後述する PIC 法により密度

n

_iと流束

γ

を集計し， i i

n

γ

u



(2.19) により，イオン流速

u

_iを求めた後，電流の定義式，



i e



e

u

u

j



en



(2.20) を用いて得られる．(2.18)式の右辺第二項は同様に PIC 法を使用して集計したイオン密度（＝ 電子密度）を用いて，圧力

p



nT

の関係から分布を求めて空間微分操作より算出する． さらに，ファラデーの法則，

E

B

_

_

_







t

(2.21) をオイラー法を用いて時間積分することにより，次タイム・ステップの磁束を求めること ができる．

19 / 83

2.2.4. クーロン衝突

プラズマ中の粒子は，バックグランド粒子とのクーロン衝突[18] により速度や軌道を変化 させる．摩擦による速度の減速現象を Slowing-Down Collision，軌道の変化をピッチ角散乱 という．ここでは、それぞれの現象について説明する．

2.2.4.1. Slowing-Down Collision

本シミュレーションにおける Slowing-Down Collision は，軌道計算対象のイオンが電子に 衝突し，運動エネルギーが減衰していく過程を再現している．詳しい導出過程は卒業論文[15] 参照． Slowing-Down Collision の衝突周波数



は，

 

 





















































S B S S S B S S S S S S S S

T

k

V

m

x

x

x

T

k

m

m

m

m

e

e

n

x

V

m

m

m

e

e

n

2

,

2

1

4

1

1

4

2 2 3 2 3 2 2 0 2 2 3 2 2 0 2 2 sl















(2.22) である．この



_slは，温度

T

_Sの Maxwell 分布したプラズマ中に存在するテスト粒子（イオン） の速度の e-folding time の逆数で，slowing-down collision frequency と呼ぶ．

また，(2.22)式で与えられる slowing-down collision frequency は

V

2



2

k

_B

T

_S

m

_S

の

場合は

 

x



1



であり，反対に

V

2



2

k

_B

T

_S

m

_S

で

は

 

32

3

4

x

x







である．よって，(2.22) 式は 2 つの極限を持ち，





























































































 S S S S S S S S S S S S S S

m

T

V

T

m

m

m

m

e

e

n

mV

m

T

V

m

m

m

e

e

n

2

2

1

3

2

1

,

2

1

2

8

2 2 3 2 2 0 2 2 2 2 2 3 2 0 2 2 sl



















(2.23) となる．

20 / 83 すなわち，プラズマの熱速度(

v

_th



2

T /

_S

m

_S )より充分に大きな速度をもつテスト粒子 の



_sl

は

，自身の速度の 3 乗に逆比例する．反対に熱速度より充分に小さいテスト粒子は熱 速度の 3 乗に逆比例した一定値に近づく．

2.2.4.2. ピッチ角散乱

ピッチ角散乱とは，荷電粒子同士（イオン同士）の衝突による散乱であり，そのピッチ 角は磁場に対して判定を行う．粒子の速度ベクトルを

v

，磁場を

B

とし，



₀を



1

cos

1



,

cos

_{p p} 0



















B

v

B

v

(2.24) と定義する．ピッチ角



_pの変化は













2 1 p 2 o p 0 n



1





1

















(2.25) で表す．ここで，



_pはピッチ角散乱するときの衝突周波数，



は擬似乱数を発生させる時 間刻みである．



の添え字の n は新しい速度方向を，0 は古い速度方向を示す．



はモンテ カルロ法を用いて決定する．

5

.

0

0



R



のとき「－」．

1

5

.

0



R



のとき「＋」． ここで、

R

は 0 から 1 の一様乱数である．ピッチ角散乱を考慮した粒子軌道計算をする には，この磁場に沿った座標系から円筒座標系に変換させる必要がある．次に円筒座標系 （磁場の向きを極にとった場合）に変換した各方向成分を示す．なお，旋回移送



_gは Fig. 2.7 のように，0~2



の間で乱数的に一様に与えている．なお，

_v

は絶対速度である．

21 / 83 Fig. 2.7 旋回移送角の図 [17] ．























































































g p g p p

sin

sin

cos

sin

cos

cos

0

sin

0

1

0

sin

0

cos

1

0

0

0

cos

sin

0

sin

cos





























v

v

v

v

v

v

r z (2.26)





















































sin

cos

sin

sin

sin

cos

sin

cos

cos

cos

cos

sin

sin

sin

cos

sin

sin

sin

sin

cos

cos

sin

cos

sin

cos

g p g p p g p p g p g p p

v

v

v

v

v

v

v

v

v

v

v

















z r (2.27) 2 2

sin

z r r

B

B

B







(2.28) 2 2

cos

z r z

B

B

B







(2.29) 2 2 2

sin

z r

B

B

B

B







 



(2.30) 2 2 2 2 2

cos

z r z r

B

B

B

B

B













(2.31)

22 / 83 (2.27) 式を初期値として扱うことにより，ピッチ角散乱を考慮した粒子軌道計算を行な うことができる．また，ピッチ角散乱の場合の衝突周波数



_pは，

 

 

2 z 2 2 r

v

v

v

































 



    













v

,

2

v

,

1

2

1

1

v

1

8

2 2 3 2 2 0 2 2 p

T

k

m

x

dx

x

d

x

x

m

e

e

n

B β β α α β α β (2.32) であり，衝突相手の粒子のパラメータとして

n

_βは密度，

e

_βは電荷，



_β

 

x

はエラー関数を 表し，衝突元の粒子のパラメータとして，

e

_は電荷，

m

_は質量，

v

_は相対速度を表す． また，(2.32)式で与えられる衝突周波数は前節の Slowing-Down 周波数と同様に，衝突元 の粒子の速さに対する極限を有する．以下に，特に重要な

_













  

m

T

k

_B

2

v

2 の場合を記 す．   











T

k

v

m

x

x

m

e

e

n

B α α β α β

2

,

1

3

4

v

1

8

2 2 3 2 2 0 2 2 p

















(2.33) ※卒業論文時のモデルでは，トロイダル方向（θ 方向）磁場が存在しない FRC プラズマに 対して導出したものを使用しており，ピッチ角を計算する段階で磁場ベクトルに対する角 度判定が正しく行えていなかった．これにより，Spheromak プラズマに適用して計算した場 合，粒子が部分的に密集する高密度域や非常に低密度な領域が発生する現象が起こった．

23 / 83

2.2.5. PIC (Particle In Cell) 法

本研究では，計算領域をセルに分割し，格子点上の物理量を求めるシミュレーション技 法を採用している．セル内部に含まれる超粒子の位置や速度情報を格子点に反映させ，粒 子の情報を流体の情報に変換するのが目的である．ここでは，その手法の一つである PIC (Particle In Cell) 法について紹介する．本研究では円筒座標系で PIC 法を用いており，Fig. 2.8 のようにセルで空間を区切り，各セル内にある超粒子の重み（情報）を格子点に変換する． PIC 方法自体は卒業論文[15]参照してほしい． Fig. 2.8 円筒座標系における PIC の概念図 [17] ．

2.2.5.1 勾配補正フィルタ

卒業論文[15]_{段階にて解決していなかった PIC 集計における一番の問題点として，境界セ} ルの取り扱いが挙げられる．Figure 2.9 に一例として，FRC プラズマの理論的に導出した平 衡密度分布とシミュレーションで PIC 集計を用いて得た初期密度分布（平衡状態に相当） を示す． なお，Fig. 2.9 は中央面分布である．中央面分布とは，プラズマの中心部（R, Z 二次元平 面の Z = 0）を R 軸に沿って各種物理量を抽出した分布である．

24 / 83 Fig. 2.9 理論値及び PIC 集計による中央面密度の比較図． Figure 2.9 にて，R = 1.0 境界は密度がほぼ 0 のため目立たないが，装置軸（R = 0.0）にて 理論値と集計値に大きな差が発生している．これは，PIC の境界セルの問題であり，Fig. 2.8 における

r

j_1

,

r

j

,

r

j_1 の関係で説明できる． 一次元的に考えると境界以外では，とあるメッシュ点

 

r

_j の値はその両隣の点

 

r

_j1 と

 

r

j1 によって囲まれた二つのセル内粒子の重みで決定される．しかし，R = 1.0 境界では装 置壁，R = 0.0 では装置軸をそれぞれ仮定しているため，粒子による重み情報は一つのセル 分となってしまう．よって，特に高密度かつセルが極めて小さい軸境界値が大きく影響を 受ける． この問題に対して，本研究では 2 つの対処を試行した．まず，一つ目は平滑化である． ・平滑化の場合 一次元 7 点平滑化式(2.34)[19]_{をすべての境界に対して適用し，境界値の低下抑えようとし} た．



3

6

2

15

1

20

15

1

6

2 3



64

1

     















_{i i i i i i i} i

f

f

f

f

f

f

f

f

(2.34)

25 / 83 ・勾配補正の場合 平滑化の場合と同様に，一次元分布上で求めたい境界の値を

f

_iとし，その隣接値を 2 1

,

  i i

f

f

とする．それぞれを 2 階微分の項までテイラー展開すると， 2 2 2 1

2

1

x

x

f

x

x

f

f

f

_{i i}



















 (2.35)









2 2 2 2 2 2 2

2

2

2

2

1

2

x

x

f

x

x

f

f

x

x

f

x

x

f

f

f

i i i





































 (2.36) となり，fi₂ 4fi₁より，

x

f

f

f

x

f

i i i















 

2

4

3

_{1 2} (2.37) が求まる． 境界値と隣接値の間には勾配がないと仮定すると，この(2.37)式は 0 になる．よって境界 値は，

3

4

_1



_2



i i i

f

f

f

(2.38) となる．この(2.38)式をすべての境界に適用して勾配補正フィルタとした．なお，装置軸と 装置壁では符号が反転するのみである． ※なお，卒業論文時のモデルでは，境界値を隣セルで強制的に上書きしていたため勾配無 しではあるが，極端な分布になってしまい，数値ノイズが発生していた．

26 / 83 ・各フィルタの効果

Fig. 2.10 各フィルタ適用時の中央面密度比較図．

Figure 2.10 に各フィルタ適用時の中央面密度比較図を記す．Figure 2.10 の赤線は Fig. 2.9 と同様に一切の補正を行っていない生の PIC 集計結果を表している．青線は境界値のみに 一次元 7 点平滑化式を適用した場合であり，緑線はさらに隣値に適用した場合を示す．最 後に，黒線は境界を勾配補正フィルタによって求めたものである． まず，境界値のみを平滑化した場合は，平滑化前よりは改善しているが，隣接値より勾 配が発生しており，好ましくない．次に，この勾配を低減するために隣値も平滑化したの が緑線であり，勾配そのものは減少しているが，依然として目立ってしまう．最後に，境 界を勾配補正フィルタによって求めたものは，Fig. 2.9 の理論値とほぼ同じであり，好まし い結果となった．（Figure 2.10 の各線が重なって判別しづらくなるため，理論値線は省略し た．） よって，本シミュレーションでは勾配補正を採用した．

27 / 83

2.2.6. 平滑化

本研究では，ハイブリッド・シミュレーションを行っているため，数値振動が発生しや すい．よって，各物理量を求めた後で平滑化し，高周波成分を除去している．詳しい説明 は卒業論文[15] に記載してあるため，ここでは，B-spline 2 次元平滑化式（節点追加方式）に よる平滑化について紹介する程度に留める． Fig. 2.11 残差二乗和が最小となる spline 関数 [17] ． X, Y 平面上の格子点



x

_i

,

y

_j



;

i



1

,

2

,



,

n

_x

,

j



1

,

2

,



,

n

_y，観測値

f

_i,_j



f



x

_i

,

y

_j



， 観測誤差



_i,_j





x

_i





y

_j，残差二乗和の許容値 2 t



，及び X, Y 方向の初期節点列 s n







1

,

2

,



,

;



1

,



2

,



,



l_s，通常は

n

s



n

l



2

， 1

min

 

i

，

i

x





 

i i n_s



max

x



，

 

，

min

1 j j

y





 

j j l_s



max

y



が与えられたとき，節点を適応的に追加していき，残差二乗 和が与えられた残差二乗和の許容値 2 t



以下となるような双 M 次の B-spline 平滑化式を求め る． 節点列を更新しながら次第に節点の個数を増加させ，残差二乗和が許容値 2 t  以下になる まで反復する．これが，Fig. 2.11 のような平滑化結果となる．

28 / 83

2.3. 中性化

本節では，CT が中性化セル内を通過する際に，中性化セル内中性ガスと荷電交換して中 性化していく過程のモデルを説明する．また，荷電交換によって中性流体がイオン化し， 発生する「低速イオン」についても解説する．

2.3.1. 軸方向 NBI モデル

初期モデルでは，実際の実験装置と同様に中性粒子を分布させた計算領域中に CT 入射 （移送）し，CT 内のイオンを軌道計算しつつ，中性粒子との荷電交換反応をモンテカルロ 法により再現することで中性化率を計算しようと試みた．しかし，このモデルでは CT 内の イオンに移送速度を追加させて計算するため，非常に長い計算領域（CT 中性化セルの長さ に比例）が必要であり，それに伴う計算負荷の増加に計算環境が追い付かないという問題 が発生した． それを解決するために，Fig. 2.12 のようにシミュレーションではセル内の中性粒子自体に 軸方向の並進成分速度を持たせる NBI（Neutral Beam Injection）モデルを採用した．このシ ミュレーション・モデルと実験モデルは粒子間の相互作用において等価であるため，実験 モデルと同じような結果が得られる．イメージとしては，静止している CT に対して中性化 セルが入射してくるのを想像すれば良い．

29 / 83

2.3.2. 中性流体

本モデルでは，CT 中性化セル内の中性粒子を流体として取り扱う．これは，個々の粒子 として軌道計算する場合より計算負荷が軽く，電離していないため粒子としての軌道を詳 しく分析する必要がないからである． Fig. 2.13 静止 CT に対する中性流体入射のモデル図． Figure 2.13 のように，計算領域内で静止している CT に対して一定のセル幅で区切られて いる中性流体が入射速度に合わせて入射する． セル幅は，入射速度を v [m/s]と衝突判定時間



t

[

s

]

を乗算し，さらに衝突判定回数 N を 使用して調整している．例として，

]

m/s

[

10

2

]

km/s

[

200

5







v

]

s

[

10

04

.

1

10

04

.

1

01

.

0

01

.

0

8 6  















t



N = 1 とすると，セル幅は

2

.

08



10

3

[

m

]

となる． なお，衝突判定の間隔は磁場・電場の分布更新と同じ間隔で行うため，

0

.

01



とした．ま た，中性化流体のセル数とメモリの空き状況より衝突判定回数を 1 とした．（メモリが不足 する場合は，判定回数を増やしてセル数を減らす処理を行う．）

30 / 83 中性流体の密度に関しては，中性化セル内のガス圧から算出しており，1 [mTorr]と仮定す ると，

]

[m

10

32

.

3

[mTorr]

1

,

[J/K]

10

380658

.

1

[K]

300

[Pa]

10

133.322

]

[m

3 19 23 -3 3   

















p

p

n

n となる．

2.3.3. 中性化判定

静止している CT 内イオンは入射してきた中性流体の各メッシュ間に突入する．この際， 中性化は CT イオンの超粒子的重みを減らすことによって再現している． とあるイオン超粒子（i 番目）の重み（粒子数に相当）を Wiとすると，

 

_{ }

v

n

t

W

t

t

W

n i i







d

d

(2.39)



t

t



W

 

t



n

v

t



W

_i







_i

1



_n





(2.40) のように時間変化を解くことができる．ここで，nnは中性流体の密度であり，



v

は反応 速度係数である．（反応速度係数は次節で説明する．）なお，前節と同じく，衝突判定の間 隔は磁場・電場の分布更新と同じ間隔で行うため，

0

.

01



とした． 同様に，荷電交換によって中性流体密度が低下していくのを，

v



i n n

d

d

n

n

t

n





(2.41)



n

t

t



t

n

t

t

n

_n

(





)



_n

(

)

1



_i

(

)



v



(2.42)

31 / 83 で再現している．ここで，niは CT イオンを PIC 集計した密度であり，



v

は(2.39)式と同 じ反応速度係数である．

2.3.4. 反応速度係数

本研究では，水素の荷電交換に注目して中性化を再現している．物理的には，CT の温度 低下による再結合反応での中性化が考えられるが，それは先行研究[8] によって支配的ではな いことが判明しているため，今回はモデルに考慮しない．また，同様に Fig. 2.14 のような 水素への電子衝突による電離や励起も荷電交換に比べて支配的な反応ではないため，考慮 していない． Fig. 2.14 水素の各種反応式 [8] ． 荷電交換反応の反応速度定数を算出するのにあたって，本モデルでは衝突断面積より簡 易的に導出している．

32 / 83 Fig. 2.15 NIFS データベースの荷電交換断面積データ [21] ． Figure 2.15 は自然科学研究機構・核融合科学研究所（NIFS）の原子分子数値データベー スから取得した水素の荷電交換断面積データである．まず，この図より，プラズマの温度 が 10 [eV]程度なら荷電交換断面積



が 19 10 [m2] (1015[cm2])であることが読み取れる． 次に，プラズマ内イオンと中性流体が衝突する際の相対速度 v は，CT 入射速度 200 [km/s] から 5 10 0 . 2  [m/s]と見積もることができ，これは 100 [eV]のエネルギーに相当する． 荷電交換断面積



と相対速度 v より反応速度係数



v

は，

/s]

[m

10

0

.

2

10

0

.

2

10

3 14 5 19  











v



となる．

33 / 83

2.3.5. 低速イオン

荷電交換によって CT 内の高速イオンは中性化し，ガス状態に変化する．これとは逆に， 中性化セル内の中性ガス（シミュレーションでは流体として取り扱い）はイオン化する． この際，中性化セル内の中性ガスは室温であるので，発生するイオンも CT 内高速イオンに 比べて低温であり，熱速度が相対的に低い．よって，CT 内高速イオンと区別するため，「低 速イオン」と本研究では称する． Fig. 2.16 低速イオン生成のモデル図． 低速イオンは，高速 CT イオンが中性流体と荷電交換した際，同量（同じ粒子数）発生す る．よって，本シミュレーションでは判定時間毎の CT イオン超粒子重み低下量を PIC 集計 し，空間分布に変換している．この空間分布に基づいて，新しく低速イオンを超粒子とし て Fig. 2.16 のように計算領域に均等配置し，軌道計算を行っている．低速イオンの配置メ ッシュは PIC メッシュと同じであり，均等配置のため，重みの調整によって発生量の空間 的大小を調整している．（重みの定義式は(2.13)式参照．） 低速イオンの効果としては電場への影響が期待される．低速イオンは高速イオンより速 度が低いため，2.2.4.で説明したクーロン衝突の衝突周波数が高く，それは電場の衝突項増 加につながる．電場の定義式は(2.18)式参照． また，軸方向 NBI モデルの関係以上，相対的に見ると低速イオンは計算領域を入射速度 200 [km/s]の速さで突破していることになる．（一部は CT の磁場に補足されるが，大部分は 中性化セル内に残ると考えられるため．）よって，背景の電子流体と連動して電場が形成さ れ，CT の全体的な減速及び膨張という効果も期待される．

34 / 83

2.4. 高速並列演算環境

本研究では，水素イオンに相当する，重みを持った超粒子約 7 億個の粒子軌道計算を行 う．一般家庭で使用する PC 一台では，数マイクロ秒の実時間に相当するシミュレーション の計算に数か月かかるため，多数の PC を並列接続して演算する高速並列演算環境を新たに 構築した． 一般的に，並列計算機は 2 つの種類，密結合クラスターと疎結合クラスターに分類され る．この 2 つの大きな違いはそのメモリ・モデルにある．複数のプロセッサが共有される 大容量メモリ・バスに接続された共有メモリ・モデルと，プロセッサ毎にメモリが接続さ れている分散メモリ・モデルである[22] ． スーパー・コンピュータと呼ばれる計算機は、主に共有メモリ・モデルであり、CPU を 数百から数千個接続した並列計算機である。逆に，本研究で用いる並列演算環境は分散メ モリ・モデルであり，複数の PC を LAN によって接続し，仮想的に一台の計算機として利 用する．なお，利用する PC は 4 台である．

2.4.1. ハードウェア

並列演算環境を構築するにあたって，当研究では市販の自作 PC パーツを使用することを 重視した．故障した時のメンテナンス性に加え，実際に構築する時の設定のし易さを配慮 すると，メーカー製の既製 PC よりも価格などを含めて利点が大きい． 使用したパーツの一覧を以下の Table 2.2 に示す．また，完成した並列環境を設置したラ ックを Fig. 2.17 に記す． Table 2.2 使用パーツ一覧． 製品名 メーカー 型番 数量

CPU AMD FX-8170 3.9 [GHz] 8 Core 4

Motherboard ASUS M5A88-M 4

PC Case Scythe GX-3901-BK 4

RAM Silicon Power PC10600 DDR3 16 [GB] (4[GB]x4) 4

HDD WD S-ATA 500 [GB] 7200 [RPM] 4

35 / 83

Fig. 2.17 並列演算環境の全体図． （新並列環境は奥のラック最上段 4 台である．）

2.4.2. ソフトウェア

並列演算環境の土台となる OS は，長時間稼働の安定性を配慮して，数ある Linux の中の RHEL（Red Hat Enterprise Linux）Clone である CentOS 5 を選択した．CentOS は有料である RHEL のコードを無料公開したものであり，大企業や公立研究機関での利用実績がある．な お，アメリカのフェルミ国立加速器研究所が Linux クラスターを構築する上で自主開発した Scientific Linux（旧 Fermi Linux）も RHEL Clone である．

次に，前節で述べたメモリ・モデルによって，要求されるプログラムやソース・コード が異なってくる．分散メモリ・タイプの並列プログラミングとして例を挙げるなら，PC 間 の通信を制御する，MPICH や OpenMPI 等の MPI (Message Passing Interface) がある．本研究 では，メモリ使用量の軽減や互換性の関係上，MPI-2 規格に対応した OpenMPI を使用する． なお，MPI とは MPI フォーラムで規格化された並列コンピューティングを利用するため の標準化された規格である．この MPI は特定のソフトウェアがあるわけでなく，並列処理 用のメッセージ・パッシング規格[23, 24]_{である．並列コンピューティングにおいて数値計算} の過程に各ノード間を様々なデータが送受信されるわけだが，その送受信の手順を定めた ものが MPI である．