粒子－流体混合モデルによるFRCプラズマの自発的トロイダル回転現象に関する研究

平成２１年度 修 士 論 文

粒子−流体混合モデルによる FRC プラズマの

自発的トロイダル回転現象に関する研究

指導教員 高橋 俊樹 准教授

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

電気電子工学専攻

池田 達也

目

次

第 1 章 序論

・・・・・・・・・・・・・・・・・・・・・・・・・・・1 1-1 核融合 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・1 1-2 磁場反転配位プラズマ ・・・・・・・・・・・・・・・・・・・・・3 1-3 ドップラーシフト計測 ・・・・・・・・・・・・・・・・・・・・・6 1-4 研究目的 ・・・・・・・・・・・・・・・・・・・・・・・・・・・8

第 2 章 プラズマの回転理論

・・・・・・・・・・・・・・・・・・・9 2-1 従来のプラズマ回転説 ・・・・・・・・・・・・・・・・・・・・・9 2-1-1 粒子損失によるプラズマ回転説 ・・・・・・・・・・・・・・・9 2-1-2 電場短絡によるプラズマ回転説 ・・・・・・・・・・・・13 2-2 磁束減衰に伴うプラズマの回転説 ・・・・・・・・・・・・・18 2-2-1 正準角運動量からの予測 ・・・・・・・・・・・・・・・18 2-2-2 運動方程式からの予測 ・・・・・・・・・・・・・・・・19 2-2-3 ステルマー領域からの予測 ・・・・・・・・・・・・・・20

第 3 章 粒子軌道シミュレーションモデル

・・・・・・・・・・24 3-1 計算モデル ・・・・・・・・・・・・・・・・・・・・・・・・24 3-2 FRC プラズマの平衡状態 ・・・・・・・・・・・・・・・・・・25 3-2-1 MHD 平衡 ・・・・・・・・・・・・・・・・・・・・・・・25 3-2-2 理想 MHD 方程式 ・・・・・・・・・・・・・・・・・・・25 3-2-3 Grad-Shafranov 方程式 ・・・・・・・・・・・・・・・・・・25 3-3 イオンの初期分布 ・・・・・・・・・・・・・・・・・・・・・・27 3-4 粒子軌道計算 ・・・・・・・・・・・・・・・・・・・・・・・・29 3-4-1 磁束の減衰と場の関係 ・・・・・・・・・・・・・・・・・・29 3-4-2 場の計算法 ・・・・・・・・・・・・・・・・・・・・・・・30 3-4-3 クーロン衝突 ・・・・・・・・・・・・・・・・・・・・・・32 3-4-3-1 slowing-down collision ・・・・・・・・・・・・・・・・・32 3-4-3-2 ピッチ角散乱 ・・・・・・・・・・・・・・・・・・・・36 3-4-4 PIC 法 ・・・・・・・・・・・・・・・・・・・・・・・・・38 3-4-5 平滑化 ・・・・・・・・・・・・・・・・・・・・・・・・42

第 4 章 ドップラーシフト計測の妥当性

・・・・・・・・・・・46 4-1 実験結果とシミュレーションの比較 ・・・・・・・・・・・・・46 4-2 重水素イオンと不純物イオンの流速の比較 ・・・・・・・・・・49 4-2-1 オームの式を簡略化した場合のイオン流速 ・・・・・・・・49 4-2-2 ローレンツ力及び圧力勾配を考慮した場合のイオン流速 ・・50 4-3 考察 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・52 4-3-1 測定点による違い ・・・・・・・・・・・・・・・・・・・52 4-3-1-1 反磁性ドリフト ・・・・・・・・・・・・・・・・・・53 4-3-1-2 イオン密度分布 ・・・・・・・・・・・・・・・・・・55 4-3-2 電場による違い ・・・・・・・・・・・・・・・・・・・・56

第 5 章 結論

・・・・・・・・・・・・・・・・・・・・・・・・・58

参考文献

・・・・・・・・・・・・・・・・・・・・・・・・・・・・59

謝辞

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・61

第 1 章 序論

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

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

力線のカゴの形や電磁石の配置によってさまざまなタイプがある。その代表的 なものに「ミラー型」「ヘリカル型」「トカマク型」の三つのタイプがある。こ の中でもトカマク型は将来の核融合炉に最も有力とされている。2016 年頃の完 成を目指し、フランスのカダラッシュに建設することになった ITER（国際熱核 融合実験炉）でも採用されている。 本研究では、ミラー型に分類される磁場反転配位と呼ばれる閉じ込め方式を 使ったプラズマに関する研究である。 1-2 磁場反転配位プラズマ 磁場反転配位（FRC : Field-Reversed Configuration）プラズマは、Fig. 1-4 に示 されるような磁場構造をしたプラズマである。FRC は開いたミラー磁場中に、 閉じた磁力線領域が形成されており、この二つの領域の境界を separatrix と呼ぶ。 また、中心領域に磁場がゼロとなる場所（磁気中性点 : O-point）が存在し、そ の付近に高温のプラズマを閉じ込めることができる。開いた磁力線によっての みプラズマを閉じ込めるミラー方式では、開放端からの粒子損失率が大きく、 またベータ値も低い。しかし、このような閉じた磁力線領域を持つ FRC は、磁 場閉じ込め方式の中で最も高いベータ値を持つ。ベータ値はプラズマ圧と磁気 圧との比を表し、一般的には （プラズマ圧）／（外部磁気圧） (1.1) で定義される。直線的な磁力線に沿って磁場の変化がない場合、 const. 2 ₀ 2 B p (1.2) Fig. 1-3 磁場閉じ込め方式[1]

