逆磁場テータピンチ生成における電子挙動シミュレーション

目次

第 1 章 序論

1

1-1 核融合 1

1-2 磁場反転配位プラズマの特徴 4

1-3 FRC プラズマの生成法 7

1-3-1 回転磁場法(Rotating Magnetic Field:RMF) 7

1-3-2 異極性スフェロマク合体法(Counter Helisity Spheromak Meraging) 8

1-3-3 逆磁場テータピンチ(FRTP)法 9 1-4 研究目的 10

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

11 2-1 シミュレーション条件及び各種パラメータ 12 2-2 電子の初期分布 15 2-3 粒子軌道計算 17 2-3-1 電子の運動方程式 17 2-3-2 場の計算法 17 2-3-3 slowing-down collision 20 2-3-2 場の計算法 21 2-4 FRC プラズマの平衡状態 24 2-4-1 MHD 平衡 24 2-4-2 理想 MHD 方程式 24 2-4-3 Grad-Shafranov 方程式 25 2-5 PIC 法 27 2-6 モンテカルロ法 31 2-7 衝突断面積 32 2-8 正規分布 32 2-9 ボックスミューラー法 34 2 -10 有限差分法 37 2-11 反復法（SOR 法） 38

第 3 章 電子挙動シミュレーション結果と考察

41 3-1 シミュレーション結果 41 3-1 反射の導入 52 3-3 反射導入後のシミュレーション結果 53 3-4 G-S 方程式から求めた FRC の平衡状態と無反射モデルとの比較 63 3-5 考察 66

第 4 章 まとめと展望

67 4-1 まとめ 67 4-2 今後の展望 67

参考文献

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・68

謝辞

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・70

第 1 章 序論

