極限的高ベータプラズマの準平衡計算と高周波揺動解析

0

令 和 二 年 度 修 士 論 文

極限的高ベータプラズマの準平衡計算と高周波揺動解析

指導教員 高橋 俊樹 准教授

群馬大学大学院理工学府 理工学専攻

電子情報・数理教育プログラム

杉木 悠人

1

目次

第 1 章

序論

2

1.1

研究背景

. . .

2

1.2

磁場反転配位プラズマ

. . .

3

1.3

FRC プラズマの

シミュレーション方法

. . .

5

1.4

研究目的

. . .

6

第 2 章

外部磁束減衰を考慮した準平衡計算

7

2.1

研究目的

. . .

7

2.2

理論と計算方法

. . .

7

2.3

物理量の規格化

. . .

12

2.4

計算結果

. . .

13

2.5

まとめ

. . .

19

第 3 章

磁気リコネクション解析のための 3 次元フル

パーティクルシミュレーション

20

3.1

平衡計算の結果

. . .

20

3.2

シミュレーションコードに

ついて

. . .

21

3.3

分散関係の導出

. . .

22

3.4

計算の条件

. . .

29

3.5

計算結果

. . .

31

3.6

まとめ

. . .

40

第 4 章

まとめ

42

謝辞

43

参考文献

44

2

第１章 序論

1.1 研究背景

現在の日本では Fig.1-1 へ示すように、石油などの化石燃料を用いた火力発電が主な発電 方法である。しかし日本の鉱産資源の埋蔵量は乏しく、燃料の大部分は輸入に頼っているた めエネルギー自給率は他国に比べ低い水準にある。また広く知られているように化石燃料 の埋蔵量には限度があり、地球温暖化の影響となる温室効果ガスを排出することや、SOx・ NOx と呼ばれる大気汚染の原因となる物質が発生することもデメリットの一つである。 2011 年まで原子力発電は主流な発電方法の一つであったが、東日本大震災によって発生し た東京電力福島第一原子力発電所の事故の影響により原子力発電への不安は拭いきれず、 現在では発電量は大きく減少している。また、温室効果ガスの発生量や生活環境への影響が 少ない水力発電も期待できるが、ダムの設置コストの問題や建設の際に森林を破壊するこ とになるため生態系の破壊を起こしかねない。これらの影響を考えると積極的にダムを建 設することはできないといえる。 Fig 1-1. 日本の発電電力量の推移[1] このように従来の発電方式は様々な問題を抱えている。そのため世界的に新エネルギー の開発がすすめられている。新エネルギーの中でも太陽光発電や風力発電などの再生可能 エネルギーは温室効果ガスなどの影響も少ないクリーンな発電方式である。しかし Table 1-1 へ示すように発電量が少ないことや天候などに大きく影響を受けることが問題点として

3 挙げられる。そのためまだ安定した電力供給は望めない現状である。 そこで、新エネルギーの中でも核融合発電に注目する。原子力発電が「核分裂」反応によ り発生する核エネルギーを用いる発電方式である一方、「核融合」反応を利用する発電方式 である。原子力発電では反応によって発生する中性子や放射線が問題となるが、核融合反応 では比較的中性子の発生が少なく安全であるといえる。この核融合反応は太陽でも行われ ており、太陽は核融合反応によって熱と光を生み出している。しかしこの核融合反応を地上 で行うには、高密度・高温なプラズマ状態が必要であるという問題点もある。プラズマ状態 とは原子が陽イオンと電子に電離した状態である。核融合により安定した電力を得るため には、高密度で高温のプラズマを核融合反応が起こるまで定常的に維持する必要がある。当 研究室では核融合発電を目指した核融合プラズマのシミュレーション研究を行っている。 本研究では後述する FRC プラズマというプラズマを対象とし、FRC プラズマの諸特性を解 明するための研究を行う。

1.2 磁場反転配位プラズマ(FRC プラズマ)

当研究の対象となる FRC プラズマについて示す[3]。FRC とは磁場反転配位(Field–Reversed Configuration )の略である。 FRC プラズマはその名前の通り、プラズマ内部で磁場が反転す ることが特徴的である。また、開いた磁場の中に磁力線が閉じている領域があり、その境界 を separatrix という。磁場が反転するためセパラトリクス内部中央面上に O-point 、Z 軸上 に X-point という磁気中性点と呼ばれる点がある。FRC プラズマの構造を Fig. 1-2 へ示す。 O-point は平面上で見ると点であるが、3 次元の円筒座標系を考えると点の集合となり、 Field-null circle という装置軸を囲む円とみなせる。 この FRC プラズマの特徴として閉じ込めの指標となる 𝛽 値が、他の閉じ込め方式よりも 高くなるという特徴がある。𝛽 値が高いということは、同じ外部磁気圧を与えた場合に高 いプラズマ圧を得ることができる。また同じプラズマ圧を得ようとした場合に外部磁場が 小さくて良いため外部コイルへ流す電流が少なくて良い。そのため経済的で高効率である、 Table 1-1. 原子力発電所 1 基の年間発電電力量との比較[2] 発電種類 原子力発電 太陽光発電（住宅用） 風力発電 1 基あたりの設備容量 100 万 kW 3.5kW 1,000kW 設備利用率 80% 12% 20% 必要な基数 1 基 約 190 万基（軒） 約 4,000 基