となり、FRC の場合、磁気中性点が存在するので 0 0 2 0 2 ex 2 2 p B p B (1.3) が成り立つ。ここで、B は外部磁場の大きさを表し、_ex p は磁気中性点での圧力₀ を表す。つまり、磁気中性点でのベータ値は、 1 2 ₀ 2 ex 0 B p (1.4) となり、外部磁気圧はそのまま磁気中性点でのプラズマ圧力となる。ベータ値 が高い場合、ベータ値が低いものと比較すると、ある一定の外部磁場のときに、 より高い圧力のプラズマを閉じ込めることができる。つまり、FRC は他の閉じ 込め方式より弱い磁場でプラズマを閉じ込めることができるので、電力の消費 量が少なくて済み、経済的な閉じ込め方式であると言える。 高ベータプラズマである FRC プラズマは、アドバンスド核融合炉心プラズマ としての可能性を持っている。前節で触れた ITER で、実用化を検討されている 核融合反応は、重水素と三重水素による反応（D-T 反応）である。Fig. 1-5 に核 融合反応断面積の衝突エネルギー依存性を示す。D-T 反応は他の反応に比べ、反 応断面積が大きく、ピークが高い。そのため、他の反応より容易に核融合反応 が起きる。しかし、14.06MeV の高速中性子が発生するため、核融合炉の構造材 が放射化されてしまう。このため、核融合炉材料の耐久年数が短く、炉壁の交 換、放射化した材料の廃棄などに高いコストが必要とされる。それに対し、ア ドバンスド核融合反応として、重水素とヘリウム 3 による反応（ 3 He D 反応） がある。この反応は D-T 反応と比べ反応断面積が小さく、核融合反応が起こり にくいが、高速中性子が発生せず、代わりに 14.67MeV の陽子が発生する。陽子 は電荷を持つため、磁場によって捕捉されるので、外部磁場を制御することで 核融合炉の損傷を防ぐことができる。また、磁場によって陽子をエネルギー変 換機に導き、陽子の持つエネルギーを電気エネルギーに変換することが可能で ある。 3 He D 反応は炉壁の交換、放射性物質の廃棄にコストがかからず、陽子 から直接エネルギーを得られるため、経済的な反応である。 このように、高ベータプラズマである FRC プラズマは、現在実用化が検討さ れているものよりも、経済性が高く魅力的である。しかし、プラズマの配位維 持時間が短く、FRC を実用化するためには配位維持時間を約 10 倍から 100 倍に 伸ばす必要がある。現在、FRC プラズマの配位維持時間が短い理由は明らかに されていない。しかし、可能性として （1）磁場揺動が大きくこれにより粒子が輸送されるため粒子閉じ込めが極 めて悪い。 （2）異常抵抗によりプラズマ電流が減少し、閉じ込め磁場が減衰する。

（3）セパラトリクス内部で回転不安定性が起こり、配位が崩壊している。 などが考えられている。FRC プラズマの配位維持時間を伸ばすには、これらの 可能性を検証し、その現象が起こるメカニズムを調べ、改善する必要がある。 10-4 10-3 10-2 10-1 100 101 100 101 102 103 [b a rn ] E [keV] D-T D-D D-3He

O-point

開いた磁場 閉じた磁場

X-point

Separatrix

z

r

Fig. 1-4 FRC の磁場構造 Fig. 1-5 核融合反応断面積

1-3 ドップラーシフト計測 前節で、FRC プラズマの配位維持時間が短い理由の一つとして、回転不安定 性があると説明した。しかし、この回転不安定性を誘起する回転起源について は未解明である。そこで、回転不安定性を解消し、配位維持時間を伸ばすため に、回転そのものを調べる必要がある。 FRC プラズマの回転を調べるためには、プラズマを構成している重水素イオ ンの流速分布を計測する必要がある。そのための計測方法がドップラーシフト 計測である。ドップラーシフトの概念図を Fig. 1-6 に示す。また、Fig. 1-7 はド ップラーシフト計測の模式図である。FRC プラズマ内のイオンは電子と衝突し、 励起される。そして、基底状態に戻るとき、可視光領域付近の特定波長の線ス ペクトルを放射する。この放出されたスペクトルを光学系により検出すること で、流速分布を計測することができる。イオンの流速は以下の式で表される。 c c u 0 0 0 i (1.5) ここで , ₀,cは、それぞれドップラーシフト後の波長、ドップラーシフト前の 基準波長、光速である。スペクトルの基準波長は既知であるので、計測により 求められたドップラーシフトの大きさよりイオンの流速を求められる。 しかし、ドップラーシフト計測では重水素イオンの流速を計測することはで きない。なぜなら、核融合の節で説明したように、プラズマ内は非常に高温で あるため、プラズマ内で重水素イオンは完全に電離しており、発光せず、スペ クトルを計測できないからである。そこで、ドップラーシフト計測では、重水 素イオンではなく、プラズマ内に存在する不純物イオンを用いて計測している。 その方法を以下に示す。 まず、不純物イオンについての運動方程式を立てる。 im im im im im im im im im im u u E u B R u t p n q t (1.6) ここで、 _im,q_im,u_imは、それぞれ不純物イオンの質量密度、電荷、流速である。 また、 はバックグランド粒子種を意味し、ここでは電子と重水素イオンであ る。この方程式の時間発展、電場、圧力勾配、粘性テンソルの項を無視し、 方 向のみについて考える。r 方向と z 方向の流速がないと考えると、 u_im u_imの 項はゼロになる。また、u_im Bの項も影響がほとんどないと考えられ、無視で き、摩擦力の項が支配的となる。摩擦力は u u R m n (1.7)

で表されるので、(1.6)式は 0 D im imD im im e im ime im im D e, im m n u u m n u u R (1.8) と書ける。(1.8)式を不純物イオンの流速についてまとめると、 imD ime D imD e ime im u u u (1.9) となる。電子と重水素イオンの密度・温度が等しいと仮定すると、(1.9)式は 2 1 D 2 1 e D 2 1 D e 2 1 e im m m u m u m u (1.10) と書くことができる。ここで、電子の質量は重水素イオンの質量に比べ、十分 小さいとし、 12 e m とm_e12u_e を無視すると、 D im u u (1.11) という近似式が得られる。 このような考え方によって、ドップラーシフト計測では、不純物イオンの流 速を計測し、それを重水素イオンの流速とみなしている。 スペクトルの基準波長 ドップラーシフト Fig. 1-6 ドップラーシフトの概念図

1-4 研究目的 前節で、ドップラーシフト計測では不純物イオンの流速を、重水素イオンの 流速とみなしていると説明し、その考え方を示し、(1.11)式を導いた。しかし、 その計算の途中で、様々な近似を行なっている。そのため、実際の重水素イオ ンの流速と不純物イオンの流速が異なっている可能性がある。 本研究の目的は、ドップラーシフト計測において、本当に不純物イオンの流 速を重水素イオンの流速とみなしてよいのかを、シミュレーションを用いて確 認し、ドップラーシフト計測の妥当性を評価することである。 θ z 視線方向 Fig. 1-7 ドップラーシフト計測の模式図

2 章 プラズマの回転理論