1-1 核融合 今日，世界の人口の増加や発展途上国の経済成長に伴って，エネルギー消費 量が急激に増加しており，化石燃料の枯渇が問題となっている．今世紀の後半 からエネルギー不足が深刻化すると予測されており，一刻も早く，化石燃料に 代わるエネルギー源の開発が求められている． Fig. 1-1 各資源の確認可採埋蔵量と利用可能な年数 人類究極のエネルギーと言われる核融合エネルギーは，資源の枯渇の心配が なく，大規模な発電ができ，地球環境への負荷も少ない．そのため，実現すれ ばエネルギー問題を解決することができる． 核融合とは軽い原子核同士が融合して，より重い原子核になる反応である． 原子核同士がある程度接近すると，原子核の間に働く引力（核力）が静電的な 反発力（クーロン力）に打ち勝って融合する．この際，衝突前の総質量よりも 衝突後の総質量の方が軽くなる．これを質量欠損という．この質量欠損分は， アインシュタインの導き出した式 2 m c E  から生成粒子の運動エネルギーとな る．このエネルギーを電気エネルギーとして利用するのが核融合エネルギーで ある．利用可能な核融合反応として考えられているものは以下の通りである． ) MeV 03 . 3 ( p ) MeV 01 . 1 ( T D D   ) MeV 06 . 14 ( n ) MeV 52 . 3 ( He T D  4  ) MeV 67 . 14 ( p ) MeV 67 . 3 ( He He D 3  4  ) MeV 45 . 2 ( n ) MeV 82 . 0 ( He D D  3 

4.8MeV He T n Li6   4 n He T 2.5MeV) n( Li7    4 ここで， D は重水素， T は三重水素， p は陽子， n は中性子を表す． 核融合反応を起こすには，クーロン力により反発しあう原子核同士を核力の 働く距離まで近づけなければならない．このときに必要な速度は１千㎞/s 以上で， 原子核を 1 億度以上に加熱することで得られる．また，発電するには核融合反 応を定常的に起こさなければならない．定常的に反応を起こすには高密度にす る必要がある．しかし，原子核は正の電荷を持っており，高密度にするために は閉じ込め領域全域において電荷中性になっていなければならない．原子核と 電子がバラバラになった状態，つまりプラズマ状態がこれに相当する．これら のことから，核融合反応を定常的に起こすには，高温・高密度のプラズマを作 り，反応が起こるまで閉じ込める必要がある． 核融合においてプラズマを閉じ込める方法は，大きく分けて 2 種類ある．1 つ は，燃料をつめた小球（ペレット）に四方から大出力のレーザーや粒子ビーム 等を照射して圧縮し，核融合反応を起こす慣性閉じ込め方式である．小球の外 層は氷状，中側はガス状になっている．この燃料小球に強力なレーザーを当て ると，表面は熱せられて外側に膨張する．その反作用で中のガスが圧縮され高 温プラズマになり，圧縮されたプラズマは，ほんの一瞬だけその場所にとどま ろうとする「慣性の力」で閉じ込められる．このようにして，中心部で核融合 が起こる． Fig. 1-2 慣性閉じ込め方式[1] もう 1 つは，磁場閉じ込め方式である．プラズマ中の電子はマイナス，原子 核はプラスの電荷を帯びており，電子や原子核は磁力線の周りを回転しながら 運動（ラーマ運動）する性質がある．そこで，コイルに電流を流した電磁石を 立体的にうまく配置し，磁力線のカゴを作り，その中に高密度のプラズマを閉

力線のカゴの形や電磁石の配置によってさまざまなタイプがある．その代表的 なものに「ミラー型」「ヘリカル型」「トカマク型」の三つのタイプがある．こ の中でもトカマク型は将来の核融合炉に最も有力とされている．2016 年頃の完 成を目指し，フランスのカダラッシュに建設することになった ITER（国際熱核 融合実験炉）でも採用されている． Fig. 1-3 磁場閉じ込め方式[1] Fig. 1-4 国際熱核融合炉（ITER）模式図[2] 本研究では，ミラー型に分類される磁場反転配位と呼ばれる閉じ込め方式を 使ったプラズマに関する研究である．

1-2 磁場反転配位プラズマの特徴 磁場反転配位（FRC : Field-Reversed Configuration）プラズマ[3]は，Fig. 1-5 に 示されるような磁場構造をしたプラズマである．FRC は開いたミラー磁場中に， 閉じた磁力線領域が形成されており，この二つの領域の境界を separatrix と呼ぶ． また，中心領域に磁場がゼロとなる場所（磁気中性点 : O-point）が存在し，そ の付近に高温のプラズマを閉じ込めることができる．開いた磁力線によっての みプラズマを閉じ込めるミラー方式では，開放端からの粒子損失率が大きく， またベータ値も低い．しかし，このような閉じた磁力線領域を持つ FRC は，磁 場閉じ込め方式の中で最も高いベータ値を持つ．ベータ値はプラズマ圧と磁気 圧との比を表し，一般的には   （プラズマ圧）／（外部磁気圧） (1．1) で定義される．直線的な磁力線に沿って磁場の変化がない場合， const. 2 ₀ 2    B p (1．2) となり，FRC の場合，磁気中性点が存在するので 0 0 2 0 2 ex 2 2 p B p B _{  }   (1．3) が成り立つ．ここで，B は外部磁場の大きさを表し，_ex p は磁気中性点での圧₀ 力を表す．つまり，磁気中性点でのベータ値は， 1 2 ₀ 2 ex 0     p B (1．4) となり，外部磁気圧はそのまま磁気中性点でのプラズマ圧力となる．ベータ 値が高い場合，ベータ値が低いものと比較すると，ある一定の外部磁場のとき に，より高い圧力のプラズマを閉じ込めることができる．つまり，FRC は他の 閉じ込め方式より弱い磁場でプラズマを閉じ込めることができるので，電力の 消費量が少なくて済み，経済的な閉じ込め方式であると言える． 高ベータプラズマである FRC プラズマは，アドバンスド核融合炉心プラズマ としての可能性を持っている．前節で触れた ITER で，実用化を検討されている 核融合反応は，重水素と三重水素による反応（D-T 反応）である．Fig. 1-5 に核 融合反応断面積の衝突エネルギー依存性を示す．D-T 反応は他の反応に比べ，反 応断面積が大きく，ピークが高い．そのため，他の反応より容易に核融合反応 が起きる．しかし，14．06MeV の高速中性子が発生するため，核融合炉の構造 材が放射化されてしまう．このため，核融合炉材料の耐久年数が短く，炉壁の 交換，放射化した材料の廃棄などに高いコストが必要とされる．それに対し， アドバンスド核融合反応として，重水素とヘリウム 3 による反応（ 3 He D 反応）

にくいが，高速中性子が発生せず，代わりに 14．67MeV の陽子が発生する．陽 子は電荷を持つため，磁場によって捕捉されるので，外部磁場を制御すること で核融合炉の損傷を防ぐことができる．また，磁場によって陽子をエネルギー 変換機に導き，陽子の持つエネルギーを電気エネルギーに変換することが可能 である． 3 He D 反応は炉壁の交換，放射性物質の廃棄にコストがかからず，陽 子から直接エネルギーを得られるため，経済的な反応である． このように，高ベータプラズマである FRC プラズマは，現在実用化が検討さ れているものよりも，経済性が高く魅力的である．しかし，プラズマの配位維 持時間が短く，FRC を実用化するためには配位維持時間を約 10 倍から 100 倍に 伸ばす必要がある．現在，FRC プラズマの配位維持時間が短い理由は明らかに されていない．しかし，可能性として （1）磁場揺動が大きくこれにより粒子が輸送されるため粒子閉じ込めが極 めて悪い． （2）異常抵抗によりプラズマ電流が減少し，閉じ込め磁場が減衰する． （3）セパラトリクス内部で回転不安定性が起こり，配位が崩壊している． などが考えられている．FRC プラズマの配位維持時間を伸ばすには，これら の可能性を検証し，その現象が起こるメカニズムを調べ，改善する必要がある． Fig. 1-5 FRC の磁場構造

O-point

開いた磁場 閉じた磁場

X-point

Separatrix

z



r

10-4 10-3 10-2 10-1 100 101 100 101 102 103  [ ba rn] E [keV] D-T D-D D-3He Fig. 1-6 核融合反応

1-3 FRC プラズマの生成方法

FRC プラズマには大きく三種類の生成法が存在する．

1-3-1 回転磁場法(Rotating Magnetic Field:RMF)

外部から与えた回転磁場がプラズマ中の電子に周方向駆動力を与えること で，プラズマ電流を発生・維持する方法が回転磁場法による生成法である[4-6]． Fig. 1-7 に示すようにイオンサイクロトロン周波数より高い周波数の回転磁場を プラズマ柱に印加することにより，電子のみ回転磁場に追従し，プラズマ中に 周方向電流を駆動するものである． Fig. 1-7 回転磁場によるプラズマ電流駆動の原理．

1-3-2 異極性スフェロマク合体法(Counter-Helicity Spheromak Meraging) 合体生成は Fig. 1-8 のように互いに逆向きのトロイダル磁場（磁気ヘリシテ ィ）を有するスフェロマクプラズマ 2 個を軸対称合体させて，トロイダル磁場 （磁気ヘリシティ）を打ち消し，熱エネルギーに変換することにより～１の FRC を形成する[7-9]．逆向きヘリシティのスフェロマク同士の磁力線は再結合 によってトロイダル方向に半周程度引き伸ばされた輪になり，それが解消・緩 和する間にトロイダル磁場エネルギーに変換される． （a）向きの違う 2 つのスフェロマク （b）合体した直後の FRC （c）合体後加熱されて成長した FRC Fig. 1-8 スフェロマク合体を用いた FRC 生成

spheromaks

=0

=0

1-3-3 逆磁場テータピンチ(FRTP)法 逆磁場テータピンチ法とは，最も典型的な FRC プラズマの生成方法で，数ア ルヴェン時間での生成で，θピンチ放電による電子電流の駆動と磁気再結合に よる配位形成（磁場反転）とプラズマ加熱（径方向，軸方向圧縮・衝撃波加熱） を同時に実現するものである[10-11]． Fig. 1-9 のような円筒形のテータピンチ装置内にガスを封入し予備電離プラズ マを生成し，バイアス磁場 B0を重畳する．その後，逆向き磁場を急速に印加し， 急速な磁場変化により発生する誘導電場 Eθ=(r/2)(dBz/dt)によりプラズマ電流 Iθ を駆動する．この電流と印加した磁場とのローレンツ力によって径方向の衝撃 波加熱と断熱圧縮加熱により第一段階のプラズマ加熱が起こる．同時に装置端 部で磁気再結合が起こり，閉じたポロイダル磁場が形成される．その後，閉じ た磁力線の張力による軸方向圧縮加熱およびポロイダル磁場のジュール加熱に より第二段階のプラズマ加熱が行われ，FRC プラズマが生成される．この生成 の過程は，数アルヴェン時間でおこる．このような生成時には，プラズマのダ イナミックな状態が続き自己組織化現象などの結果として FRC の生成とプラズ マ加熱が同時に進行する．

（1）真空引き，ガス詰め （2）予備電離 （3）バイアス磁場印加 (4)主圧縮磁場印加 （5）径方向圧縮 （6）軸方向圧縮 Fig. 1-9 FRTP 法による FRC 生成 1-4 研究目的 前節で，プラズマの生成方法について述べた．現行の FRC プラズマシミュレ ーションでは，FRC の平衡を初期条件として計算している事が多く，生成直後 のプラズマの状態を知ることができない．よって，より正確なシミュレーショ ンを行う上でコイル電流等実験系を反映し，電離過程などを考慮した CT 生成の 粒子シミュレーションにより生成過程や生成後の CT の特徴を明らかにする必 要性がある． 本研究の目的は FRC 生成の粒子シミュレーションを通し，電子の挙動から FRC プラズマ生成の過程と生成直後の FRC の状態を再現することである．

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

本研究では FRTP 法による FRC プラズマ生成の電子軌道シミュレーションを 行う．FRTP 法で用いる実験装置の例として日本大学の実験で使用されている NUCTE (Nihon University Compact Torus Experiment) -Ⅲを挙げる．Fig. 2-1 は NUCTE-Ⅲの装置図である[10-11]． FRTP 法では外部コイルにバイアス電流をかけておき，短い立ち上がり時間で ピーク電流まで電流値を変化させることにより多数の電離を起こし，生じるθ 方向の電流によってプラズマを閉じ込める．そこで外部コイルとプラズマ電流 をそれぞれ計算することで，実験系を反映したモデルとして電子シミュレーシ ョンを行う．ここでは，その際に用いる計算モデルやシミュレーション方法に ついて説明する．

2-1 シミュレーション条件及び各種パラメータ シミュレーションに使用するパラメータは Table 2-1 のとおりである． バイアス電流からピーク電流へと変化させる時間をシミュレーション時間と し，この間をシミュレーションする．このときのステップ数と外部コイルにか ける電流値の関係を Fig. 2-2 に示す．また外部コイルにかけたバイアス電流によ って生成される初期の磁束と磁場を Fig. 2-3 に示す． 軌道計算時間１ステップ は電子のサイクロトロン周波数で規格化し，1.14[nsec]程度として計算を行った． FRC は軸対称かつ反転対称であるため，装置中央面近傍において Z 方向に一 様かつθ 方向に対称である．これらを考慮し装置中央面 r 成分のみの 1-D モデル として計算することができる．なお装置壁までを計算領域とし，計算領域から 出た粒子は損失したものとして扱う． Table 2-1 シミュレーションに使用したパラメータ 装置半径 : 0.20 [m] 外部コイル半径 : 0.22 [m] バイアス電流 : 10 [kA] ピーク電流 : -500 [kA] シミュレーション時間 : 5 [sec] 初期電子温度・密度 : 10 [eV] ・ 1.0×1018 [m-3] 中性ガス密度 : 3.3×1020 [m-3] 初期粒子数 : 50，000

(a)

(b)

Fig. 2-3 初期磁束・磁場の径方向分布 (a)初期磁束 (b)初期磁場

シミュレーションの概要として Fig. 2-4 に計算コードのフローチャートを示す． Fig. 2-4 計算コードのフローチャート

外部コイルによる磁束計算

電子の軌道計算

・モンテカルロ法を用いて乱数による

電離衝突の判定を行う

磁場・電場を計算

プラズマ電流による磁束計算

PIC法を用いて電子密度n

_e

・流速u

_e

を集計

(PIC:Particle In Cell)

プラズマ電流を算出する

t+Δt

計算パラメータの入力

粒子の初期条件の設定

START

t <

simulation

time

STOP

2-2 電子の初期分布 プラズマイオンを粒子的に扱い，マクロなプラズマを解析するためには多数 の粒子を同時に計算させる必要がある．正確な解析のためには全プラズマ粒子 の軌道計算を行う方法が最も良いが，現実の実験室のプラズマでは 1cm3_あたり 1010～1014個の粒子が存在するため，それらをすべてシミュレーションすること は現実的ではない．そこで，今回行う粒子シミュレーションでは，取り扱う粒 子は実際の粒子の電荷質量比を一定に保ったまま，多数の粒子の電荷と質量を まとめた超粒子として扱う．本研究では，5 万個の粒子を超粒子として配置する． 粒子の初期配置は Fig. 2-5 に示すように計算領域内に一様乱数を用いて電子を 配置し，速度については Fig. 2-6 のようなガウス分布にそうような正規乱数によ って与えた．なおイオンは静止しているものと仮定し，密度のみ考慮していく． イオンの初期密度に関しては電子の初期分布から集計された密度と同一である ものとし，電離が発生し新たに電子が増えた位置で電子と同様に密度が増加す るものとする．

Fig. 2-5 粒子の初期配置

Fig. 2-6 ガウス分布

v

f

2-3 粒子軌道計算 2-3-1 電子の運動方程式 粒子の軌道計算には運動方程式



 





_

 i e ei e e e e e q m t m v v B E  v d d (2.1) を用いている．ここで，eは電子で，i はイオン粒子である．_eiは電子とイオ ン粒子との衝突周波数であり，

 



_                      _{2 2 3} 0 2 2 v 1 4 x m m m q q n (2.2)

 

x 2 xet tdt 0



       T m x 2 v2  (2.3) で与えられる． また，(2.1)式より軌道を計算するには磁場 B ，電場 E が必要である．これらを どのように求めるかは，次節で説明する． 2-3-2 場の計算法 場の計算方法について説明する．外部コイル・プラズマ電流それぞれからの 磁束・電場・磁場を計算する． まず，外部コイルに流れる電流より磁束を計算する．ファラデーの電磁誘 導の法則 t       E B (2.4) より，両辺を面積分すると





_



        E dS B dS t (2.5) となる．また，円筒断面積はdS2rdrであるため，磁束 の定義式は S B d 2 1 _ 

_

  (2.6) と与えられる．ビオサバールの法則から



_   C d 3 0 | | ) ( 4 r s s r s B   となるので，以上から外部コイルに流れる電流から求まる磁束_exは

       d rr r r r t I r c c c ex



   2 0 2 2 0 cos 2 cos 4 ) ( (2.7) ここでr は求めたい地点から外部コイルまでの距離，r は磁束を求めたい位置で_c ある．もとめた磁束_exから，磁場B_zexは r r B ex zex _  1  さらにファラデーの法則から(2.5)式の左辺にストークスの定理を適用し，右辺 に(2.6)式を代入すると，(2.5)式は







         rE_ t d 2 2 E l (2.8)   rE t    (2.9) 以上から電場E_exは t r E ex ex       1 (2.10) のように求められる． またプラズマ電流によって作り出される磁束_pを計算する．PIC 法を用いて 電子密度 ne，流速 ueを集計し電流 j を求めると    eneue j  (2.11) となる．ここで集計された電流 j を用い，アンペールの式より _ B j 0  (2.12) 磁束を計算すると ) ( ) ( 1 0j r r r r r                (2.13) となるので後述する有限差分法を用いて i 番目のメッシュの磁束を求めると

 

                       2 / 1 2 / 1 , e 0 2 2 / 1 1 p 2 / 1 1 p p 1 1 i i i i i i i i r r j r r r      (2.14) のように計算することができる．反復法(SOR 法)を用いて収束計算を行う．プラ ズマ電流より求まる磁束_pから磁場B ・電場_zp E_pはそれぞれ r r B_zp p   1  (2.15) t r E_p p       1 (2.16) と求めることができる． このように求めた磁場と電場を足しあわせ， zp zex z B B B   p ex E E E_  _  _ として代入し，電子の軌道計算を行う．

2-3-3 slowing-down collision

プラズマ中の粒子は，バックグランド粒子とのクーロン衝突により速度や軌 道を変化させる．ここでは，Slowing-Down collision とよばれる摩擦による速度 の減速現象について説明する．

プラズマ運動論を用いて，slowing-down collision の特性を議論する．そこでま ず，Landau’s Formula による衝突項を持つ，Vlasov 方程式

 

f

 

t f

 

t m m m L t f t S S s s S s S S S , , ' 1 1 d , d d ' 3 ' ' u v u v g gg g 1 u v v _                        









'



' ' S_, S S S S f f C



 (2.17) u v g          ln , 8 min max 2 0 2 ' 2 ' k k e e LSS s s    によって支配される分布関数 fS

 

v,t は Boltzmann の H 定理を満たしているこ とを示そう．

 

,t 0 fS v であるので，H 関数を

 

t f

 

t f

 

t H S S S , ln , dv v  v 



(2.18) と定義すると，次のような式を得ることができる．

 



 



 

 

f

 

t f

 

t m m t f m L t f t t f t H t S S S S S S S S S S S S S , , 1 1 , d d , d d 1 , d d d 3 , v v u v g gg g 1 v v u v v v v ' ' ' ' ln ln                             







 (2.19) ここで，vとu，Sと とを交換することによって，(2.19)式と同等の式を得るS' ことができる．その２つの式を和と取り，両辺を２で割ると

 

 

 

_            





f t m t f m L t H t S S S S S S S S , ln 1 , ln 1 d d 2 1 d d ' ' ' ' , u u v v u v

 

t f

 

t f m m S S S S , , 1 1 3 _{v u} v v g gg g 1 _' '                         (2.20) となる．ここで，

 

f

 

t m t f m S S S S SS ln , 1 , ln 1 u u v v X ' ' ' _      とすると，(2.20)式の被積分関数は

 

 

   

t f t f m m t f m t f m S S S S S S S S , , 1 1 , 1 , 1 3 v v u v g gg g 1 u u v v ' ' ' ' ln ln                                  

   

' ₃ ' ' , , S _{SS SS} S t f t f X g gg g 1 X u v _          

   







2



' 2 ' 2 3 ' ' 1 , , S _{SS SS SS} S t f t f g X g X g X u v    

 

, 

 

, 1₃







2 0  ' _' SS S S t f t f g X g u v (2.21) となる．したがって(2.20)と(2.21)式から

 

0 d d _ t H t (2.22) であることがわかる．すなわち，(2.18)で定義された H-関数は常に減少してい く．これを H 定理と呼ぶ．(2.21)式より，dH

 

t /dt0となるのは任意のS,S'とv,u に対して





 

f

 

t m t f m S S S S , 1 , 1 u u v v u v q ' ' ln //        (2.23) の場合のみである．(2.23)式は次式のようにも表せる．

 

 

u



v u

 

v u



u v v         , , ln 1 , ln 1 A t f m t f m S S S S ' ' (2.24) この条件(2.24)式は shifted-Maxwell 分布

 





_         2 2 exp , S S B S S S T k m A t f v v V (2.25)

において，T_S T_S'及びV_S V_S'の場合，満たされる．結局，H 関数(2.18)式は (2.17)式で表されたクーロン衝突によって時間とともに減少し，tにおいて v v v     S_' S S_' S T T T , の shifted-Maxwell 分布となって定常値になる． (2.17)式の衝突項は積分を含んでいるため，分布関数を決定するためには，連 立積分方程式を解かなければならないため，実際の応用には不便である．そこ で，プラズマの分布関数は時間的に Maxwell 分布に近づくことを用い，field particle の分布関数が Maxwell 分布からあまり大きくずれていない場合を想定し， (2.17)式の簡単化を行う．まず，(2.17)式において

 

_              2 ' 2 3 ' ' 2 ' exp 2 ' , u T k m T k m n t f S B S S B S S S  u (2.26) とおく．次に，







u v v v u v u v u v u v 1 _           3  (2.27) であることを利用して(2.26)を(2.17)に代入して，uに関して （vの方向に極 を持つ極座標系に変換して）積分すると，

 

_

                              ' 3 3 2 2 0 2 2 2 ' 8 , dt d S S S B S S B S S m T k m T k m e e n t f       2 2 v v v vv v vv v 1 v v 

 

t f T m S , v v v _         B k (2.28)

となり，Linearized collision Operator が得られる．ここで，

 

 

2 v ' S B S t T k m x x x t t e x 2 d d d 2 x 0    



 _     (2.29) である．(2.28)式は積分を含まないので，衝突現象を記述する基礎方程式とし て多く用いられる．ここで，(2.28)式において分布関数f

 

v,t を

 

t



 

x



f v,  vV t e (2.30) とする．_V

 

_t はテスト粒子の速度を表し，その方向に x 軸をとる．(2.28)式の 両辺に_{v をかけて，両辺を速度空間で積分すると， x}

 





 

 

 

t t V t V t d d d d d _{ } 

_

v_x v_x  v_y  v_z v 左辺





                              S j S S B S S B S S m T k m T k m e e n _{   }  _  3 2 j x 3 j x j 2 x x v v v v v x ' v v v 2 8 02 2 2 2 v v v 右辺

 



V t



 

 



j x y z



T m z S , , d d d              z y x y x B j j v v v v v v k v v   

   



_          S S S S t V V x m m m e e n 3 2 2 0 2 2 1 4    (2.31) となる．よって，次の式が得られる．

 

t V

 

t V t sl d d _ (2.32) ここで，

 



        S S S S x V m m m e e n _    _{2 2 3} 0 2 2 sl 1 1 4

 



_                S B S S S B S S S S T k V m x x x T k m m m m e e n 2 2 1 4 2 2 3 2 3 2 2 0 2 2 ,    (2.33) である．この_slは，温度T の Maxwell 分布したプラズマ中に存在するテスト_S

粒子の速度の e-folding time の逆数で，slowing-down collision frequency と呼ぶ． また，(2.33)式で与えられる slowing-down collision frequency はV2 2k_BT_S m_S

の場合は

 

x 1であり，反対にV2 2k_BT_S m_S では

 

32 3 4 x x    である．よ って，(2.33)式は 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.34) すなわち，プラズマの熱速度(v_th  2T /_S m_S )より充分に大きな速度をもつテ スト粒子の_slは，自身の速度の３乗に逆比例する．反対に熱速度より充分に小 さいテスト粒子は熱速度の３乗に逆比例した一定値に近づく．