4

という利点がある。𝛽 値は(1-1) 式のように表される。 𝛽 =プラズマ圧

外部磁気圧 (1-1)

FRC プラズマの装置の例として Fig.1-3 に示す日本大学 NUCTE-Ⅲ(Nihon University Compact Torus Experiment)がある。これは逆磁場シータピンチ(FRTP)法[3]_{でプラズマを生}

成する装置である。 Fig. 1-2. FRC プラズマの構造[4] Fig. 1-3. NUCTE-Ⅲの装置概形 [5] この FRC プラズマの問題点として配位維持時間の短さが挙げられる。これは異常抵抗な どにより磁束の減衰が大きいことや、回転不安定性などによる配位の崩壊など原因として 考えられる。配位時間を伸長させるためにビーム入射やプラズマ同士を合体させる実験や シミュレーションが行われている。当研究室の松崎はプラズマの移送合体をシミュレーシ ョンした。結果を一部示す。

5 Fig. 1-4. MHD シミュレーションによる FRC の衝突・合体 [6]

1.3 FRC プラズマのシミュレーション手法

先に示した松崎の結果は MHD シミュレーションという手法で行っている。FRC プラズマ のシミュレーション方法は以下に示すように三種類ある。 フルパーティクル イオンと電子について、運動方程式から軌道を計算し、電場・磁場 の変化を計算する手法。計算コストは大きい。 ハイブリッド イオンは粒子として計算するが、電子は流体と近似して計算する。 フルパーティクルシミュレーションより計算コストは少ない。 MHD プラズマ全体を電磁流体と近似し、巨視的な物理量の変化を計算す る手法。計算コストは二つよりも少ない。 FRC プラズマの場合では、先に述べたように磁場が 0 となる点が存在する。この点におい て粒子の旋回半径は大きくなり、運動論的な動きをする。そのため FRC プラズマではフル パーティクルシミュレーションかハイブリッドシミュレーションが好ましいとされている。

6

1.4 研究目的

本研究では FRC プラズマを基にした計算を行っている。第 2 章においては装置外部のコ イルの磁場が時間変化する場合の内部のプラズマの変化について述べる。この研究では日 本大学での実験を基に、従来の計算手法では計算できなかった外部磁場の変化を含んだ計 算を行い、外部磁場が定常である場合と比較してどのように変化が異なるかを確認するこ とを目的とした。また、第 3 章ではフルパーティクルシミュレーションによる高密度プラ ズマの揺動解析を行う。ここでは磁気リコネクション解析を目指したフルパーティクルシ ミュレーションコードを開発し、そのコードの妥当性を検証するとともに高密度プラズマ における波動伝播を確認することを目的とする。

7

第 2 章 外部磁束減衰を考慮した準平衡計算

2.1 研究目的

現在、広くシミュレーション研究がおこなわれており、フルパーティクル・ハイブリッド などの粒子シミュレーションやプラズマ全体を電磁流体と近似する MHD という手法がある ことを序論にて述べた。これらのシミュレーション方法では、プラズマの変化を計算する手 法であり、装置の外部コイルから与えられる磁場は一定である[7],[8]_{と考えて計算を行ってい} た。外部磁場が一定であればセパラトリクス半径は一様に減衰すると考えられる。しかし板 垣が行った NUCTE - III での実験[7]_{は、外部磁場が時間変化することにより、セパラトリク} ス半径が一定に近い値となる結果が得られた。外部磁場の時間変化とセパラトリクス半径 の時間変化を示す。 Fig. 2-1. 外部磁場(上)とセパラトリクス半径(下)の時間変化[9] 当研究では外部磁束の時間変化を考慮した計算手法を開発した。この計算手法では外部 磁束の時間変化が計算できる。この新しい計算手法による結果と、従来の計算手法から求め られたセパラトリクス半径の時間変化などを比較し、外部磁束減衰の影響を調べることを 目的とした。

8

2.2

理論と計算方法

当研究では日本大学にある装置 NUCTE-Ⅲ(Nihon University Compact Torus Experiment)を 装置モデルとした。装置概形と計算に用いたパラメータをそれぞれ Fig.2-2、Table 2-1 へ示 す。セパラトリクス圧力と装置壁の磁束、圧力は後述する圧力関数のパラメータである。 Fig. 2-2. NUCTE-Ⅲの装置概形[10] Table 2-1. NUCTE-Ⅲの装置パラメータと計算パラメータ 装置半径 0.17 [m] 𝑟w 装置長 2.00 [m] 2 𝑍M 外部磁場 0.50 [T] 𝐵ex イオン温度 100 [eV] 𝑇i 電子温度 50.0 [eV] 𝑇e セパラトリクス圧力 4.88 [ - ] 𝑝s 装置壁での圧力 1.00 × 10−7 _{[ - ] 𝑝} w 装置壁での磁束 -1.00 [ - ] 𝜓w 従来のシミュレーションなどでは式 (2-1) に示す Grad-Shafranov 方程式を解くことによ って平衡状態の磁束を求める。 (2-1) ここで 𝜓 は磁束関数であり、𝜇0 は真空の透磁率である。d𝑝/d𝜓 は圧力関数の微分である。 この方程式では圧力を磁束の関数とみなしており、圧力関数を磁束の多項式とみなすと以 下のように示すことができる。また 𝜓 < 0 となるセパラリクス外部では指数関数的に減少 していくものと考える[11]_{。そのような仮定の下で圧力関数は式 (2-2)のように表現できる。}