2-1 従来のプラズマ回転説 FRC の回転起源は、いまだ明らかにされていない。トロイダル方向（円筒座 標系のθ方向）の回転は、楕円変形に対して不安定であり、これは n=2 モード 回転不安定性として知られている。テータピンチ生成された FRC プラズマの支 配的な崩壊過程であり、回転起源を明らかにすることは配位維持の観点からも 非常に重要である。過去の研究論文によく引用されているメカニズムは二つあ り、一つは粒子損失説、もう一つは電場の短絡説である。ここでは、この二つ の説について簡単に説明する。 2-1-1 粒子損失によるプラズマ回転説 FRC プラズマは高ベータプラズマであり、したがって磁場の変化の特性長に 対するイオン旋回半径が大きい。この結果、セパラトリクス内イオンの旋回中 心が、輸送過程によってセパラトリクス外部に移動し、開いた磁力線に沿って 端損失する可能性が高い。また特に輸送過程を経なくとも、生成時の運動方向 によって端損失が起こりうる。その様子を Fig. 2-1 に示す。セパラトリクス近傍 の流体要素に着目して考えてみよう。もし、プラズマ生成時に等方的であると すれば、θ方向の速度成分が正である粒子数と負である粒子数は等しい。正で あるイオンは、ローレンツ力によってセパラトリクス内部に閉じ込められるが、 負であるイオンは旋回中心がセパラトリクス外部になる。これらのイオンは、 開いた磁力線によって閉じ込め装置外部に導かれる可能性が高い。このため、 θ方向の速度成分によって選択的な損失が起こる。一方、セパラトリクス内部 に閉じ込められたイオンは、θ方向に対して正の方向（Fig. 2-1 の場合、イオン 反磁性ドリフト方向）へ運動しているものとなる。つまり、コアプラズマは回 転することとなる。 いわゆるロケット効果による回転機構は、理論的にも指摘されてきた。特に、 速度空間粒子損失（VSPL: Velocity Space Particle Loss）の存在が示され、この結 果としてコアプラズマの回転が起こりうることも、FRC 研究者にはよく理解さ れている。ここで、VSPL について少し説明しよう。 無衝突で軸対称プラズマの場合、運動の恒量は全エネルギーT と正準角運動 量 P であり、 ) , ( 2 1 2 2 2 z r q m T v_r v v_{z (2.1)}

) , ( zr q mr P v (2.2) である。これらの物理量を使うと容易に、 2 2 2 2 2 2 2 ) , ( 2 1 2 1 ) , ( mr z r q P m m z r q T v_r v v_z v _(2.3) の不等式を得る。(2.3)式は、運動の恒量（定数）と位置を表す変数で記述される ことに注意する。つまり、実効的なポテンシャルエネルギーに関する式とみる ことができる。したがって、右辺はステルマーポテンシャルと呼ばれる。(2.3) の不等式は、軸対称プラズマにおける粒子の運動可能領域を記述している。も しイオン運動が確率論的（ストキャスティック）であり、運動可能領域を隈無 く動く場合、エルゴード的であると呼ばれる。逆に、周期的な運動の場合は、 二つの運動の恒量以外にも保存量が存在することになり、運動可能領域の一部 のみを動きうる。エルゴード的な性質があれば、運動可能領域の形状を調べる ことによって、粒子損失を議論することができる。代表的な運動可能領域の形 状を模式的に Fig. 2-2 に示す。閉じた運動可能領域(a)の場合、イオンはセパラト リクス内部のみを運動可能であり、軌道が外部に及ぶことはほとんどない。し たがって、端損失することはあり得ない。一方で、(b)のように運動可能領域が 閉じ込め領域内で閉じていない場合、イオン運動がストキャスティックであれ ば端部から損失してしまう。このように、ミラー端部において(2.3)の不等式を満 足すれば、運動可能領域が閉じ込め領域内で閉じないことになり、損失する可 能性がある。生成直後にイオンの速度分布は等方的であると仮定して、VSPL に よるロケット効果でどれだけの流速が得られるかを、簡単に評価してみよう。 静電場のない軸対称 FRC プラズマの中央面（z=0）において、密度分布に応じ て擬似乱数を用いて粒子を生成する。速度分布は Maxwell 分布をしていると仮 定し、さらに運動方向は等方的であるとする。後の章で説明する平衡計算から 密度や磁場の分布を計算し、さらに Maxwell 分布を想定して粒子を生成する。 ここで、プラズマ温度は一様であると仮定する。粒子を生成した位置と速度に よって、二つの運動の恒量である全エネルギーT と正準角運動量P が定まる。 さらに、ミラー端部におけるステルマーポテンシャルの極小値を計算し、(2.3) が成立するかを確認する。もし、(2.3)の不等式が成り立てば、イオンはミラー端 部で運動可能であるため端損失すると判断する。そうでないイオンは、閉じ込 め領域内に捕捉される。数値計算の結果を Fig. 2-3 に示す。これによれば、磁気 中性点においても約 3 割の粒子が損失しうることがわかる。残された７割の粒 子によって 0.3 アルフベン速度程度の流速が生じうることが示されている。セパ ラトリクス近傍では、端損失粒子の割合は急激に増加し、これによってセパラ トリクス上で約 0.5 アルフベン速度のトロイダル流速が得られる。

この結果から、粒子損失によるプラズマの回転は実験を説明しうるといえる。 なぜなら、実験で得られている流速も、ほぼ同様の 0.3 アルフベン速度程度であ るからである。生成後、徐々にプラズマ回転速度が増加するのは、粒子損失が 徐々に起こるからであると考えることもできる。Fig. 2-3 で得られている回転速 度は、生成直後の磁場や密度分布から想定される最大の回転速度であるといえ る。磁場や密度分布が時間変化する場合は、詳細な検討が必要である。 最近のハイブリッドシミュレーションによる計算結果の報告では、磁束減衰 に伴う粒子損失がプラズマ回転の原因と指摘されている。Fig. 2-4 は、磁束(a)と セパラトリクス内部のイオン粒子数(b)の時間発展を示している。磁束が減衰す るに伴い、セパラトリクス近傍の粒子の案内中心が開放端磁場領域に移動し、 結果的にセパラトリクス内部のイオン粒子数が減少している、という説明であ る。さらに、Fig. 2-5 はこの結果、セパラトリクス内部に存在するイオンの角運 動量が時間とともに増大していることがわかる。80 アルフベン時間で、0.3 アル フベン速度の回転が得られている。Beloba 等の計算結果は、時間とともに増加 する流速分布を得ている点では、実験を再現できている。一方で、磁束が減衰 している系であるにもかかわらず、全イオンの角運動量が保存することには疑 問が残る。また、全角運動量の保存から、セパラトリクス外部のイオンが持つ 角運動量は、内部のそれと符号を逆にする。つまりセパラトリクスを境界とし て、回転方向が逆転することになる。NUCTE-III の実験では、これまでのところ、 開放端磁場領域の回転速度が逆転することは計測されておらず、実験を完全に 説明できていない。 Fig. 2-1 セパラトリクス近傍のイオン運動方向による端損失

Fig. 2-2 運動可能領域の形状。(a)閉じた運動可能領域と(b)開いた運動可能領域。 開いた運動可能領域の場合、プラズマイオンはミラー端部にまで運動可能であ るため、端部からの損失が起こりやすい。 Fig. 2-3 ステルマー領域計算による回転速度の評価。青の実線は、運動可能領 域が閉じている粒子密度である。ここでn₀は磁気中性点の密度である。赤の実 線は、運動可能領域が閉じているイオンの平均速度（トロイダル方向）である。 Fig. 2-4 捕捉磁束の最大値(a)とセパラトリクス内部のイオン粒子数(b)（初期を 1 に規格化）の時間発展。Belova 等のハイブリッドシミュレーションの結果。

Fig. 2-5 角運動量(a)と流速(b)の時間発展。L_inはセパラトリクス内部のイオンが 持つ角運動量であり、L_totalは計算粒子すべての角運動量の総和である。 2-1-2 電場の短絡によるプラズマ回転説 FRC のコアプラズマのイオン反磁性方向への回転メカニズムには複数の説が 提唱されている。その一つとして有力視されているものに、電場の短絡による 周辺プラズマの回転が挙げられる。この説によれば、粘性トルクによりセパラ トリクス内部にまで回転が及ぶ。 FRC の開いた磁力線は，実験装置によってはミラー端部の導電壁に接触する ことがある。Fig. 2-6 を参照してもらいたい。スクレイプオフ領域と呼ばれるセ パラトリクス周辺領域のプラズマ電子は、開いた磁力線に沿って動きやすく、 導電壁にまで到達する。もし、磁力線間に電位差が生じていれば、電子は速や かに動き電位差はなくなるであろう。 Fig. 2-6 ミラー端部の導電壁に接触する磁力線とスクレイプオフプラズマ

一般的な電場短絡による回転理論は、イオン流体の運動方程式を用いて、以 下のように議論を進められる。 粘性と摩擦力を無視したイオン流体方程式は、 i i i i i i i i Zen p t u u E u B u (2.4) プラズマの生成直後は流れがなく平衡状態にあるとすると、イオン流体の径方 向の力のバランスは、(2.4)式から 0 i i i r p E en Z _{r (2.5)} によって表される。つまり、イオン圧力勾配のある FRC プラズマが流れのない 平衡状態にあるためには、径方向電場が存在することになる。今、この状況下 で開いた磁力線が導電壁に導かれているとき、電子は径方向電場をキャンセル するように動く。結果的に電場はなくなるが、その代わりにイオントロイダル 流速が、 0 i i i i r p B u en Z _z (2.6) を満足するように現れる。このトロイダル流速は、 r p B en Z u z i i i i 1 (2.7) であり、これはイオン反磁性ドリフト速度である。 これに対して、Steinhauer はスクレイプオフプラズマの回転を次のように議論 した。電場を決定するのは慣性力が無視できる電子のダイナミクスからであり、 (2.5)式のイオン流体の運動方程式からではない。磁力線に沿った方向の、電子の 運動方程式は 0 e || e s p E en である。電場を静電ポテンシャルで記述すると、 0 e e s p s en で、電子温度は磁力線に沿って一様であると仮定すると、 0 ln e e _n e kT s ゆえに、 const. ln _e e _n e kT (2.8) が得られる。(2.8)式を径方向に微分することで、

0 e e e r n en kT E_{r (2.9)} を得るが、この径方向電場によって、スクレイプオフプラズマはE Bドリフト する。特に中央面近傍では、磁場は z 成分のみであるから、 0 2 i z z r B B E u _(2.10) となり、このイオン流速も(2.7)式同様、イオン反磁性ドリフトの方向になる。(2.9) 式で得られた電場を、(2.10)式に代入すると、 r p B en u z e e i 1 (2.11) となる。(2.7)式と(2.11)式を比較すると、イオン圧力勾配と電子圧力勾配の違い があるだけである。FRC プラズマでは電子温度が低いため、(2.11)式で得られる 回転速度のほうが、若干小さい。 ところで、(2.7)および(2.11)式にしたがって、トロイダル方向への回転をする のは、スクレイプオフプラズマであり、セパラトリクス内部のコアプラズマで はない。プラズマの粘性に起因するトルクが発生し、セパラトリクス内部のプ ラズマをスピンアップするというのが、電場短絡によるプラズマ回転説の主張 である。 軸対称性を仮定し、トロイダル方向の流速のみがある場合、粘性力による流 速の時間発展は、 r u r r r t u , i 2 ci i i 10 3 nT (2.12) によって表される。スクレイプオフプラズマの回転がコアプラズマまでを回転 させうるか、FRC プラズマの圧力や磁場分布を考慮して計算してみよう。 Fig. 2-7 に示すように、セパラトリクス上で流速が時間とともに増加する場合 を考える。日本大学の FRC 実験装置 NUCTE-III（Nihon University Compact Torus Experiment-III）における実験結果と比較するために、トロイダル流速の加速度 をほぼ一致させ、赤色の実線で示すように 35 アルフベン時間に 0.4 アルフベン 速度まで比例増加すると仮定する。もし、粘性トルクがセパラトリクス内部を 回転するのに十分であれば、または、粘性によってセパラトリクス上の運動量 が十分に拡散すれば、セパラトリクス内部の流速分布は剛体回転に近づいてい るはずである。一方で、粘性の効果が十分でなければ、セパラトリクス近傍に 局所的に流速が現れ、その他の内部領域で流速はゼロになる。 (2.12)式を、有限差分法を用いて数値的に解いた結果が Fig. 2-8 である。容易

にわかるように、流速分布は r に比例して増加している。この流速分布は、剛体 回転分布である。したがって、周辺の回転は、粘性によって十分に内部まで伝 搬することが示された。つまり、電場の短絡によるスクレイプオフプラズマの 回転が、FRC の回転起源となりうることがわかる。 Fig. 2-7 粘性トルクによる FRC コアプラズマの回転を確認するためのモデル。 左の図は、セパラトリクス周辺のスクレイプオフプラズマが回転していること を模式的に示している。右はセパラトリクス上の流速が時間に関して比例増加 するモデルを示す。 Fig. 2-8 セパラトリクス上の流速が時間発展するときのセパラトリクス内部の トロイダル流速分布。プラズマの粘性によってトロイダル運動量が拡散する。 黒の実線は 10 アルフベン時間の分布であり、同様に青は 20、赤は 30 アルフベ ン時間の流速分布を示す。セパラトリクスの位置は，r/r_w 0.35になる。

しかしながら、セパラトリクス内部流速の時間発展はスクレイプオフプラズ マ回転の時間発展に大きく影響されるので、(2.7)式と(2.11)式で与えられる反磁 性ドリフト速度の時間発展をよく調べる必要がある。NUCTE−III の実験では、 セパラトリクス内部の流速は Fig. 2-7 に示したモデルのように時間的に増加する ことがわかっている。電場短絡による回転説の妥当性が成り立つためには、ス クレイプオフプラズマの反磁性ドリフト速度時間発展も、同様の傾向を示す必 要がある。テータピンチで生成した FRC プラズマは、閉じ込めが悪く、徐々に 磁束が減衰し、セパラトリクス半径も小さくなる。圧力勾配も小さくなり、磁 場も時間とともに弱くなる。減衰する FRC プラズマのセパラトリクス上の反磁 性ドリフト速度を Fig. 2-9 に示す。これは、電気抵抗に関する異常係数を 10 と して磁束の減衰を再現したモデルにおける計算結果である。詳細については、 次章に述べる。Fig. 2-9 における流速の観測点はセパラトリクスであり、時間と ともに装置中心軸方向へ移動している。Fig. 2-9 から、反磁性ドリフト速度は、 徐々に減少していることがわかる。磁場や密度が減少すると、反磁性ドリフト 速度は増大するので、これらの減少よりも、セパラトリクス上の圧力勾配の減 少の影響が大きいことを示している。周辺プラズマの回転は、内部のプラズマ 回転の駆動源である。Fig. 2-9 の計算結果は、内部回転の駆動力が時間とともに なくなることを示しており、NUCTE−III の実験で得られている流速の時間発展 を説明することができない。 電場の短絡によるスクレイプオフプラズマの回転は、電子のダイナミクスの 時間スケールで決まるため、磁力線が導電壁に導かれる限り生成直後から開始 するはずであり、その速度は反磁性ドリフト速度である。この時、セパラトリ クス内部の回転は、Fig. 2-8 の結果からもわかるように、粘性力によって速やか に駆動される。しかし、実際の実験結果によるとトロイダル流速は緩やかに立 ち上がっており、電場の短絡が回転起源であるという主張にこの点でも疑問が 残る。

Fig. 2-9 セパラトリクスにおける反磁性ドリフト速度の時間発展。電気抵抗は 古典理論に異常係数 10 をかけ、磁束減衰時間が NUCTE-III の実験とほぼ一致す るようにしている。 2-2 磁束減衰に伴うプラズマの回転説 2-1 節では、FRC プラズマの回転に対する従来説について説明した。本節では 新たな説として提唱している「磁束減衰に伴うプラズマの回転説」について、 正準角運動量、運動方程式、ステルマー領域の３つの観点から説明する。ただ し、これらはあくまでも１粒子における予測であり、多数の粒子が存在するプ ラズマにおいて同様のことが言えるかは詳細な数値シミュレーションを行なわ ないと明確にできない。数値シミュレーションについては第 3 章に記述する。 2-2-1 正準角運動量からの予測 軸対称で衝突がない場合、正準角運動量とは粒子の持つ運動の恒量であり、 次の式で表される。 const. q mr p v (2.13) ここで、 m 、qはそれぞれ粒子の質量、電荷であり、 は磁束を示す。(2.13)式 において磁束の減衰を考えたとき、減衰量を とすると(2.13)式は次のように 書くことができる。 v r q m (2.14) 磁束は減衰するため 0 であり、このことから rv 0であることがわかる。 v r は粒子の持つ角運動量であり、磁束に反比例して値が増加する。したがって、

磁束の減衰に伴いトロイダル方向の速度が増加、つまり回転速度が増加するこ とが考えられる。ただし、正準角運動量はプラズマ中の粒子とバックグラウン ド粒子との衝突がないという条件下で保存するため、実際に頻繁に衝突が起こ っているプラズマ内で同じことが言えるかは分からない。 2-2-2 運動方程式からの予測 FRC プラズマは、周方向の反磁性プラズマ電流が閉じた磁力線を生成し、こ れによってプラズマ自身が閉じ込められる。プラズマ電流は電気抵抗によって 減衰するため、円筒断面を貫く磁束もこれに伴い減衰する。このとき磁束減衰 よって誘導電場Eが生じる。したがって、クーロン衝突がないときの粒子の運 動方程式は ) ( d d E B v v q t m (2.15) で表される。FRC 中の粒子は、磁場構造の特異性から Fig. 2-10 に示すように、 主に 3 つの軌道に分けられる。

(a) gyration (b) figure-8 (c) betatron

Fig. 2-10 FRC プラズマにおける粒子軌道の種類 中でもベータトロン粒子と呼ばれる軸を取り囲むように運動する粒子は、旋 回半径が大きく、運動方向は常にθ方向に同一の方向となる。これに反して gyration や figure-8 は旋回しながら、θ方向に向きを交互に変化させている。ベ ータトロン粒子の場合、電場によって常に加速されることで、θ方向の運動量 が増加する。このためプラズマの回転を引き起こすことが考えられる。

2-2-3 ステルマー領域からの予測 Fig. 2-10 に示したように FRC プラズマ中の粒子軌道は、gyration、figure-8、 betatron の３種に分けられる。それぞれの粒子がどのような軌道を描くかはステ ルマー領域計算を用いることで判断できる。ステルマー領域とは、粒子が運動 できる範囲を指しており、 q p r mH 2 (2.16) により定義される。ここで、Hは粒子の持つ全エネルギー（静電ポテンシャル がないとき）を示す。例えば、gyration 運動する粒子において(2.16)式を用いる と Fig. 2-11 のような図が描ける。この粒子はHが小さいことが特徴である。 Fig. 2-11 gyration 粒子のステルマー領域 赤線で示された r についての２つの範囲が、(2.16)式の条件を満たす範囲、つま り粒子の存在できる領域である。どちらの範囲においても p q 0となる点 を含んでおり、この点においてv 0であるということが (2.15)式からわかる。 赤線で示された範囲にそれぞれv 0である点を有しているため、その範囲内を 旋回運動する粒子であることが示唆される。したがって、Fig. 2-11 は gyration 粒 子であることは明らかである。 figure-8、betatron 運動する粒子に関しても(2.16)式を用いることで同様に判別す ることができ、それぞれの軌道におけるステルマー領域は Fig. 2-12、Fig. 2-13 のようになる。figure-8 運動する粒子は、gyration に比べてHが大きく、運動の 範囲がより広範囲になり、この結果 Fig. 2-11 では二つの範囲に分断されていた ものが、Fig. 2-12 では連結していることが理解できる。これは、ラーマー半径が 大きいために、向きの違う磁力線領域まで運動が可能となり、よって旋回方向 が逆転し結果的に 8 の字を描くことを示している。betatron 運動する粒子の大き

な特徴は、高エネルギーであることに加えて角運動量が大きいことである。こ の結果、 p の値がq の最大値（つまり磁気中性点での値）を上回り、運動可 能領域において、 v の符号が一定である。 Fig. 2-12 figure-8 粒子のステルマー領域 Fig. 2-13 betatron 粒子のステルマー領域 次に磁束が減衰したときの粒子軌道について gyration を例に考える。(2.16)式 において、磁束の減衰を考えると、粒子のステルマー領域は Fig. 2-14(a)→(c)の ように変化していく。gyration を表していたステルマー領域が磁束の減衰に伴い 次第に figure-8 のステルマー領域に変化していく様子がわかる。同様に figure-8 運動する粒子は、磁束の減衰に伴い betatron 運動へ変化する。粒子によっては、 元々gyration 運動していたものが最終的には betatron 運動へ変わるものもあるだ ろう。betatron 運動する粒子については、2-2 節で述べたように磁束減衰によっ て生じた電場によって加速される(Fig. 2-15)ので、全体的に見てもプラズマの回 転速度が増加することが予測できる。実際に磁束減衰を考慮して単一粒子の軌 道計算を行なった結果については Fig. 2-16 に示す。この図からもわかるように、 gyration 粒子のドリフト方向はイオン常磁性の方向であり、実験結果とは異なる

方向である。しかしながら、figure-8 粒子は反磁性方向であり、磁束減衰によっ て軌道変化が起こり、トロイダル回転が起こることが粒子軌道から説明しうる ことがわかる。 (a)gyration (b)figure-8 (c)betatron Fig. 2-14 磁束の減衰に伴う粒子の軌道変化

Fig. 2-15 betatron 粒子のトロイダル流速及び運動エネルギーの時間発展 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 y x gyration→figure-8 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 y x figure-8→betatron Fig. 2-16 磁束減衰による粒子の軌道変化

第 3 章 ハイブリッドシミュレーションモデル

本研究では、重水素イオンと不純物イオンの軌道計算を行ない、PIC 法を用い て超粒子の情報をもとに流体の情報に変換し、流体情報を場の計算に反映させ、 その場の中でまた粒子を動かすハイブリッドシミュレーションを行っている。 ここでは、その際に用いる計算モデルやハイブリッドシミュレーション法につ いて説明する。 3-1 計算モデル

本研究は日本大学の実験で使用されている NUCTE (Nihon University Compact Torus Experiment) -Ⅲをモデルとしている。Fig. 3-1 は NUCTE-Ⅲの装置図である。 この装置のパラメータは以下の通りである。また、計算に用いるプラズマのパ ラメータも示す。 装置半径 : 0.17[m] 外部磁場 : 0.5[T] 電子温度・密度 : 124[eV] ・ 2.5 1021[m-3] 重水素イオン温度・密度 : 124[eV] ・ 2.5 1021[m-3] 不純物イオン温度・密度 : 100[eV] ・ 2.5 1016[m-3] 本研究では、不純物イオンは 4 価の炭素を使用している。

3-2 FRC プラズマの平衡状態 3-2-1 MHD 平衡 核融合を起こすには、プラズマを高温、高密度にして容器から離してプラズ マを保持しなければならない。その質量の運動または状態変化を示す方程式に おいて時間的変化の項がゼロであるとき、その物質は平衡状態にある。磁気閉 じ込め方式による平衡状態はプラズマの圧力勾配とローレンツ力がつりあった 状態である。この節ではプラズマの平衡状態を表す式を理想 MHD 方程式から導 く。 3-2-2 理想 MHD 方程式 理想 MHD 方程式とは電気抵抗が無視できる高温プラズマの状態を表した方 程式である。理想 MHD 方程式は以下のように記述される。 0 v t (3.1) B J v v v p t (3.2) 0 p t v (3.3) 0 B v E (3.4) t B E (3.5) j B ₀ (3.6) 0 B (3.7) ここで、E,B,jはそれぞれ電場強度、磁場強度、電流密度を表している。また、 0 , , , ,v p はそれぞれプラズマの質量密度、速度、圧力、比熱比、真空透磁率 を表している。 3-2-3 Grad-Shafranov 方程式 プラズマの平衡は一般的に Grad-Shafranov 方程式を解くことにより求められ る。この節では(3.1)~(3.7)の理想 MHD 方程式から Grad-Shafranov 方程式を導出 する。 平衡状態では時間に関する項は 0 となるので t 0とし、また静的な平衡と 考えるのでu 0とする。この条件により理想 MHD 方程式は次のように整理す ることができる。 B j p (3.8) j B 0 (3.9) 0 B (3.10) ここで、(3.8)式の両辺をBで内積すると

0 B p (3.11) となる。FRC プラズマは軸対称な配位を持っているので、円筒座標系(r, ,z)を 用いる。(3.11)式から等圧面上に磁力線が沿うことから、等圧面上に乗る磁力線 の作る面を考えることができる。この面を磁気面と呼び、その磁気面を表す関 数を で表す。磁束関数 の定義は rrB_z r 0 d (3.12) で与えられる。式(3.84)の両辺を r で微分すると r r B_z 1 (3.13) 磁場Bをベクトルポテンシャルで表すと z z r r z r A A A z r r e e e e e e A B 1 (3.14)

ここで、e_r cos ,sin 、e sin ,cos であるから r r d d d d e e e e , であることに注意すると、(3.14)式は次のように表せる。 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 (3.15) モデルは円筒座標系で軸対称なので 0 である。また、磁場の θ 成分は平 衡状態のとき存在しないのでB 0である。(3.15)式を整理すると z r rA r r z A e e B 1 (3.16) となる。(3.16)式のｚ成分より、 rA r r B_z 1 (3.17) となる。(3.13)式と(3.17)式を比べると rA (3.18) であることがわかる。また、(3.16)式のｒ成分と(3.18)式より z r Br 1 (3.19) となる。 次に、(3.9)式の θ 成分をとり、(3.13)と(3.19)式で書き表すと、 r B z B J r z 0 2 2 1 1 z r r r r r (3.20) となる。一方、(3.8)式の r 成分をとると

r p B J _z となり、(3.13)式を代入することによって、 d dp r j (3.21) となる。よって、(3.20)式と(3.21)式より次のような方程式を得ることができる。 d d 1 2 0 2 2 p r z r r r r (3.22) (3.22)式は FRC における Grad-Shafranov 方程式であり、境界条件と圧力勾配を与 えることによって Grad-Shafranov 方程式を解き、磁束関数 を求めることがで きる。 3-3 イオンの初期分布 プラズマイオンを粒子的に扱い、マクロなプラズマを解析するためには多数 の粒子を同時に計算させる必要がある。正確な解析のためには全プラズマ粒子 の軌道計算を行う方法が最も良いが、現実の実験室のプラズマでは 1cm3_あたり 1010∼1014個の粒子が存在するため、それらをすべてシミュレーションすること は不可能である。そこで、粒子シミュレーションでは、取り扱う粒子は実際の 粒子の電荷質量比を一定に保ったまま、多数の粒子の電荷と質量をまとめた超 粒子として扱う。本研究では、重水素イオンと不純物イオンの 2 種類を扱って おり、それぞれ約 80 万個の粒子を配置し、軌道計算をしている。 FRC は軸対称かつ反転対称であるため、r-z 平面の 1/4 を計算することで、解析 できる。そこで、Fig. 3-2 に示すように粒子を配置し、軌道計算を行う。r 方向、 z 方向ともに 129 分割し、その一つ一つのグリッド点上に様々な速度ベクトルを 持つ粒子を配置する。速度については Fig. 3-3 のような Maxwell 分布を仮定し、 それを n 分割して与えている。また、速度は r、θ、z の３方向に関して与える必 要がある。つまり、Maxwell 分布に沿った分布を作るために、１グリット点ごと に 3 n 通りの速度を持った粒子を配置することになる。これらに分布関数 i 2 2 2 2 3 i i i 2 exp 2 T m T m n f (vr v vz ) (3.23) を用いて求めた重みを持たせる。ここで、n_i,m_i,T_iはそれぞれイオンの密度、質 量、温度を示す。密度と分布関数の関係は、 v d f n (3.24) であるから、分布関数 f に微小速度dvをかけたものは、密度の次元をもってい

る。したがって、超粒子の重みを とすると、 V f dv d (3.25) となる。 0 129 1 w

r

mir

z

0 129 1 w

r

mir

z

v

f

0 3 3 m T_i Fig. 3-2 粒子の配置方法 Fig. 3-3 Maxwell 分布

3-4 粒子軌道計算 粒子の軌道計算には運動方程式 u v E B v v m q t m d d (3.26) を用いている。ここで、 はテスト粒子で、 はバックグランド粒子である。 はテスト粒子とバックグランド粒子との衝突周波数であり、 3 2 2 0 2 2 1 4 v x m m m q q n (3.27) t t e x 2 x t d 0 _T m x 2 2 v (3.28) で与えられる。本研究のシミュレーションでは、テスト粒子は重水素イオンと 不純物イオンを用い、バックグランド粒子は電子と重水素イオン、不純物イオ ンを用いている。 また、(3.26)式より軌道を計算するには磁場B、電場E、バックグランド粒子 の流速 u が必要である。これらをどのように求めるかは、次節で説明する。 3-4-1 磁束の減衰と場の関係 FRC プラズマの磁束は時間の経過とともに減衰していく。そこで、シミュレ ーションでは、上述した平衡状態の計算により初期の磁束を求めた後、時間刻 みごとに磁束を減少させる必要がある。ここでは、磁束の減衰式と場の計算方 法について説明する。 まず、磁束の減衰式を求めよう。ファラデーの電磁誘導の法則 t B E (3.29) より、両辺を面積分すると S B S E d d t (3.30) となる。また、円筒断面積はdS 2 rdrであるため、(3.12)の定義式は S B d 2 1 (3.31) と書くことができる。(3.30)式の左辺にストークスの定理を適用し、右辺に(3.31) 式を代入すると、(3.30)式は

E E l r t d 2 2 (3.32) E r t (3.33) と書くことができる。また、電場はオームの式 e e e 1 p en j B u E (3.34) より求めることができる。ここで、u_e,n_e,p_eはそれぞれ電子の流速、密度、圧力 である。また、 は電気抵抗であり異常係数を含んでいる。実験予測される電気 抵抗は、クーロン衝突をもとに計算された古典的電気抵抗に比べて十分大きく なるため、異常係数と呼ばれる係数をかけ電気抵抗を人工的に増倍している。 磁束減衰の計算では、ローレンツ力u_e B及び圧力勾配 p については無視し、_e (3.34)式を j E (3.35) と簡略化する。この(3.35)式を(3.33)式に代入することで、 j r t (3.36) という磁束の減衰式が得られる。シミュレーションでは、この式に電流密度を 代入し、磁束の変化量を求め、これを前のステップの磁束に加えることで、磁 束の減衰を再現している。そして、磁束の時間発展が実験結果を再現するよう に、異常係数を設定している。 3-4-2 場の計算法 次に、場の計算方法について説明する。磁場、電流密度、電場といった場の 値は、磁束の減衰に変化するため、磁束と同様に時間刻みごとに計算していく 必要がある。 まず、磁場の計算方法についてである。磁場はすでに求めた(3.13)式と(3.19) 式に、減衰後の磁束を代入することで求めることができる。 次に、電流密度の計算方法についてである。(3.6)式のアンペールの法則より r B z B j r z 0 1 (3.37) となる。この式に磁場を代入することで、電流密度を求めることができる。な お、FRC における電流はトロイダル方向のみに存在するため、 j ,_r j_zについては

考えない。 最後に、電場の計算方法についてである。電場は(3.34)式のオームの式により 求める。まず、(3.34)式の右辺第一項のローレンツ力u_e Bは、後述する PIC 法 により密度n と流束_i γを集計し、 i i n γ u (3.38) により、イオン流速u を求めた後、電流の定義式 _i e i e u u j en (3.39) を用いて、電子流速u を求め、磁場との外積を計算することで得られる。右辺_e 第二項は(3.37)式から得られた電流密度を代入することで得られる。右辺第三項 の電子圧力は、電子圧力の定義式 e e e n T p (3.40) により計算する。ここで、T は電子温度である。1 ステップ前の電子温度を_e T_eoと すると、(3.40)式の電子温度は t T t T T d d eo o e e (3.41) となる。ここで tは 1 ステップの時間である。電子温度の変化量は、電子流体 に対するエネルギー方程式 ei e e e e e d d 2 3 Q T n t T n u (3.42) より、 e e e ei e 3 2 d d u T n Q t T (3.43) で与えられる。ここで、Q は熱生成項であり、エネルギー保存則 _ei e e i ie ei Q i Q u u R により得られる。熱生成項Q は後述する PIC 法により集計し求める。 _ie

3-4-3 クーロン衝突 プラズマ中の粒子は、バックグランド粒子とのクーロン衝突により速度や軌 道を変化させる。摩擦による速度の減速現象を Slowing-Down collision、軌道の 変化をピッチ角散乱という。ここでは、それぞれの現象について説明する。 3-4-3-1 slowing-down collision ここではプラズマ運動論を用いて、slowing-down collision の特性を議論する。そ こでまず、Landau’s Formula による衝突項を持つ、Vlasov 方程式

t f t f 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 t ' ' ' S_, S S S S f f C (3.44) u v g , ln 8 _min max 2 0 2 ' 2 ' k k e e LSS s s によって支配される分布関数 fS v,t は Boltzmann の H 定理を満たしていること を示そう。 0 , t fS v であるので、H 関数を t f t f t H S S S , ln , dv v v (3.45) と定義すると、次のような式を得ることができる。 t f t f 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 t (3.46) ここで、vと 、u SとS'とを交換することによって、(3.46)式と同等の式を得る ことができる。その２つの式を和と取り、両辺を２で割ると t f 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 v u v g gg g 1 ' ' t (3.47) となる。ここで、 t f m t f m S S S S SS ln , 1 , ln 1 u u v v X ' ' ' とすると、(3.47)式の被積分関数は 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 t ' 3 ' ' , , SS SS S S t f t f X g gg g 1 X u v t 2 ' 2 ' 2 3 ' ' 1 , , S _{SS SS SS} S t f t f g X g X g X u v 0 1 , , S' _{3 SS'} 2 S t f t f g X g u v (3.48) となる。したがって(3.47)と(3.48)式から 0 d d t H t (3.49) であることがわかる。すなわち、(3.45)で定義された H-関数は常に減少していく。 これを H 定理と呼ぶ。(3.48)式より、dH t /dt 0となるのは任意のS,S'とv,uに 対して t f m t f m S S S S , 1 , 1 u u v v u v q ' ' ln // (3.50) の場合のみである。(3.50)式は次式のようにも表せる。 u v u v u u v v ln , , 1 , ln 1 A t f m t f m S S S S ' ' (3.51) この条件(3.51)式は shifted-Maxwell 分布 2 2 exp , _S S B S S S T k m A t f v v V (3.52) において、T_S T_S'及びV_S V_S'の場合、満たされる。結局、H 関数(3.45)式は(3.44)