2-4 FRC プラズマの平衡状態 研究目的において FRC プラズマは，平衡状態を仮定しシミュレーションされ ることが多いと述べたがこの節ではプラズマ平衡の求め方について説明する． 2-4-1 MHD 平衡 核融合を起こすには，プラズマを高温，高密度にして容器から離してプラズ マを保持しなければならない．その質量の運動または状態変化を示す方程式に おいて時間的変化の項がゼロであるとき，その物質は平衡状態にある．磁気閉 じ込め方式による平衡状態はプラズマの圧力勾配とローレンツ力がつりあった 状態である．この節ではプラズマの平衡状態を表す式を理想 MHD 方程式から導 く． 2-4-2 理想 MHD 方程式 理想 MHD 方程式とは電気抵抗が無視できる高温プラズマの状態を表した方 程式である．理想 MHD 方程式は以下のように記述される．

 

0      v   t (2.35)



v



v J B v _{  }       _{ }   p t  (2.36)





0       _{ }   _ p t v (2.37) 0   v B E (2.38) t       E B (2.39) j B₀   (2.40) 0    B (2.41) ここで，E,B,jはそれぞれ電場強度，磁場強度，電流密度を表している．また， 0 , , , ,    v p はそれぞれプラズマの質量密度，速度，圧力，比熱比，真空透磁率 を表している．