9 (2-2) ここで以下の条件を与えて係数を計算する。 ・セパラトリクス上で ・装置壁において ・セパラトリクス上で関数が連続 ・セパラトリクス上で一階微分、二階微分が連続 これらの条件から係数を決定すると式 (2-3) のようになる。ただし係数 は自由度を持 たせた係数である。 (2-3) また、これを磁束について微分すると d𝑝/d𝜓 が得られる。 (2-4) この関数から求めた圧力 と電流 の中央面分布を Fig. 2-3 に示す。 Fig. 2-3. 圧力と電流密度の中央面分布 ここで、磁束関数について言及する。磁束 とは式 (2-5) のような関係がある。今後計 算では磁束関数を用いて計算を行う。

10 (2-5) 慣習的に磁束関数を「磁束」と呼ぶため、断りがない場合、磁束関数のことを単に「磁束」 と呼ぶことにする。 Grad-Shafranov 方程式を解く際に、装置壁における磁束の値を境界条件として与える。境 界条件の分布を Fig. 2-4 に示す。 Fig. 2-4. 装置壁における磁束の値 示したパラメータと境界条件を用いて計算した平衡状態の磁束の 𝑟 − 𝑧 平面分布を Fig. 2-5 に示す。 Fig. 2-5. 平衡状態の磁束の平面分布 次に、従来の磁束の時間変化の計算方法について述べる。磁束の時間変化はファラデーの 電磁誘導の法則から

11 (2-6) と計算される。ここで A は異常係数という係数であり、古典的抵抗率 𝜂 に乗ずることで異 常抵抗を表している。プラズマの実験では異常に抵抗値が大きくなることが確認されてお り、クーロン衝突による抵抗の他に要因があると考えられている。シミュレーションでは係 数を用いて異常な抵抗値を表している。複数の異常係数から磁束の時間変化を計算した結 果を Fig. 2-6 に示す。単調な時間変化をしていることが確認できる。 Fig. 2-6. 磁束の時間変化と異常係数の関係 ここまで、従来の平衡計算と磁束の時間変化の計算方法について述べた。次に準平衡計算 について説明する。従来のシミュレーション方法では外部磁束の変化を考慮した計算を行 うことができないため、当研究では各ステップで近似的に平衡が成り立っていると仮定し て平衡計算を行う「準平衡計算」という手法を用いて計算を行った。各ステップの平衡計算 には圧力関数のパラメータと境界条件を時間変化させて利用する。以下に時間変化の計算 方法を示す。 (1) 圧力関数のパラメータ 圧力関数のパラメータにはセパラトリクスでの圧力 𝑝s 、装置壁での圧力 𝑝w 、装置壁での 磁束 𝜓w がある。圧力のパラメータはエネルギー閉じ込め時間 𝜏E を時定数として指数関数 的に減衰すると考え、𝜓w は外部磁束と同じように減衰すると考える。外部磁束減衰の時定 数を 𝜏F として、各パラメータの初期値を 𝑝so , 𝑝wo , 𝜓wo とすると

12 (2-7) (2-8) (2-9) と表すことができる。ここでエネルギー閉じ込め時間とは、対流損失や放射損失などによっ て決まる値[12]_{であり、エネルギー閉じ込め時間の値によって磁束の時間変化がどのように} 変化するかを調べる。 (2) 境界条件 平衡状態での境界条件は Fig. 2-4 に示した分布になる。外部磁束の減衰により境界条件も 変化する。外部磁束減衰の時定数を としているので、平衡状態の境界条件を とす ると、時刻 𝑡 における境界条件 は (2-10) のように表すことができる。 外部磁束の分布を確認するために、外部磁束とプラズマの磁束を Grad-Shafranov 方程式 から計算した。まず、外部磁束の計算は(2-1)の右辺を 0 とすることにより計算できる。 (2-11) 次にプラズマの磁束は、壁での境界条件を 0 として計算を行う。準平衡計算における境 界条件の変化は外部磁束の変化に対応し、圧力関数のパラメータの変化はプラズマの磁束 の変化に対応する。外部磁束とプラズマの磁束の 平面分布を Fig. 2-6 に示す。 Fig. 2-7. 外部磁束(左)とプラズマの磁束(右)の平面分布

2.3 物理量の規格化

数値計算を行う際に、物理量を規格化して無次元量へ変換して計算を行う。Table 2-1 へ 規格化定数を示す。 は壁での磁束の値、 は真空の透磁率であり、 と はそれぞ れ装置半径と装置長である。規格化時間は を用いて表す。

13 Table 2-1. 物理量の規格化定数 磁束 電流 磁場 時間 圧力 装置長 装置半径

2.4

計算結果