式で表さ れたク ーロン 衝突に よって 時間 とともに 減少し 、 t において v v v _' ' S S S S T T T , の shifted-Maxwell 分布となって定常値になる。 (3.44)式の衝突項は積分を含んでいるため、分布関数を決定するためには、連 立積分方程式を解かなければならないため、実際の応用には不便である。そこ で、プラズマの分布関数は時間的に Maxwell 分布に近づくことを用い、field particle の分布関数が Maxwell 分布からあまり大きくずれていない場合を想定し、 (3.44)式の簡単化を行う。まず、(3.44)式において 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 (3.53) とおく。次に、 u v v v u v u v u v u v 1 3 t (3.54) であることを利用して(3.53)を(3.44)に代入して、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 t f T m S , v v v k_B (3.55)

となり、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 (3.56) である。(3.55)式は積分を含まないので、衝突現象を記述する基礎方程式として 多く用いられる。ここで、(3.55)式において分布関数 f v,t を x t f v, v V te (3.57) とする。V はテスト粒子の速度を表し、その方向に x 軸をとる。(3.55)式の両t 辺にv をかけて、両辺を速度空間で積分すると、 _x t t V t V t d d d d d v z y x x v v v v 左辺

S j S S B S S B S S m T k m T k m e e n 2 3 j x 3 j x j 2 x x _v v v v v x ' v v v 2 8 ₀2 2 2 2 v v v 右辺 z y x j t V T m z S , , d d d _{x y z} 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 (3.58) となる。よって、次の式が得られる。 t V t V t sl d d (3.59) ここで、 S S S S x V m m m e e n 3 2 2 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 , (3.60) である。この _slは、温度T の Maxwell 分布したプラズマ中に存在するテスト粒_S

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