2-4-3 Grad-Shafranov 方程式 プラズマの平衡は一般的に Grad-Shafranov 方程式を解くことにより求められ る．この節では(2.35)~(2.41)の理想 MHD 方程式から Grad-Shafranov 方程式を導 出する． 平衡状態では時間に関する項は 0 となるので t0とし，また静的な平衡と 考えるのでu0とする．この条件により理想 MHD 方程式は次のように整理す ることができる． B j  p (2.42) j B₀   (2.43) 0    B (2.44) ここで，(2.42)式の両辺を B で内積すると 0   p B (2.45) となる．FRC プラズマは軸対称な配位を持っているので，円筒座標系(r,,z)を 用いる．(2.45)式から等圧面上に磁力線が沿うことから，等圧面上に乗る磁力線 の作る面を考えることができる．この面を磁気面と呼び，その磁気面を表す関 数をで表す．磁束関数 の定義は



 rrB_z r 0 d  (2.46) で与えられる．式(3．84)の両辺を r で微分すると r r B_z   1  (2.47) 磁場 B をベクトルポテンシャルで表すと



r r z z



z r A A A z r r e e e e e e A B                         _ 1 (2.48)

ここで，e_r 



cos,sin



，e_ 



sin,cos



であるから

r r d d d d e e e e       ,  であることに注意すると，(2．48)式は次のように表せる． z r z r r z r A A r r A r A z A z A A r e e e B       _                                     1 1 (2.49) モデルは円筒座標系で軸対称なので



 