準平衡計算を行った結果を示す。エネルギー閉じ込め時間 と外部磁場の時定数 の 両方の値を変化させ、それらの時定数によりどのようにセパラトリクス半径が変化するか を計算して求めた。この計算方法は平衡状態から まで準平衡計算してセパラトリクス 半径の増加率を計算した。Figure 2-8 にセパラトリクス半径増加率の結果を示す。増加率は 平衡状態の のセパラトリクス半径との比から計算している。図中の破線は →∞とし た時のセパラトリクス半径の増加率である。このとき、 → → というパラメー タが得られる。理論上はこのパラメータを用いて計算を行ったときのセパラトリクス半径 増加率となると考えられる。 この結果から の値が大きいとセパラトリクス半径増加率は小さくなることが確認で きる。これは外部磁束の変化がゆるやかな場合、単純な減衰に近い状態となるためだと考え られる。また、 が大きいとき、セパラトリクス半径増加率の変化が小さくなっているこ とが確認できる。これはプラズマ内部の変化が小さく、セパラトリクス半径増加率が飽和す るためだと考えられる。 時定数 と の関係が得られた。実験の外部磁場の時間変化からの時定数を計算する。 先ほど示した Fig. 2-1 から外部磁束はおよそ で半減していることが確認できる。ここ から時定数を計算すると という値が得られる。外部磁束はこの時定数を用い て時間変化を計算する。この時定数を用いてセパラトリクス半径が一定近い条件をもとめ ると Fig. 2-9 に示すように とした場合においてセパラトリクス半径は一定に 近いということが確認できた。 これらの時定数を用いて準平衡計算を行った結果を Fig. 2-10 に示す。この結果ではセパ ラトリクス半径は一定となることが確認できたが、計算途中において急激な値の変化が発 生していることが確認できる。この値の変化の前後において平面分布を確認した。この結果 を Fig. 2-11 に示す。

14

Fig. 2-8. エネルギー閉じ込め時間によるセパラトリクス半径増加率の変化

15 Fig. 2-10. 磁束の最大値とセパラトリクス半径の時間変化 Fig. 2-11. 不連続な値の変化前後における磁束の平面分布 ここで準平衡計算の数値計算手法について考察する。このような平衡・準平衡の計算は反 復法で収束計算を行っている。従来は SOR 法という手法で計算を行っていた。簡単に説明 するとガウス=ザイデル法をもとにした計算手法であり、加速パラメータ を

16 の範囲で設定することで収束を加速させる手法である。 としたとき、ガウス=ザイデ ル法と一致する。この加速係数による歪みが磁束の収束に影響していると考え、ガウス=ザ イデル法を用いて準平衡計算を行った。この計算による磁束の変化を Fig. 2-11 に示す。約 2 倍程度の時間まで計算できているが、やはりこの場合も不連続な変化をしていることが確 認できた。 FRC プラズマではセパラトリクス半径が十分な大きさでないと磁場が反転しないため、 準平衡計算ではセパラトリクス半径の値によって近似的な平衡が計算できないことが確認 できた。 Fig. 2-12. ガウス＝ザイデル法による準平衡計算の結果

2.5 まとめ

本研究では、外部磁束減衰の効果を取り入れた準平衡計算という計算手法を開発し計 算を行った。近似的に平衡が成り立っているとした場合、セパラトリクスにおける圧力の値 によっては平衡が計算できなかった。また、エネルギー閉じ込め時間 𝜏E と外部磁束の時定 数 𝜏F との関係によりセパラトリクス半径がどのように時間変化するかを確認した。 今回はセパラトリクス半径が一定という条件のもと、エネルギー閉じ込め時間 𝜏E を計算 した。本研究では実験におけるエネルギー閉じ込め時間を見積もることが可能になる。 また、この計算手法は MHD 的な平衡を計算する手法であり、時間変化する場を平衡状態と 近似して計算を行うには準平衡計算という手法では限界がある。磁場が 0 となる領域が存 在する FRC プラズマでは粒子の効果が大きいため、粒子シミュレーションにおいて外部磁 束減衰の効果を取り入れることも重要だと考えられる。

17

第 3 章 3 次元フルパーティクルシミュレーション

での磁力線方向の波動解析

3.1 研究背景

序論にも示したように、プラズマ同士を合体させて配位時間を伸長させようという実験、 シミュレーションが行われている。当研究室の松崎が行ったシミュレーションは MHD とい う手法で巨視的な変化を計算していた。ここでプラズマ同士が合体する瞬間を Fig. 3-1 へ示 す。このプラズマ同士が衝突する部分において、閉じた磁力線が繋ぎ変わる「磁気リコネク ション」いう、現象が発生する[13][14]_{。磁気リコネクションの概略図を Fig. 3-2} に示す。X-point 付近では磁場が小さいため、電子の旋回半径が大きくなり運動論的な扱いをしなけれ ばならない。また、磁気リコネクションが発生する際には電子の動きが重要であると考えら れているため、電子の軌道計算ができるフルパーティクルシミュレーションに注目した。さ らに従来のコードは円筒座標系の 平面を計算領域とする 2 次元シミュレーションを 行っていたが、磁気リコネクションのような 3 次元的な現象の解析には適さない。このよ うな理由から電子の軌道を計算できる 3 次元フルパーティクルシミュレーションコードを 開発した。 このシミュレーションコードの計算の妥当性を確認するために簡単なモデルで、密度・磁 場を一様に設定したモデルを用いて、揺動を加えた時の波動伝播を計算するベンチマーク テストを行う。また、ベンチマークテストを行うと同時に分散関係の密度依存性を確認する ことを、本研究の目的とした。 Fig. 3-1. MHD シミュレーションによる 2 つの FRC の衝突・合体[6]