の場合は x 1であり、反対にV2 2k_BT_S m_S では 32 3 4 x x である。よ って、(3.60)式は 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 Q Q , (3.61) すなわち、プラズマの熱速度(v_th 2T /_S m_S )より充分に大きな速度をもつテス ト粒子の _slは、自身の速度の３乗に逆比例する。反対に熱速度より充分に小さ いテスト粒子は熱速度の３乗に逆比例した一定値に近づく。

3-4-3-2 ピッチ角散乱 粒子の速度をv 、磁場に平行な粒子の速度成分をv とし、Fig. 3-4 のようにす_// る。 を p // cos v v 1 cos p 1 (3.62) と定義する。ピッチ角 _pの変化は 2 1 p 2 o p o n 1 1 (3.63) で表す。ここで、 _pはピッチ角散乱するときの衝突周波数、 は擬似乱数を発 生させる時間刻みである。 の添え字の n は新しい速度方向を、o は古い速度方 向を示す。 はモンテカルロ法を用いて決定する。 5 . 0 0 R のとき「−」 1 5 . 0 R のとき「＋」 ここで、Rは 0 から 1 の一様乱数である。ピッチ角散乱を考慮した粒子軌道計 算をするには、この磁場に沿った座標系から円筒座標系に変換させる必要があ る。次に円筒座標系に変換する手順を示す。Fig. 3-5 のように磁場と水平方向と のなす角を 、r z平面上で磁場に垂直な速度成分をv とする。まず、Fig. 3-5 より sin cos cos sin // // v v v v v v z r (3.64) となる。次にv v 平面から見た図を Fig. 3-6 に示す。v とv のなす角を と すると、 cos sin v v sin sin v v p p θ (3.65) と表すことができる。また、Fig. 3-4 より p cos v v_// (3.66) となる。(3.64),(3.65),(3.66)式より sin cos sin v cos cos v v sin sin v v cos cos sin v sin cos v v p p p p p z θ r (3.67)