0である．また，磁場のθ 成分は 平衡状態のとき存在しないのでB 0である．(2.49)式を整理すると

 

_z r rA r r z A e e B                     1 _{ (2.50)} となる．(2.50)式のｚ成分より，

 

rA_ r r B_z   1 (2.51) となる．(2.47)式と(2.51)式を比べると   rA (2.52) であることがわかる．また，(2.50)式のｒ成分と(2.52)式より z r B_r     1  (2.53) となる． 次に，(2.43)式の θ 成分をとり，(2.47)と(2.53)式で書き表すと， r B z B J r z        ₀                      1 1 2₂ z r r r r r   (2.54) となる．一方，(2.42)式の r 成分をとると r p B J _z     となり，(2.47)式を代入することによって，   d dp r j  (2.55) となる．よって，(2.54)式と(2.55)式より次のような方程式を得ることができる．     d d 1 2 0 2 2 p r z r r r r               (2.56) (2.56)式は FRC における Grad-Shafranov 方程式であり，境界条件と圧力勾配を 与えることによって Grad-Shafranov 方程式を解き，磁束関数 を求めることがで きる．

2-5 PIC 法 本研究では，計算領域をセルに分割し，格子点上の物理量を求めるシミュレ ーション技法を採用している．セル内部に含まれる超粒子の位置や速度情報を 格子点に反映させ，粒子の情報を流体の情報に変換するのが目的である．ここ では，その手法の一つである PIC(Particle In Cell)法について説明する．本研究で は円筒座標系で PIC 法を用いているが，簡単のため．まずはデカルト座標で説 明する．デカルト座標系における PIC 法の概念図について Fig. 2-7 に示す． Fig. 2-7 は２次元平面における１つのセルを表しており，グリッド点をそれぞ れ G1～G4 と設定する．また，１つのセル全体の面積を S とし，セル内に配置 した超粒子 P によって区切られた４つの区間をそれぞれ A，B，C，D とする． このとき，グリッド点 G1～G4 に振り分けられる超粒子 P の情報は，自身によ って区切られた４つの面積によって次のように決められる． 2 1 S C n_G  (2.57) 2 2 S D n_G  (2.58) 2 3 S A n_G  (2.59) 2 4 S B n_G  (2.60) ここでn ～_G1 n_G4は各グリッド点における密度であり， は１つの超粒子が持 つ重みを示す．各グリッド点における流束 についても同様に求めることができ， 以下のように決められる． 2 1 S C p G v    (2.61) 2 2 S D p G v    (2.62) 2 3 S A p G v    (2.63) 2 4 S B p G v    (2.64) ) )( (x x y ₁ y A  _{i j}  ，B (x_i1x)(y_j1y)， ) )( (x_{i 1} x y y_j C _   ，D(xx_j)(yy_j)，S xy

ここで，v は超粒子 P が持つ速度のトロイダル成分である．(2.57)～(2.64)式_p はセル内に１つの超粒子が存在する場合での計算方法であるが，実際には多数 の超粒子がセル内に存在する．したがって，セル内のすべての粒子においても 同様の計算をする必要があり，N 個の粒子が存在する場合は次のようになる．



  N i i G S C n 1 2 1  (2.65)



  N i i G S D n 1 2 2  (2.66)



  N i i G S A n 1 2 3  (2.67)



  N i i G S B n 1 2 4  (3.68)



  N i i pi G S C 1 2 1 v   (3.69)



  N i i pi G S D 1 2 2 v   (3.70)



  N i i pi G S A 1 2 3 v   (3.71)



  N i i pi G S B 1 2 4 v   (2.72) 全粒子の軌道計算後に(2.65)～(2.72)式を用いることで，粒子の情報を流体情報 へと変換できる．以上がデカルト座標系における PIC 法の概要である．

Fig. 2-7 デカルト座標系における PIC 法の概念図 次に本研究で用いている円筒座標系における PIC 法について説明する．円筒 座標系における PIC 法の概念図については Fig. 2-8 に示す．基本的にはデカルト 座標系での計算と同じであるが，デカルト座標系が面積比を用いているのに対 し，円筒座標系では体積比を用いている．超粒子 Q によって区切られた区間を それぞれ E，F，G，H とすると，超粒子 Q の情報は各グリッド点に次のように 振り分けられる．



  N i H V G n 1 2 1  (2.73)



x ,

i

y

j





x

i1

,

y

j





x

i

,

y

j1





x

i1

,

y

j1



y

x

D

C

B

A

x



y



y

x

1  j

y

j

y

i

x

x

_i1



  N i H V H n 1 2 2  (2.74)



  N i H V E n 1 2 3  (2.75)



  N i H V F n 1 2 4  (2.76)



  N i pi H V G v 1 2 1   (2.77)



  N i pi H V H v 1 2 2   (2.78)



  N i pi H V E v 1 2 3   (2.79)



  N i pi H V F v 1 2 4   (2.80) ) )( (r_j r z z_j E  _12  2  ，F (r_j12r2)(z_i1z)， ) )( (r2 r 2 z ₁ z G  _{j i}  ，H(r2r_j2)(zz_i)，V (r_j12r_j2)(z_i1z_i) Fig. 2-8 円筒座標系における PIC の概念図



z ,i rj





z ,i1 rj





zi,rj1





zi1,rj1



r

z

H

G

F

E

j r j

r

1  j

r

z

r

2-6 モンテカルロ法 モンテカルロ法は，フォン・ノイマンとユラムによって1945 年頃から用いられはじ めたものであり，当初は「決定論的な数学の問題を乱数を用いて解くこと」が目的であ った．しかし，最近では精度の高い乱数が得られることにより，それを用いたシミュレ ーションが盛んとなり，現在では，それをモンテカルロ法と呼ぶ．モンテカルロ法は， 確率論的モデルに基づいた確率過程を問題とするのである．モンテカルロ法では，現象 のモデルを構築し，その様子を調べようというものである．例としては，形状の複雑な 領域の面積の計算などがある．シミュレーションでは，乱数を発生しながら確率過程を 計算機上で試行していく．その最も基本的なアルゴリズムを以下に示す． Fig. 2-9 モンテカルロ法のアルゴリズム[16] コンピュータ上で発生させた乱数は先にも述べたように，一様に理想的な乱数ではな く，疑似乱数と呼ばれるものである．疑似乱数というのは，理想的な一様乱数を模倣し たものであるので，やはり完全に一様なものではない．そこで，性質の良い乱数をいか に多数発生させることができるかどうかが，モンテカルロ・シミュレーションにおいて 良質の結果を得るためのポイントとなる． あるイベントが起こる断面積をとすると，そのイベントが t の間に起こっ たとする判定方法は，0~1までの一様な擬似乱数を用いて

 

t ξ n v v  (2.81) と計算することができる．ここでnはプラズマ密度，vは中性粒子とプラズマ の相対速度である．ここで注意するのは t のとり方である．これを大きくする と(2-1)の左辺は 1 を超えてしまう．このとき t は長く取りすぎており，1 を超 えた分だけ起こる回数を過小評価することになる．そのため，(2.81)の左辺は 1 よりも小さくなるように注意しなければならない．

2-7 衝突断面積 モンテカルロ法による粒子の衝突散乱を判定するためには，粒子の断面積が 判っていなければならない．ここで衝突断面積について説明する． 衝突断面積とは，ビーム粒子が他の粒子と衝突して散乱するために通らなけ ればならない面積のことである．この面積は Fig. 2-10 のように粒子の断面積か ら， 2 4 a   と与えられる． Fig. 2-10 衝突断面積 今回のモデルにおいて，衝突に関係する粒子は水素の中性粒子とイオンであ る．このため，水素の原子半径と水素イオンの原子半径から，衝突断面積を導 出できる． 2-8 正規分布 マクスウェル分布とは，熱力学的平衡状態において，気体分子の速度が従う 分布関数である．この分布を再現するために正規分布を用いた．正規分布は Fig. 2-11 下図のようにマクスウェル分布に非常に似た曲線を描き，近似可能である． ここではその正規分布について説明し粒子の初速度の与え方を定義する．

正規乱数とは正規分布に従った母集団を持つ乱数のことを表す．具体的には 確率密度関数が

 

 2 2 2 2 1 _       x e x f ，



x



（2.82) で表される確率分布を指す．この関数をグラフに表すと Fig. 2-12 のようにな る．このときの正規分布曲線の形状と位置はと の値によって完全に決まる． は曲線の中心を与えるものであり，  は曲線の左右への広がりを与えている． したがってを平均値， を標準偏差，2_{を分散という．この正規分布を持つ} 独立な確率変数の系列の実現値を，正規乱数列という．この正規乱数列におい て平均値に近い部分の値が多く出て，から離れるほどその出方はすくなく なる事は Fig. 2-12 から明白である． Fig. 2-12． 正規分布の確率密度関数［16］ 式（2.82)において     x Z と変換すると，

 