18

Fig. 3-2. 磁気リコネクションの概略図[13]

3.2 シミュレーションコードについて

磁気リコネクションは 3 次元的な現象であるため、新たに 3 次元のフルパーティクルシ ミュレーションコード CEDAR (Calculation of Electron motion Distributed around X-point for Analsis of Reconnection of magnetic field : 磁気リコネクション解析のための X-point 近傍に 分布する電子の挙動計算) を作成した。本研究においてはベンチマークテストを行うために 特別に初期分布や粒子配置・メッシュ数などを設定した。こちらについてはベンチマークテ ストにおける条件と共に説明する。ここではまず場の計算と粒子の軌道計算について述べ る。 ・場の計算 本研究で開発したコードは 3 次元のデカルト座標系のコードを用いている。円筒座標系 の 𝑟 − 𝑧 平面で平衡計算を行い、平衡状態の磁束 と磁場 を読み込み 3 次元座標に 変換することで初期の磁場を設定する。場の規格化時間 は電子のサイクロトロン周波数 の逆数とし、 であり、計算の時間刻みは としている。メッシュ数 は 3 方向に 65 メッシュとしている。また、各方向に計算領域端は反対側と同じであるとす る周期境界条件を導入している。 ・軌道計算 軌道計算は場の刻みの 1/10 であり、10 回軌道計算を行い、密度 と流束 を 集計し、平均をとる。軌道計算も同様に周期境界として考えている。今回の行う計算の流れ を Fig. 3-3 に示す。

19 Fig. 3-3. ベンチマークテストにおける計算の流れ 次に方程式の規格化について述べる。今回計算する物理量を Table 3-1 に示す。 Table 3-1. 計算に用いる物理量 磁場 電場 イオン流速 電子流速 イオン密度 電子密度 電流 この方程式を無次元量に置き換えて計算を行う。Table 3-1 に示した物理量を無次元量と する規格化定数を設定している。規格化定数について Table 3-2 に列挙する。本研究では密 度の値を変更して計算を行っているため、密度の規格化定数は異なる値となる。 Table 3-2. 規格化の式と規格化定数

20

3.3 分散関係の導出

ここで波動伝播を計算するため、波数 と周波数 の関係である分散関係を導出した。 本研究で得た分散関係を式 に示す。 (3-1) ここで、 は電子とイオンのサイクロトロン周波数であり、 は電子とイオンの プラズマ周波数である。それぞれの計算式を式 (3-2) に示す。 (3-2) という式で求められる。この分散関係について述べる。 まず 本研究で は対象とする プラズマ は一様密度・ 一様磁 場 とし、冷たい プラズ マ とする。 方向にはそれぞれ一様であるとする。 次に基礎方程式について述べる。式 (3-3) から式 (3-9) に基礎方程式を示す。 ・連続の式 (3-3) (3-4) ・運動方程式 ∇ (3-5) ∇ (3-6) ・マクスウェル方程式 (3-7) (3-8) ・電流の式 (3-9) ここで は光速、 は真空の透磁率である。 は素電荷であり、 はイオンの電荷数で ある。これらの物理量を、平衡を表す 0 次の量(添え字 : 0) と摂動を表す 1 次の量(添え字 :

21 1)に分けて考える。平衡量における条件を以下に示す。 𝑛e0= 𝑛i0= (constant) (一様密度), (一様磁場), この条件を用いて式を線形化する。線形化の条件として摂動項同士の積は高次の項とみ なし積は 0 とする。成分について書き下すと式 (3-10) から式 (3-15) のように表される。 ・連続の式 (3-10) ・運動方程式 (3-11) (3-12) ・マクスウェル方程式 (3-13) (3-14) ・電流の式 (3-15) ここで、波の空間・時間変化を表す式は として表すことができ る。よって各物理量の空間微分・時間微分は式 (3-16) のように計算できる。 (3-16) この関係を用いて線形化した方程式を表すと式 (3-17) のような関係式で表すことができ る。

22 (3-17) −i𝜔𝑚i𝑢i1𝑧 = 𝑍𝑒𝐸1𝑧 これを一つの変数で表示することで波動の伝播が計算できる。まず 成分の物理量につい て変形する。流速と電場についてまとめると式 (3-18) のように変形できる。 (3-18) ここで、それぞれの物理量は 1 次の電子密度 で表記できるので、これを式 (3-19) に示 す。 (3-19) また、この場合の周波数 について変形すると (3-20) の関係が得られる。 ∴ (3-20) 次に 成分について式変形を行う。まず電場と磁場の関係式から