となる。(3.67)式を初期値として扱うことにより、ピッチ角散乱を考慮した粒子 軌道計算を行なうことができる。また、ピッチ角散乱の場合の衝突周波数 _pは 2 3 2 2 0 2 2 p 1 2 1 1 1 8 dx x d x x m e e n _β β α α β α β v (3.68) である。 v v p B // v p z v r v v v // v B Fig. 3-4 磁場に沿った座標系 Fig. 3-5 磁場に沿った座標系と円筒座標系の対応関係

3-4-4 PIC 法 本研究では、計算領域をセルに分割し、格子点上の物理量を求めるシミュレ ーション技法を採用している。セル内部に含まれる超粒子の位置や速度情報を 格子点に反映させ、粒子の情報を流体の情報に変換するのが目的である。ここ では、その手法の一つである PIC(Particle In Cell)法について説明する。本研究で は円筒座標系で PIC 法を用いているが、簡単のため。まずはデカルト座標で説 明する。デカルト座標系における PIC 法の概念図について Fig. 3-7 に示す。 Fig. 3-7 は２次元平面における１つのセルを表しており、グリッド点をそれぞ れ G1∼G4 と設定する。また、１つのセル全体の面積を S とし、セル内に配置 した超粒子 P によって区切られた４つの区間をそれぞれ A、B、C、D とする。 このとき、グリッド点 G1∼G4 に振り分けられる超粒子 P の情報は、自身によ って区切られた４つの面積によって次のように決められる。 2 1 S C n_G (3.69) 2 2 S D n_G (3.70) 2 3 S A n_G (3.71) 2 4 S B n_G (3.72) ここでn ∼_G1 n_G4は各グリッド点における密度であり、 は１つの超粒子が持つ 重みを示す。各グリッド点における流束 についても同様に求めることができ、 v θ v B Fig. 3-6 旋回方向の図