2 2 2 1 z e x g    ，



z



（2.83) となり，これは平均値0，分散2 1の正規分布と考える事ができる．こ れを基準型正規分布と呼んでいる．その曲線は Fig. 2-13 のようになる．

Fig. 2-13 基準型正規分布の確率密度関数［16］ 一般にこの基準型正規分布をもつ乱数列を生成する事で平均値，分散2_の 正規分布を持つ関数は式（2.83)より     _n n z x



n1.2.3...



（2.84) と変換する事で容易に作成が可能である． 2-9 ボックスミューラー法 マクスウェル分布を正規分布を用いて表せる事を前述したが，その正規分布 をコンピュータ上で再現しなければならない．正規分布を作成する方法として， 一様乱数を生成し，ボックス＝ミューラー法を用いて一様乱数を正規乱数に変 換する方法が一般的である．ここではボックス＝ミューラー法によって一様乱 数から正規乱数の生成について説明する． ある乱数₁，₂が標準正規分布の確率変数であると与えられたとき，それぞ れx1 xdx，y2  ydyになる確率は dx x       2 2 1 exp 2 1  ， y dy     2 2 1 exp 2 1  （2.85) と表せる．また，それが同時に成立する確率は両者の積であり，p ,

 

x y をそ の確率密度関数であるとすると

 