23 (3-21) の関係が得られる。これを用いて磁場の関係式となるように変形する。まず について 運動方程式を変形して磁場の関係式とする。これを式 (3-22)に示す。 ∴ (3-22) 同様にほかの成分についても変形した結果を式 (3-23) から式 (3-25) に示す。 (3-23) (3-24) (3-25) そして と の関係を求める。式 (3-14) 式について変形すると という関係が得られるので、この関係式から と の連立方程式が導ける。 この連立方程式の係数は両式で同じ係数となっているので、

24 と表すことができる。ここで係数は ≡ ( ) ≡ ( ) と置いた。この方程式の解は となるので のとき、 のとき、 となる。これを変形すると式 (3-1) の分散関係の式が表すことができる。 (3-1) 式中の によって物理量の表記が異なる。物理量はそれぞれ で表すことができる。 まず なので式(3-21)から(3-25)にこの関係を代入して式(3-26)から(3-32)のように表 すことができる。 (3-26) (3-27) (3-28) (3-29) (3-30) (3-31)

25 (3-32) 式(3-26)から式(3-32) に示した物理量は に揺動を与えることで揺動を求めることがで きる。今回の研究では で計算される揺動は 0 とし、 方向の揺動のみ考えることと した。 分散関係の式を改めて示すと (3-1) という式になる。プラズマ周波数は密度の値により異なるため分散関係の曲線の形も変化 する。本研究では最大密度 を と FRC プラズマを想定した密度に設定し、 計算を行う。ここで、電場の偏光の方向により R 波と L 波という二種類の分散関係がある。 この分散関係をグラフに描くと Fig. 3-4. のように表される[13]_{。式 (3-1) の ± において、+と} した場合は R 波、− とした場合は L 波となる。 Fig. 3-4. R 波と L 波の分散関係のグラフ[15] ここで となる解 は式 (3-33) のように表される。 (3-33) 本研究では密度の値を と設定している。この密度の値の R

26 波と L 波の分散関係の曲線を Fig. 3-5 に示す。縦軸と横軸の値は電子サイクロトロン周波 数 で規格化した値を用いている。 Fig. 3-5. 密度の値による分散関係式のグラフ R 波の曲線は電子サイクロトロン周波数を漸近線に持つことが確認できる。一方、L 波の 曲線においては、さらに低い周波数領域において解を持つ。イオンサイクロトロン周波数を 漸近線にもつ曲線となる。例として の場合の低周波部分の L 波の分散関係 のグラフを示す。

27 Fig. 3-6. イオンサイクロトロン周波数付近における分散関係式のグラフ 計算においては は整数 という条件のもとに を決定する。以上の分散関係 の結果から分散関係は 4 通りの解を持つことが確認できる。

3.4 計算の条件

(i) 対象とするプラズマ 今回の研究では粒子を等間隔に配置する。計算では 𝑥, 𝑦 方向に一様であると想定してい るため、𝑥, 𝑦 方向には 3 格子点と設定し、セル数を減らして計算を行う。これによりセル内 の粒子数は増やして計算することができる。あるセルの格子点 𝑥𝑛 から𝑥𝑛+1 までに粒子を配 置したときの粒子位置 𝑋 は式 (3-34) から計算できる。 (3-34) ここで ℎ はセル幅、𝑀 はセル内の一方向に配置する粒子の数である。本研究では粒子を 3 方 向に 16 個配置した。1 セル内の粒子の数は163 _{個である。また、粒子は、粒子複数個分の} 重みをもつ超粒子として扱っている。超粒子が持つ重みは密度により計算できる。メッシュ 内の密度を としてセル内の体積を とすると式 (3-35)によって重みが計算できる。 (3-35) また今回は 𝑥, 𝑦 方向における波動伝播について考える。そのため電子密度 𝑛e1 によって与 えられる波動については 0 として考える。

28 (ii) 分散関係から求められる 𝜔, 𝑘 の値 本研究では上述した 3 通りの密度の値を用いて計算を行う。周期境界条件を用いているた め、波数 𝑘 は 𝜋 の整数倍となる値を設定する。今回はそれぞれの密度において とな る点において解を求め計算を行う。R 波の における解を用いて計算を行っている。 この理由として、他の領域における解について考えると、R 波と L 波両者とも の領 域では、密度の値によっては解が高周波となるため実験などでは観測できない波動となる 可能性がある。また、L 波の 周辺の解は波が観測できるまでの計算 step 数が膨大となる ため、電子サイクロトロン周波数に近い値となる解をもつ領域に注目した。このときの解の 値とプラズマ周波数とサイクロトロン周波数について Table3-3 と 3-4 に値を示す。角周波 数と周波数についてそれぞれ値を示してある。 Table 3-3. 密度の値による角周波数の値 (単位: ) Table 3-4. 密度の値による周波数の値(単位: ) における条件を用いるため、式(3-26)から(3-32) に示した初期分布は、式(3-36)から(3-42) のように表される。 (3-36) (3-37) (3-38)

29 (3-39) (3-40) (3-41) (3-42) 例として 1.00 × 1020 _{における分布を示す。規格化値で分布を示している。} Fig. 3-7. 𝑛0= 1.00 × 1020 のときの物理量の初期分布