x y dxdy



x y



dxdy p      _{ }  2 2 2 1 exp 2 1 ,  （2.85）

である．ここで，x と y を極座標に変換すると xrcos，yrsinであ り， 2 2 2 r y x   ，







 



          drd drd rshin r drd r r drd y r y x r x dxdy                2 2 cos cos sin sin cos （2.86） である．したがって正規分布関数をする二つの乱数を₁ _r cos_，     ₂  sin としたときrr rdrかつ dとなる確率は

 

    drd r r drd r p        2 2 1 exp 2 1 , （2.87） と表せる．ここで，逆関数法を用いて_r，_ を₁，₂で表す． 1  ，₂が一様乱数であるなら，  

 

                                           



 



2 0 2 0 2 0 2 0 2 1 2 1 exp 1 2 1 exp 2 1 exp 2 1 2 1 exp 2 1 , r r r r r r dr r r dr d r r drd r p P             （2.89）

 

 

                  2 1 2 1 2 1 exp 2 1 2 1 exp 2 1 , 0 0 2 0 0 2 2                        





 



  d d r drd r r drd r p P （2.90） 両者から_r，_ を求めると



1 1



log 2  _r    （2.91） 2 2 1 _  _  （2.92） これを₁，₂の式に代入して



1



2

1  cos 2log1  cos2

  _{r }    ， （2.93）



1



2

2  sin 2log1  sin2

  _{ }    （2.94） が得られる．ここで1₁も一様乱数であるため，11 1と置き換えて 2 1 1 2log cos2    （2.95） 2 1 2 2log cos2    （2.96） となる．これがボックスミューラー法であり，一様乱数₁，₂から正規乱数₁， 2  を得るための方法である．

2-10 有限差分法 ここでは差分法について説明する．差分法で微分方程式を解くプロセスは， その方程式が持つ物理的イメージを対応させやすいので，数値計算の視点から もとの微分方程式の意味を理解するという点で有効である． 偏微分係数は次式のように定義できる．



  

x y x f y x x f x f x          , , lim 0 (2.97) 偏微分係数が(2．97)式のように定義されるとき十分小さな x に対しては



  

x y x f y x x f x f        , , (2.98) と近似することができる．この近似法を有限差分法(finite difference method)と 呼ぶ． 差分法の定量的な性質について議論するときは Taylor 展開を利用する．2 変数 関数 f

 