3.5 計算結果

1. 定点での時間変化と周波数解析

3 種類の条件において原点における時間変化を計算した。10 万 step の計算を行い、磁場 𝐵𝑥 と電子流速 𝑢e𝑥 の時間変化を確認した。この結果を Fig. 3-7 と Fig. 3-8 に示す。時間変化

30

とする振動が確認できるかを、時間変化の結果をフーリエ変換して周波数解析を行った。

31 Fig. 3-9. 磁場と電子流速の原点での時間変化 下の図は までを拡大した図 ここでフーリエ変換をする際に、有限区間以外で値が 0 となるような窓関数を用いて計算 を行っている。窓関数を用いることで、ピークによるサイドロープの影響を減らすことがで きる。今回はハニング窓を用いて計算を行った。磁場 のフーリエ変換の結果を Fig. 3-9 に、電子流速 のフーリエ変換を Fig. 3-10 に示す。定常成分は除去しており、ピークの 値は最大値で規格化している。磁場と電子流速の二つのフーリエ変換の結果を比較すると、 おおよそピークとして確認できた周波数の値はほぼ同じ値であった。低周波部分のピーク は分散関係の 𝜔 によるものだと考えられる。また、どの周波数成分が大きく表れているか は異なっており、電子流速の場合、電子の運動の影響を直に受けるために振動が大きく表れ ているためだと考えられる。プラズマ周波数などよりも高い周波数が確認できる。これは非 線形な効果による高周波波動であると確認できる。

32

Fig. 3-10. 磁場 𝐵𝑥 のフーリエ変換の結果

33 参考として窓関数を入れていない場合の物理量の変化についても確認した。この結果を Fig. 3-11 と Fig. 3-12 に示す。低周波部分のピークの影響が大きく、小さい周波数成分が結 果として観測しにくいことが確認できる。窓関数を導入する場合、ピークの幅が太くなり、 特定の周波数の値を検出しにくくなる傾向がある。 Fig. 3-12. 磁場 𝐵𝑥 のフーリエ変換の結果(窓関数なし) Fig. 3-13. 電子流速 𝑢e𝑥 のフーリエ変換の結果(窓関数なし)

34

2. 位相速度の計算

フーリエ変換の結果から、初期値として与えた 𝜔 の値に近い成分が確認できた。次に線 形波動の位相速度から求めた 𝜔 と理論値との比較を行った。この計算は中央面分布の時間 変化を確認することで計算できる。位相の時間変化を Fig. 3-13 に示す。この結果は 1 万 step までの結果を用いて中央面分布の変化を描画した。ただし、密度 の条件においては変化が小さいため 2 万 step までの結 果を用いている。この結果においても各 step の最大値で規格化しており、位相の変化を確 認しやすいようにしている。結果を確認すると、線形な位相に高周波の振動が混ざっている ことが確認できる。 Fig. 3-14. 中央面分布の時間変化

35 この位相の変化から同位相の点から時間変化を求めると Fig. 3-14 のような結果が得られ る。ここで、位相の変化はこの位相の時間変化の結果から位相速度を求めることができる。 この結果と理論値の比較を Table 3-5 に示す。これらの結果から密度が高いほど位相速度の 値が理論値から離れることが確認できる。また、高周波の振動は原点における時間変化でも 確認された高周波の振動であると考えられる。 Fig. 3-15. 同位相の点の時間変化

36 これらの値を用いて密度と誤差の相関を調べた。計算結果と理論値から相対誤差を求め る。相対誤差 𝜀r は式 (3-43)から計算できる。 𝜀r= | 理論値 − 計算結果 理論値 | (3-43) 密度と相対誤差の関係を Fig. 3-15 に示す。この結果から真空に近い状態であると位相速 度の誤差は少ない線形波動となることを意味している。 Table 3-5. 理論値と計算結果の値の比較 理 論 値 𝜔 [rad/ns] 65.3 29.4 4.65 𝜔/𝑘 [m/s] 2.08 × 108 9.37 × 107 1.48 × 107 計 算 値 𝜔 [rad/ns] 61.4 19.2 2.44 𝜔/𝑘 [m/s] 1.95 × 108 6.12 × 107 7.75 × 106 Fig. 3-16. 密度と位相速度の相対誤差の相関

37 密度が の場合ではおよそ理論値に近い位相速度が計算できた。一方で現在は のときの波数 を設定していたが、他の の値を用いて結果を確認した。波数の 値 から求められる を用いて 1 万 step 計算を行い、位相速度を求めた。Figure 3-17 に分散関係の曲線と𝑘̅ = 𝜋, 2𝜋 における解を示す。上で示した図と横軸は変えてある。 𝜔 の解を小さい順に 𝜔̅1, ⋯ , 𝜔̅4 とし、この曲線上に図示した。分散関係の 𝜔, 𝑘 の値と計算結 果を Table 3-6 にまとめて示す。位相の時間変化は Fig. 3-17 に示した。この結果から位相 速度が大きいほど高周波の振動の影響が現れることが確認できる。 Fig. 3-17. 分散関係のグラフと𝑘+= 𝜋, 2𝜋 における解 Table 3-6. 理論値と計算結果の値の比較 理論値 計算値

38 Fig. 3-18. 𝜔̅1, ⋯ , 𝜔̅4 を用いたときの位相の時間変化

3.6 まとめ

本研究では線形波動の伝播の計算を行った。低密度の場合では位相速度は理論値に近 い値であることが確認できたが、高密度の場合ではおよそ理論値の半分の値となっている ことが計算できた。また、位相の変化から確認できる高周波の振動について考察すると、セ ル内の総粒子数が少ないために発生した振動である可能性がある。結果の比較として、粒子 数を現在は 163 _{個配置して計算を行っているが、各方向に倍の粒子を配置し 32}3 _個の粒子 を用いて同様の計算を行った。この結果を Fig. 3-16 に示す。密度の値は 1.00 × 1020 _[m−3_] として計算を行っている。この結果はおよそ 3 万 step= 34.0 × 10−11 _{[s] までの結果を確認} している。結果を比較するとほとんど同じ値で変化していることが確認できた。これは粒子 数が多くなるとそれに伴い一つの粒子数の重みが減るためであると考えられる。 粒子数、粒子の配置方法などに起因する振動ではなく、密度の値に従って高周波の振動と なっていることから、この振動は非線形な効果による高周波な波動であると考えられる。高 密度のプラズマでは線形理論よりも複雑な現象が確認され、密度の値に従って非線形な効 果が大きく表れたと考えられる。

39

40

第 4 章 まとめ

2 章において FRC プラズマの準平衡計算について述べた。準平衡計算は近似的に平衡状 態を仮定していたため、粒子の効果などの影響が入っておらず計算が破綻してしまうこと が確認できた。本来は温度などのパラメータで決まるエネルギー閉じ込め時間 𝜏E を完全な 時定数として考えていたためこちらも温度の時間変化から 𝜏E を計算し、その値を用いて圧 力などを計算するシミュレーションなどが必要であると考えられる。この外部磁束減衰の 問題はプラズマ内部の高周波の振動などが影響している可能性も考えられる。このため、フ ルパーティクルシミュレーションにおいて、電子の動きなどを計算することでさらに詳し い解明が期待できる。 3 章においては FRC プラズマを想定したフルパーティクルシミュレーションを行った。 分散関係を求めて、波動の伝播を計算したが高周波な振動が発生することが確認できた。こ の振動は非線形な効果に起因する高周波の波動であると考えられる。高密度のプラズマで は線形理論よりも複雑な現象が大きく確認される。計算の規模を拡大することにより地切 リコネクション解析などのメカニズム解明も期待できる。

41

謝 辞

本研究を進めるにあたり、日頃よりご指導いただいた髙橋俊樹准教授に厚く御礼申し上 げます。日頃より熱心で手厚い指導を賜り、大変感謝しております。 また、お忙しいなか主査・副査を引き受けていただいた高田和正教授と橋本誠司教授に深 く感謝したします・ 本研究室の皆様には日頃から大変お世話になりなした。浦野さんは同じ FRC プラズマの 研究をされており、さまざまな面から支援を頂戴しました。この場を借りて御礼申し上げま す。 同期の茂木と角掛は 3 年間同じ研究室に所属していた仲です。日頃より少々やかましく も話をしていた日々が思い出です。研究室の後輩たちも含めたくさんの刺激を得ることが できました。ありがとうございました。

42

参考文献

[1] 経済産業省エネルギー庁エネルギー白書 HP, http://www.enecho.meti.go.jp/about/whitepaper/2018html/2-1-4. [2] 経済産業省エネルギー庁日本のエネルギー2008, http://www.enecho.meti.go.jp/category/others/tyousakouhou/kyouiku/panhu/pdf/ energy2008.

[3] M. Tuszewski, Nucl. Fusion 28, 2033 (1988).

[4] Loren C. Steinhauer, Phys. Plasmas 18, 070501 (2011).

[5] 日本大学, 核融合科学研究室 HP ; http://www.phys.cst.nihon-u.ac.jp/lab/fusion. [6] 松﨑啓:修士論文「FRC の移送・衝突に関する 2 次元 MHD/Hall MHD シミュレー

ション」, 群馬大学大学院, 2015 年 3 月. [7] E.V. Belova

et al

., Nucl. Fusion 46, 162 (2006).

[8] H. Ohtani

et al

., Phys. Plasmas 10, 0145 (2003).

[9] H. Itagaki: A Doctoral Dissertation, “Studies of stabilization technique for field-reversed configuration by using magnetized plasmoid injection”,

The University of Tokyo, December 2013.

[10] H. Itagaki

et al

.,Phys. Plasmas 21, 030703 (2014).

[11] 石塚高志:修士論文「滑らかな圧力勾配を持つ逆転磁場配位プラズマの平衡に

関する研究」, 群馬大学, 2003 年 3 月.

[12] 宮本健郎「プラズマ物理・核融合」，東京大学出版会，2004 年．

[13] A. Kuwahata

et al.

, Elect. Engi. In Japan, 187, 1 (2014).

[14] H. Ohtani, R. Horiuchi, and A. Ishizawa, J. Plasma Physics 72, 929 (2006).