x,y の場合，点



x₀, y₀



を中心とする Taylor 展開は次式で与えられる．



  

f

 

a b y y x x n b a f y b x a f n n , ! 1 , , 1



                    (2.99) 今，例としてy0の場合を考えると(2．99)式は



  

   

   

 

                2 2 ₂ 2 3 ₃, ! 3 , ! 2 , , , x b a f x x b a f x x b a f x b a f b x a f と展開されるから，これをf xについて解けば

 



  

 

  

 _               6 , 2 , , , , 2 3 3 2 2 x x b a f x x b a f x b a f b x a f x b a f (2.100) を得ることができる． なお，項数が増えていくと，離散化された値を f



ax,by



のように書く のは大変なのでこれを下付きの添え字を用いることにする．ある点



x₀, y₀



を基 準として各々の点に番号をふり，



, 0,1,2



, ₀ 0       x i x y y j y i j x_{i j} のように表す．また，対応する関数値や微分係数を

 

 

x f x y x f f y x f i j i j i j ij       , , , , ,

と表現すると(2.100)式は

 

              ①               _ 6 2 2 3 , 3 2 , 2 , , 1 , x x f x x f x f f x f_{i j i j i j i j i j}

 

x O x f fi j i j _{ }    1, , とあらわすことができる． 2-11 反復法（SOR 法） 反復法とは解を推測する事から始めて厳密解に収束するまで，あるいは厳密 解に十分近くなるまで繰り返し計算を実行する数値計算法である．本研究では Grad-Shafranov 方程式の解を求めるときに SOR 法を用いている．

SOR は Successive Over-Relaxation の略である．SOR アルゴリズムは反復法の 一種で Gauss-Seidel 法として知られる方法の精密化である．なお Gauss-Seidel 法 自体は Jacobi 法を発展させたものである．SOR 法を説明するにはこの 2 つの方 法を最初に述べることにする． Jacobi 法を説明するのにn個の未知数をもつn個の一般の連立方程式を例にあ げる． 1 1 2 12 1 11x a x a x b a   _{n n}  2 2 2 22 1 21x a x a x b a   _{n n}  n n nn n n x a x a x b a _{1 1} _{2 2}  このような連立方程式に対して第 1 番目の式をx₁の式に直すと 11 1 2 12 1 1 a x a x a b x     n n であり，これよりx₁の第 k 回目の反復値をx₁ k とすると       11 1 2 12 1 1 1 a x a x a b x k n n k k    

が導かれる．同様に第 2 番目の式を変形して，x に対する反復公式得られる．₂ そして一般にx に対する反復公式は_i i 番目の式から得られ，           i i k n n i k i i i k i i i k i i k i a x a x a x a x a b x 1   1 1   1 1  1 1  となる．これが一般的な連立方程式に対する Jacobi 法である． Jacobi 法では初期値    0  0 2 0 1 ,x , ,xn x で計算を始め，これを用いて次の近似解    1  1 2 1 1 ,x , ,xn x を計算する．しかし，  1 2 x の計算においてすでにx が得られているのに₁ 1 x₁の 初期値  0 1 x を用いることは無意味である．同様の理由でx を計算するときには₃ 1 一つ前の反復で得られた値  0  0 4 , ,xn x とともに新しい値x と₁ 1 x を用いること₂ 1 が賢明である．この考え方は解ベクトルの各々の成分に適応され Jacobi 法の修 正版，すなわち最新の近似解を用いる Gauss-Seidel 法を与える．したがって Gasuu-Seidel の反復は           i i k n n i k i i i k i i i k i i k i a x a x a x a x a b x               1 1 1 1 1 1 1 1 1 で与えられる． しかし Gauss-Seidel 法はnの数が大きくなってくると収束速度が遅くなる性質 がある．そこで収束速度を上げるための方法として SOR 法がある．この方法は 次のような考え方に基づいている．すなわち Gauss-Seidel 法を用いて，新しい値  k1 i x が生成されれば新しい値と前の値の重み付き平均を計算することによりさ らによい値が得られるというものである．ゆえに  k1 i x は

   

_

_

 k i k i k i x x x 1  1  1 (2.101) で与えられる．ここで

 

0 は k と i に無関係な任意のパラメーターである． 1   であるならば  k1 i x は Gauss-Seidel 法により得られる値であることに注意さ れたい．1のとき(2.101)式は SOR 法を定義する．

第 3 章 電子挙動シミュレーション結果と考察

本章では，第 2 章で記述した粒子軌道計算コードを用いて行ったシミュレー ション結果を示す．シミュレーション結果の妥当性を評価するために前章で記 述した平衡計算法 Grad-Shafranov 方程式から導いた装置中央部における 1D 平衡 プラズマの各種径方向分布との比較を行う． 計算が進むと場の数値が大きくなるため，粒子計算を解くアダムス法が収束 しづらくなる．このためシミュレーションタイム最後まで計算が終了していな い． 3-1 シミュレーション結果 Fig. 3-1～3-6 に各シミュレーション結果を示す． Fig. 3-1 に粒子位置の時間発展を示す．2 章で述べたように電子の初期配置は 装置半径r で規格化した計算領域内に一様に分布させている． _w

Fig. 3-1-1(a)10step から(f)60step では見かけ上大きな変化はなく，少しずつ粒子 密度が薄くなっていく様子が分かる．シミュレーションの初期では外部コイル が作る磁場が支配的であり，その磁場勾配により電子が装置外部にむかって加 速され計算領域から出てしまい粒子が損失していく．(i)90step で r=0.9 付近に粒 子が少し線状に集中していき，(j)100step では装置中央をくりぬいたドーナツ状 に粒子が径方向に圧縮されてきている．(k)500step から Fig. 3-1-2 (c)2500step で は段々と電離を起こしながら，粒子が装置中央へ圧縮されていく様子が確認で きる．そのため 500step 以降粒子損失が発生しにくくなる．なお粒子の運動方程 式は(2.1)式でといており R 方向の圧縮が起きる条件は





0 2 2    



i r ei e z e m q r v v B v r _   _となる．