卒

題

目

業

論

文

永久磁石を用いた磁気ノズル中のプラズマ挙動に関する実験及び
解析

氏

名

花屋倫生

所属教育分野プラズマ理工学

指 導 教 員

中島秀紀

教授

九州大学工学部エネルギー科学科
提出年月

平成 21 年 2 月

目次
第１章

序論 ......................................................................................................................................... 1

1.1 背景 ........................................................................................................................................... 1
1.2 原理 ........................................................................................................................................... 2
1.3 目的 ........................................................................................................................................... 3
第２章

実験 ......................................................................................................................................... 4

2.1 実験装置 .................................................................................................................................... 4
2.2 実験方法 .................................................................................................................................... 5
2.3 実験パラメータ ......................................................................................................................... 6
第３章

3D hybrid code ....................................................................................................................... 9

3.1 はじめに .................................................................................................................................... 9
3.2 基礎方程式 ................................................................................................................................ 9
3.3 計算方法 .................................................................................................................................. 10
3.4 プラズマ粒子の運動 ................................................................................................................ 11
3.5 イオン数密度、電流密度......................................................................................................... 12
3.6 電磁場の計算方法.................................................................................................................... 13
3.7 Digital Filter .......................................................................................................................... 15
3.8 Hybrid code に課せられる制限 .............................................................................................. 16
第４章

実験・シミュレーション結果 ............................................................................................... 17

4.1 実験結果 .................................................................................................................................. 17
4.2 シミュレーション結果 ............................................................................................................ 23
4.2.1 計算パラメータの算出...................................................................................................... 23
4.2.2 拡散項、初期摂動の導入 .................................................................................................. 24
4.2.3 推進効率 ........................................................................................................................... 25
4.3 実験とシミュレーションの比較.............................................................................................. 27
第５章

考察 ....................................................................................................................................... 31

5.1 実験の考察 .............................................................................................................................. 31
5.2 次回の実験についての調査 ..................................................................................................... 32
5.2.1 体系について .................................................................................................................... 33
5.2.2 計算パラメータ ................................................................................................................ 34
5.2.3 結果と考察 ........................................................................................................................ 34
第６章

結論 ....................................................................................................................................... 40

参考文献 .............................................................................................................................................. 41
謝辞 ...................................................................................................................................................... 42

第１章
1.1

序論

背景

地球に恵みを与え続けている太陽は、内部の核融合によってそのエネルギーを生み出している。この核融合エネルギーを実用化しようという試みが、５０年ほど前より行われている。その中で、磁気ノズルによるプラズマの制御技術が考えられている。
磁気ノズルを用いると、核融合により爆発的に膨張するプラズマを磁場との相互作用
により、運動方向を変えて推力を得ることができる。固体壁との相互作用がないため、高い排出速度と大きな推力を達成できる。この磁気ノズルを用いたロケットとしてレーザー核融合ロケットが考案されている。レーザー核融合ロケットの概念図を Fig.1-1 に示す。

Fig.1-1 レーザー核融合ロケットの概念図
過去に、磁気ノズルに関する設計・シミュレーションが行われてきた。Hyde は、超伝導(SCM)コイルを 2 つ用いたモデルを提案している(1)。また、2 次元 MHD シミュレーションコードを用いてプラズマ挙動と推進効率を計算し、1 つのコイルで推進効率65％という結果を得ている(2)。長峯と中島(3)は 3 次元(3D) hybrid code を用いてプラズ
1

マ挙動と推進効率を計算し、推進効率 65％という結果を得ている。Vchivkov(4)は同様のコードで推進効率の最適化を行い、推進効率 70％という結果を得ている。また、実験との比較も行い、推進効率の一致を得ている(5)。坂口は 2 つのコイルを同軸上に並べた計算体系において最適化を行い、3D hybrid code を用いて推進効率 75%という結果を得ている(6)。また、川渕により 3D hybrid code を用いて核融合プラズマの発生位置をずらしたり、２つのコイルを用いることでロケットの姿勢・推進方向制御についての研究が行われた(7)。松田は整形ターゲットに関して、SPH (Smoothed Particle Hydrodynamics)
法と 3D hybrid code を用いて推進効率を高める研究を行った(8)。一方、林田は、原子力ロケット用の磁気ノズルを考案し、3D hybrid code を用いてさまざまな磁気ノズルの体系において推進効率を調べた(9)。

1.2

原理

磁気ノズルの推進原理を Fig.1-2 に示す。

Fig.1-2

磁気ノズルの推進原理
2

プラズマは反磁性体であり、磁場によって運動方向が変化する。コイルによって作られた磁場中において、ターゲットにレーザーを照射することで発生するプラズマは大きなエネルギーをもっている。コイル方向に膨張するプラズマは Larmor 運動をし、半磁性電流が流れるように動きながら磁場を圧縮する。磁場は圧縮されることでその強度を増し、プラズマに対して直角に強い力を及ぼすようになる。これによってプラズマは減速し、やがて停止する。そして、磁場のエネルギーがプラズマのエネルギーよりも強くなると磁場が広がる。それに伴って、プラズマは後方に押し広げられる。このときに、磁場からプラズマへエネルギーが与えられる。このようにして、磁場がプラズマ粒子の運動方向を変化させ、押し出すときの反作用によって推力を得ることができる。

1.3

目的

上記のように、さまざまな体系における設計およびシミュレーションが行われてきた。
しかし、実際にシミュレーションのように、磁場によりプラズマの挙動が変化するのか、実験による確認はほとんどされていなかった。そこで、本研究では、実験を行い、実際にプラズマ挙動が変化するのかを調べることを目的とした。実験は、永久磁石を用いた磁場中で、ターゲットにレーザーを照射してプラズマを発生させ、ICCD カメラによって撮影し、その挙動を観察するというものである。また、実験と同様の体系・パラメー
タにおいて、3D hybrid code を用いてシミュレーションを行った。そして、実験結果とシミュレーション結果を比較した。比較を行い、それに対して考察し、それを踏まえて次回の実験についての考察を行った。

この論文の構成を以下に示す。
第１章

本研究の背景・原理および目的を述べた。

第２章

実験に関する説明を行う。

第３章

3D hybrid code に関する説明を行う。

第４章

実験とシミュレーションの結果を示し、比較する。

第５章

比較による考察、次回の実験に対する考察を行う。

第６章

結論を示す。

3

第２章
2.1

実験

実験装置

実験は、大阪大学レーザー研究所の極端紫外光（EUV:Extreme UltraViolet）生成施
設にて行った。EUV データベースレーザー装置の仕様は以下のようになっている。
EUV 発生用レーザー
Q スイッチ

Nd:YAG レーザー

ビーム数

1 ビーム

繰り返し

10Hz（シングルショットも可能）

波長1.06,0.53,0.36,0.27µm

パルス幅2~10ns（可変）

エネルギー

最大

2J

（10ns/1.06µm）

計測器には、ICCD カメラを２台と分光器を用いた。実験を行ったターゲットチェン
バーを Fig.2-1 に示す。

Fig.2-1 ターゲットチェンバー

4

2.2

実験方法

実験配位を Fig.2-2 に示す。Fig.2-1 のチェンバー内にて行い、真空にして実験を行った。ターゲットにレーザーを照射することでプラズマが発生する。そのプラズマを、ICCD カメラを２台用いて撮影した。Fig.2-2 のように、カメラ１は磁石とターゲットの水平方向(Z-X 方向)から撮影し、カメラ２は磁石と向かい合わせている形(Y-X 方向)で撮影した(Fig.2-3)。
撮影の手法は、ターゲットにレーザーを照射した瞬間を基準とし、一定時間遅らせて撮影する。その遅らせる時間の値を Delay とする。この Delay を決めて撮影することで、一定時間後のプラズマの形を見ることができる。そして、Delay の値を変えつつ撮影することで、プラズマ挙動の変化を観察することができる。

ICCDカメラ２
1cm ターゲット

レーザー

45°

ICCDカメラ１
磁石
分光器

Fig.2-2 実験配位

x

磁石

磁石

x

z
y
ターゲット

1cm

ターゲット

Fig.2-3 （左）カメラ１、（右）カメラ２から見たターゲットと磁石の位置関係
52.3

実験パラメータ

Table 2-1 に、この実験におけるパラメータを示す。
Table 2-1 実験パラメータ
磁石

ネオジウム

磁石の寸法

円柱型

レーザーエネルギー

1.1 J

レーザーパルス

2.3ns

ターゲット

ポリスチレン（CH）

ターゲット直径

100 µmターゲット質量5.2×10-10 kg

ターゲット位置における磁場の強さ

0.1 T

真空度

2.0×10-3 Pa

Φ16 mm×50 mm

実験では、ターゲットは CH であるので、イオンとして C と H が発生しうる。今回は、H イオンは C イオンに比べて質量が十分に小さいため、プラズマ粒子全体の運動量に対する影響は小さいと考え、C イオンのみの挙動を観察することにした。
ターゲットにレーザーを照射し、分光器によって観測された発光分光および C イオ
ンの波長を Fig.2-4 に示す。
ライン

波長 [nm]C+

274.649C+

391.898

C 2+

406.794

C+

426.7

C 2+

464.742C+

515.109C+

588.977C+

723.642

Fig.2-4 （左）発光分光（右）C イオンの波長6２台の ICCD カメラに特定の波長の光を透過するバンドパスフィルターを取り付けた。透過する波長は、カメラ１は 500±5nm、カメラ２は 480±5nm である。Fig.2-4
より、カメラ１では主に C+イオン、カメラ２では主に C2+イオンを撮影する。
また、ネオジウム磁石に関して解析した磁力線を Fig.2-5 に、磁石からの距離に関する磁束密度のグラフを Fig.2-6 に示す。この磁場解析は、解析ソフトの Amaze を用いて求めている。ネオジウム磁石の表面磁束密度は 1.3T である。これらの２つの図より、ターゲットの位置付近の磁束密度は約 0.1T であることがわかる。

Fig.2-5 磁力線

Fig.2-6 磁石からの距離と磁束密度の関係7また、磁場の有無による違いを比較するため、以下の２つのように体系を作り、それぞれ同様の実験を行った。
[１]ネオジウム磁石と同形のアルミ円柱を設置[２]どちらも設置しない

8

第３章
3.1

3D hybrid code

はじめに

この章では、次章以降の解析に使用する、３次元ハイブリッドコードについて述べる。
ハイブリッドと呼ばれる理由は、本手法がプラズマ中のイオンを粒子とし、電子を慣性を無視した流体として取り扱うことに起因する。この計算手法は、系の物理現象の時間スケールが、イオンサイクロトロン周波数で特徴付けられる時間スケールよりも小さい
現象を解析するために開発されたものである(10)。
このコードの特徴は、電子の運動（プラズマ振動、サイクロトロン運動）のタイムスケールで起こる現象の影響を除去し、時間幅や空間幅を大きく出来る点にある。具体的には荷電準中性条件を仮定し、イオンと電子の数密度があらゆるところで等しいと置き、
プラズマ振動を除去する。さらに電磁放射等の高周波現象を無視するために、ダーウィン近似を用いる。また、電磁場はプラズマの運動も含めて、自己無撞着に解かれる。
以後に、本コードの具体的な解法の詳細を示し、最後にコードの妥当性の検証につい
て述べる。

3.2

基礎方程式

計算に使用する方程式を導出するための基礎方程式は次のように与えられる。以下の方程式は直交座標系および、MKSA 単位系で記述されている。まず、粒子として取り扱うイオンの運動方程式は、
midv i
 ZeE  v i  B 
dt
dx i
 vi
dt

(3.1)
(3.2)

である。ここで v i はイオンの速度、 Z はイオンの価数、e は電気素量、 x i はイオンの位置、 m i はイオンの質量、また、 E 、 B はそれぞれ電場、磁場を表している。
流体として扱う電子の運動方程式を電場の解法に用いる。電子の運動方程式は以下の式で表される。
n e medv e
 ene E  v e  B   Pe
dt

(3.3)

ここで v e は電子の速度、 m e は電子の質量、 n e は電子数密度、 Pe は電子圧力である。
電子は慣性無しの流体としてみなすので、慣性による効果は無視し、(3.3)式の左辺は 09にできる( me  0 )。さらに電子圧力 Pe に関しては電子の状態方程式より、
Pe  neTe(3.4)

である。ここで Te は電子温度である。(3.4)式を(3.3)式に代入して Pe を消去する。
0  ene E  v e  B   neTe (3.5)

荷電準中性条件および、ダーウィン近似を用いたアンペアの法則により
  B p   0 J e  J i (3.6)

が得られる。ここで  0 は真空中の誘電率であり、J i 、J e はそれぞれイオン電流密度、
電子電流密度である。また、 B p はプラズマにより生成される磁場で、B p  B  B ext  B pick

(3.7)

とされる。ここで、 B は計算領域内部の合成磁場、 B ext は外部磁場、および B pick はピ
ックアップコイルにより生成される磁場である。ここで、ピックアップコイル磁場は、次式で示される。
  B pick   0 J pick

(3.8)

但し、 J pick はピックアップコイルに流れる電流密度である。
J i はイオンの数密度と電流密度から計算することが出来る。 J e は電子電流密度の定義により、J e  ene v e(3.9)

である。(3.6)、(3.9)式を用いて J e を消去し v e について求めると、

ve  
1  1
   B p  J i ene   0
(3.10)となる。これを(3.5)式に代入して電場 E について整理する。荷電準中性条件より
ene  Zen i を考慮し、

E


1 1
   B p  B  J i  B  ZTe ni 
Zeni   0


(3.11)

を得る。但し、(3.11)式から得られる電場はプラズマが存在する領域内でのみ使用さ
れ、プラズマが存在しない真空領域での電場は後節で述べるラプラス方程式から得られる。
プラズマによる磁場を記述する方程式として、ファラデーの電磁誘導の法則

B
   E
t(3.12)

を用いる。

3.3

計算方法

(3.1)式から(3.12)式を数値計算によって解くためには、基礎式を差分式に直し、数10

値的に安定で、物理的に妥当な空間幅、時間幅で解かなければならない。
計算の大まかな手順を以下に示す。
[１]外部磁場および、イオンの粒子位置、速度の初期設定を行う。
[２]格子点上でのイオンの数密度、電流密度を計算する。
[３](3.11)式により、電場を求める。
[４][１]で設定した外部磁場と、[３]で求めたプラズマによる電場と磁場から、(3.1)
式、(3.2)式を用いて新たなイオン粒子の位置、速度を求める。
[５](3.12)式により磁場を求める。
[２]～[５]を繰り返し、プラズマの時間的変化を追う。
３次元ハイブリッドコードに対する各方程式の差分式を以下の小節で示す。

3.4

プラズマ粒子の運動

本コードでは、いくつかのイオンを超粒子と呼ばれる集団として扱い、この超粒子10 万個を用いた。つまり、超粒子を個別に運動させることでプラズマ粒子の運動を追う。(3.2)式に、leap-frog 法(11)による時間中心差分を施すと、
x in 1  x in
 v in 1 / 2
t(3.13)となる。ここで各変数の上付き添え字は時間を表し、下付き添え字は i 番目のイオンを表す。また、 t は時間幅を表す。
(3.1)式に時間中心差分を施すと、

v in1 / 2  v in1 / 2 Ze  n 1 n1 / 2 v in1 / 2  B n 
E  v i
t
mi 
2






(3.14)となる。 n  1 / 2 はある時刻 n から t / 2 だけ以前の値であることを示す。同様に、n  1 / 2 は t / 2 だけ進んだ値であることを示す。上付き添え字が n 、 n  1 / 2 である変数は既知量である。
プラズマ粒子の運動を求めるためには、まず(3.14)式において、 v in 1 / 2 を左辺にまとめて、未知量である v in 1 / 2 を解く。次に、求めた v in 1 / 2 を(3.13)式の右辺に代入し、新
たな粒子の位置 x in 1 を求めることが可能となる。
(3.14)式を変形すると、Av n 1 / 2  S(3.15)となる。  、  、  、  を次のように定義する。

11



qt n
Bx
2 mi



qt n
By
2 mi

 

qt n
Bz
2 mi

  1 2   2   2
なお、ここでの E 、B は粒子位置におけるそれぞれの値を示している。これらを用いて、行列 A とベクトル S は以下のように構成される。
 1

A 
 



1

 

 
1 



 vixm  cE xn  v ym  v zm S   viym  cE yn  v zm  v xm  v m  cE n  v m  v m 
z
x
y 
 iz

(3.16)(3.17)

ここで、 c  qt mi 、 m  n  1 / 2 である。さらに、(3.15)式を次のように変形して、速度を求める式が得られる。

v n1 / 2  A 1S(3.18)

ここで、 A 1 は以下のように得られる。
 1 21A 1       3.5

      
1   2       1   2 

(3.19)

イオン数密度、電流密度

イオン数密度を求めるためには、空間内に広がっているプラズマを格子点上に割り付けなければならない。そのために、PIC（Particle In Cell）法を用いる。今、１つの超粒子の中に N 個のイオンが含まれている場合を考える。Fig. 3-1 に示すように、格子点内に存在する超粒子の位置によって、8 個の格子点によって囲まれた体積 V の立方体
を、 V1 、 V 2 、～ V 8 の体積を持つ直方体に分割する。格子点 i, j, k  に割り振られる超粒子に含まれるイオンの数は、

ni (i, j , k ) V2
NV(3.20)

で与えられる。すべての粒子についての和を ni についてとれば、ある格子点上の n i を定義することが出来る。
12

イオン電流密度 J i は、数密度 n i を用いて、
J i  Zen i v i(3.21)

より求めることが出来る。

Fig. 3-1 PIC 法

3.6
(1)

電磁場の計算方法
プラズマ領域での計算

電磁場を求めるためには、(3.11)式、(3.12)式を用いる。(3.12)式を片側差分して、B n1 / 2 だけを左辺に残すと、
B n1 / 2  B n t  En2(3.22)となる。また、(3.11)式は、

En 


T
1 1
1 n

  B np  B n 
J i  B n  e nin ni   0 ZeZee


(3.23)である。
ここで、 B n 、 E n は既知であるとすると、 (3.22)式から B n1 / 2 を求めることが出来る。次に、 B n1 / 2 、 J n 1 / 2 、 v n 1 / 2 、 nin 1 / 2 を(3.23)式に代入して E n1 / 2 を求め、線形外挿11より電磁場の予測値 Enpred、 B npredを次式より求める。

1
Enpred En  2En1 / 2
1
B npred B n1 / 2 (3.24)t1
  E npred2
13(3.25)

1
1
ここで、下付き添え字 pred は予測値であることを表している。 Enpred、B npred
を求めた

3/ 2
3/ 2
3/ 2
3 / 2
後、それらによるイオンの運動を解き、 v inpred、ninpred、J inpredを求める。次に、B npredを

次式により求める。
3 / 2
1
B npred B npredt  E n1 / 2
2(3.26)

3 / 2
そして、(3.23)式で E n を求めた時と同様にして、 E npredを求める。最後に次式に示す

線形内挿から次のタイムステップの値、 E n1 、 B n1 を求める。

E n1 1 n1 / 2 1 n3 / 2E E pred22(3.27)t  E n1
2

(3.28)

B n 1  B n 1 / 2 

これらの値 E n1 、 B n1 より、 v in  3 / 2 、 nin  3 / 2 、 J in  3 / 2 を求めることが出来る。
ここで、電場を求める際に用いた一連の手法は、予測子－修正子法と呼ばれる収束計算の一種である。この方法では、n タイムステップから n  1 タイムステップの値を求めるために、 n  3 / 2 タイムステップの予測値を求め内挿計算を行う。そして、求めた

n  1 タイムステップの値を再度予測値として用いる。具体的には、(3.28)式により求め1られた B n1 を次の予測値 B npredとし、(3.24)式から(3.28)式までを次式で定義する条件(12)を満足するまで、繰り返し計算を行う。

max i , j ,k Bin, j ,1k( q 1)  Bin, j ,1k( q )max i , j ,k Bin, j ,1k( q )(3.29)

ここで q は収束計算の回数を表し、 i 、 j 、 k はそれぞれ x、y および z 方向の格子番号を表している。  は収束計算を打ち切るための値である。
電場を計算するために用いる(3.23)式は、分母に数密度を含んでいるため、 n i がゼロとなる領域ではこの式は発散してしまう。さらに、数値計算の特性上 n i があまり小さい値になると、発散もしくは計算誤差が生じる。このことを避けるために、数密度がある特定の値 n v 以上の場合においてのみこのような電場の計算を行う。
(2)

真空領域での電場の計算

真空領域では、ラプラス方程式  2 E  0 を SOR（Successive Over Relaxation）法
を用いて解く。SOR における加速係数  は正方形領域の場合は次のように表される。

21  sin  NGMX (3.30)

ここで NGMX は x、y または z 方向の格子数である。 の値は通常 1.8 から 1.9 付近である。また、   1の場合には、SOR 法はガウス－ザイデル法に等しい。
14

SOR 法では、次式によって変数の収束を加速する。

E p (i, j , k )   E p (i  1, j , k )  E p (i  1, j , k ) E p (i, j  1, k )  E p (i, j  1, k ) E p (i, j , k  1)  E p (i, j , k  1) 6(3.31)

 (1   ) E p (i, j , k ) , ( p  x, y, z )次式で定義される残差  x , y , z p  E p (i  1, j , k )  E p (i  1, j , k ) E p (i, j  1, k )  E p (i, j  1, k ) E p (i, j , k  1)  E p (i, j , k  1)

(3.32)

 6 E p (i, j , k ) , ( p  x, y, z )が、あらかじめ決められた値以下になれば真空部分の電場の計算を終了する。

3.7

Digital Filter

PIC 法で格子点上の量を決める場合、例えば粒子の位置を格子上に割り振る場合には、メッシュ内に存在する粒子の数が尐ない領域ではその付近の格子に比べて数密度が小さくなり、数密度の変化が不連続になる可能性がある。このような現象を避けるために Digital Filter と呼ばれる一種の補間を施すことがある。Digital Filter は、格子上に定義される値、例えば電場、数密度等を計算する際に発生する高周波の雑音を除去する目的にも用いられる(13)。Digital Filter の式は以下のように表される。

newi , j ,k

cornerM i , j ,k  S iside, j , k  K i , j , k

M  6 S  12 K

(3.33)

 i , j ,k を格子点上に定義される値とし、 inew, j , k は Digital Filter を通した新しい i , j , k を表
corner
す。また、  iside, j , k と  i , j ,k は side term と corner term であり、それぞれ以下の式で表される。
 iside, j , k   i 1, j , k   i 1, j , k   i , j 1, k  i , j 1,k   i , j ,k 1   i , j ,k 1(3.34)

 icorner, j , k   i , j 1, k 1   i , j 1, k 1   i , j 1, k 1
  i , j 1,k 1   i 1, j ,k 1   i 1, j ,k 1
  i 1, j ,k 1   i 1, j ,k 1   i 1, j 1,k(3.35)

  i 1, j 1,k   i 1, j 1,k   i 1, j 1,k

ここで、 M 、 S 、 K の値をそれぞれ 12、4、1 とすると(3.33)式は次式となる。
1
4

inew
i , j ,k 
, j ,k 

1 side1
 i , j ,k   icorner, j ,k12
48(3.36)

数密度、イオン電流密度および電場を計算した後に、それぞれの値にこの Digital
15

Filter を作用させる。

3.8

Hybrid code に課せられる制限

ハイブリッドコードにおいて、差分化の時間幅 t 、空間幅 x は次の条件を満たす必要がある。
CFL（Courant-Friedrichs-Lewy）条件
t 

x
VA

(3.37)

これは、 t , x  0 で解が収束する条件である。 V A はアルフベン速度であり、VA 

B0 0 n i mi

(3.38)

で定義される。また、付加的な数値的安定条件は

 ci t  0.2(3.39)である。ここで、  ci はイオンサイクロトロン周波数である。
メッシュ幅の最小値を求める条件として
x  c  pi

(3.40)がある。この c  pi より x が小さい場合にもハイブリッドコードは発散することなく計算を行えるが、その解が正しいものである保証は無い。ここで、  pi はイオンプラズ
マ周波数である。

16

第４章
4.1

実験・シミュレーション結果

実験結果

実験によって撮影された画像を Fig.4-1 から Fig.4-5 までに示す。Fig.4-1、Fig.4-2
ではネオジウム磁石を設置した、つまり磁場がある時に ICCD カメラ１にて、Delay
0ns から 2µs までの範囲で撮影した画像である。Fig.4-3 は同様に磁場がある時にICCD カメラ２にて、Delay250ns、400ns、750ns の時に撮影した画像を並べている。
一方、Fig.4-4 はネオジウム磁石の代わりにアルミ円柱を設置した、つまり磁場が無い
時にそれぞれ ICCD カメラ１、カメラ２にて撮影した画像であり、Delay は 10ns、25ns、50ns である。また、Fig.4-5 はネオジウム磁石もアルミ円柱も設置していないときに、
ICCD カメラ１にて、Delay 0ns、20ns、50ns において撮影した画像である。
Fig.4-1、Fig.4-2 を見てみると、磁石の表面付近にプラズマ粒子が見られるが、これは発生したプラズマ粒子の一部が高速で磁石表面のほうに移動したためだと考えられる。Fig.4-1、Fig.4-2 のように、磁場があるときには Delay 2µs においてもプラズ
マ粒子が残っているように見られたが、アルミ円柱を設置した時、また、何も設置してないときには、Delay 100ns 以降においてはプラズマ粒子はほとんど見られなかっ
た。

17

磁石表面
ターゲット位置

Fig.4-1ネオジウム磁石を設置し、カメラ１にて撮影した画像(15mm×15mm)
上から

Delay 0ns、200ns、500ns

18

磁石表面
ターゲット位置

Fig.4-2ネオジウム磁石を設置し、カメラ１にて撮影した画像(15mm×15mm)
上から

Delay 1µs、1.5µs、2µs

19ターゲット位置

Fig.4-3 ネオジウム磁石を設置し、カメラ２にて撮影した画像(20mm×20mm)
上から

Delay 250ns、400ns、750ns

20アルミ表面
ターゲット位置

Fig.4-4 アルミ円柱を設置して撮影した画像（左）カメラ１、15mm×15mm（右）カメラ２、20mm×20mm
上から Delay 10ns、25ns、50ns

21

ターゲット位置

Fig.4-5

何も設置せず、カメラ１にて撮影した画像(15mm×15mm)
上から

Delay 0ns、20ns、50ns

22

4.2
4.2.1

シミュレーション結果
計算パラメータの算出

実験結果との比較のため、3D hybrid code を用いてシミュレーションを行った。 実
験の体系に近似するため、計算パラメータを算出した。
イオン質量は、ターゲットが CH であるため、C と H の質量の平均をとって、(12+1)/2=6.5 とした。
イオン価数は C の価数が 4、H の価数が 1 より、その平均をとって 2.5 とした。
プラズマ半径は実験結果の Fig.4-5 の Delay 0ns における画像より 1.5mm とした。
また、プラズマエネルギーは、レーザーエネルギーの 10%となると仮定し、レーザー
エネルギーは 1.1J なので、 1.1 0.1  0.11[ J ]とした。

以上において算出した値を含めた計算パラメータを Table 4-1 に示す。
Table 4-1

計算パラメータ

イオン質量6.5 amu

イオン価数

2.5

時間ステップ

68 ps

プラズマ質量5.0×10-10kg

プラズマ半径

1.5mm

プラズマエネルギー

0.11J

電子温度

0eV

計算領域

120mm×120mm×120mm

メッシュ数

120mm×120mm×120mm

粒子数100000磁石によって形成される初期磁場の配置は Amaze を用いて作成した。
ここで、シミュレーションの時間の基準について述べる。シミュレーションにおいてプラズマ粒子は、座標(x,y,z)=(0,0,0)より発生するものとなっている。Table 4-1 より、プラズマ半径は 1.5mm である。また、プラズマ発光 Fig.4-5 の Delay 0s のときの画像より、プラズマ膨張速度は約 40km/s と求められた。これより、
1.5  10 3
 40 ns
40  10 3

つまり、t=40ns のときが、シミュレーションにおけるプラズマが発生してからの初23

期状態となる。

4.2.2

拡散項、初期摂動の導入

ここでは、実験の条件により近似するため、磁場拡散の導入、初期摂動の導入を行った(14)。
まず、磁場拡散の影響を考慮し、これを拡散項の形で前章で述べたコードに組み込む改良(15)を行った（次式下線部）。
コード中で改良を行った基礎方程式は、電子の運動方程式に拡散項を加えたものである。
dv ene (E  v e  B J )  Pe
dt

(4.1)


T
1 1
1

(  B)  B  J i  B  e ni  Jni  0 Ze
Ze
e


(4.2)

ne me

よって、電場を求める式は、E

さらに、運動保存則から、イオンの運動方程式は、
dv i Ze(E  v i  Β   J )
dt
mi

(4.3)

dx i
 vi
dt

(4.4)

me
 ei
e 2 ne(4.5)となる。ここで、  は、

 と表される。 ei は電子-イオン衝突周波数である。Zakharov 等は ei を、電子サイクロ
トロン周波数  ce に比例する実効衝突周波数 eff によって置き換えることで実験結果を説明できると述べている(16)。

 eff   ce(4.6)

Zakharov は、この  の最適値をおよそ 0.3 と定めており、今回はこれに従った。
次に、初期摂動について述べる。Fig.4-3 の Delay250ns の時の画像を見ると、その形は等方的ではなく、尐し歪んでいるように見られる。これより、この形に近似する
ため、初期摂動を加えることにした。画像を参照して、周方向（X-Y 方向）にモード
数 m=3、それらの大きさがプラズマ初期半径の 20%となるように、初期プラズマを分
布させた。

24

4.2.3

推進効率

シミュレーションによって求められるプラズマ挙動に関して、推進効率ηを以下の式により評価した。

 mv （z方向運動量） 100[%]
 m v （全初期運動量）
z

(4.7)

0

これより、ηはプラズマ粒子全体の z 方向速度が大きいほどその値は大きくなる。
ηが大きいほど性能が高いと言える。

以上のことを踏まえて、3D hybrid code を用いたシミュレーション結果を Fig.4-6に示す。計算の結果、t=2.08µs のとき、η=49.5%となった。

25

磁石

Fig.4-6

シミュレーション結果（100mm×100mm）（左）Z-X 方向、（右）Y-X 方向
上から

t=40ns、1.06µs、2.08µs

26

4.3

実験とシミュレーションの比較

ネオジウム磁石を設置した実験結果とシミュレーション結果を比較するため、実験によって得られた画像とシミュレーションによる計算結果を並べた図を Fig.4-7、
Fig.4-8、Fig.4-9 に示す。比較のため、シミュレーションによって得られた図は拡大されている。また、4.1 節で述べたように、時間の基準を合わせるため、
Delay  t-40ns
この条件に注意して比較を行った。

27

(4.8)

Fig.4-7 実験とシミュレーションの比較（Z-X 方向）（左）実験による画像(15mm×15mm)（右）解析結果(25mm×25mm)

28

上から Delay 0s、200ns、500ns

上から t= 40ns、176ns、516ns

Fig.4-8 実験とシミュレーションの比較（Z-X 方向）（左）実験による画像(15mm×15mm)（右）解析結果(25mm×25mm)

上から Delay 1µs、1.5µs、2µs

上から t= 0.992µs、1.468µs、2.012µs

29

Fig.4-9 実験とシミュレーションの比較（Y-X 方向）（左）実験による画像(20mm×20mm)（右）解析結果(20mm×20mm)

30

上から Delay 250ns、400ns、750ns
上から t=244ns、380ns、720ns

第５章
5.1

考察

実験の考察

ここでは、第４章で示した実験結果とシミュレーション結果、及びこれらの比較により考察されることを述べる。
（1） 実験における磁場の有無による比較
Z-X 方向に関して、磁場が有るときの Fig.4-1、Fig.4-2 の画像と、磁場が無いときのFig.4-4、Fig.4-5 のカメラ１の画像を用いて比較する。Fig.4-1、Fig.4-2 より、磁場があるときには、プラズマ粒子は磁場によって +Z 方向に押されているように見られ、Delay2µs においてもプラズマ粒子が残っているように見られた。一方、Fig.4-4、Fig.4-5より、磁場が無いときには Delay100ns 以降においては、プラズマ粒子はほとんど撮影されなかった。これより、プラズマ挙動は磁場の有無によって大きく変化することが分かった。
Y-X 方向に関して磁場が有るときの Fig.4-3 の画像と、磁場が無いときの Fig.4-4 のカ
メラ２の画像を比較しても、同様のことが分かった。また、Fig.4-3 の Delay250ns、400nsの画像より、磁場があるときにはプラズマ挙動の形は等方的に広がっているとは容易に言えず、尐し歪んでいるように見られる。これより、レイリー・テーラー不安定性が発生している可能性があると思われる。
レイリー・テーラー不安定性とは流体力学的なものであり、例えば密度の大きい流体
が密度の小さい流体に上から接し、重力加速度が下向きにあって平衡状態であるとき、その境界は不安定であり、擾乱が生じるとその乱れが大きくなる現象である。今回の場
合、磁場が密度の小さい流体、プラズマが密度の大きい流体の役割となり、磁場によるプラズマ膨張の減速のため、加速度の向きがプラズマから磁場になる時に発生しうる。
この不安定性が起きると、プラズマの形の歪みが時間とともに徐々に成長し、プラズマが磁場を突き破ってしまい、推力が小さくなってしまう。よって、レイリー・テーラー
不安定性が存在する可能性を考慮しなければならない。
（2） 実験とシミュレーションの比較
Fig.4-7、Fig.4-8 の Z-X 方向に関して考察する。
Delay 200ns ではシミュレーションでは磁石表面のほうに粒子は無いものの、ターゲットの位置付近においてはその挙動は類似している。
Delay 500ns において、実験による画像では粒子はターゲットの位置付近に多く見ら31

れるが、シミュレーションでは粒子はターゲットと磁石の間でも見られ、若干の違いが見られる。これについて、シミュレーションでは初期プラズマ粒子は Z 方向では等方的に広がることを仮定としている一方で、実験ではカメラ１の位置から見ると、レーザーのターゲットへの照射方向は+Z 方向であり、１方向なので、それによってプラズマ粒子の初期運動エネルギーが+Z 方向に尐し大きくなっていると考えられ、この相違点によって違いがあると思われる。
Delay 1µs においては実験とシミュレーションにおいてほぼ一致しているように見られる。
Delay 1.5µs、2µs において、実験による画像ではターゲットの位置付近に留まっている粒子が多く見られるが、シミュレーションではターゲットの位置から Z 方向+0.5cm
から+1cm のところに多く留まっているのが見られた。しかし、全体的にみたプラズマ
粒子の分布は実験とシミュレーションともに類似しているように見られる。
次に、Fig.4-9 の Y-X 方向に関しての考察を行う。プラズマ粒子の広がり方や形の大きさはほぼ一致しているように見られる。しかし、実験による画像では、粒子全体でみた分布はターゲット位置から左のほうにずれているように見られる。この原因として、
カメラ２の位置から見ると、ターゲットへ照射するレーザーの方向は右から左、つまり-X 方向であるため、それによって-X 方向にずれたと考えられる。
また、実験結果において、Delay250ns の時には形は歪んでいるように見られたが、Delay400ns、750ns の時の画像を見ると、その歪みは成長しているようには見られなかった。よって、レイリー・テーラー不安定性は発生していなかったと思われる。同様に、シミュレーション結果からも不安定性が発達しているようには見られなかった。しかし、今回の結果のみで、発生していないと明確には言えないため、今後の実験や解析によって明らかにしていく必要がある。
以上の考察より、実験結果とシミュレーション結果では、プラズマ粒子全体の分布について多尐のずれがあるものの、その広がり方や分布した形は類似しているように見られる。従って、実験とシミュレーションはその結果において大きな違いはなく、シミュ
レーションによって実験体系の検討、及びその結果の予測が十分にできるものだと考え
られる。

5.2

次回の実験についての調査

ここでは、次回の実験の計画を立てるために、どのようなパラメータで行うかをシミュレーションにより検討する。今回の実験結果とシミュレーションとの比較により、実際のプラズマ挙動はシミュレーションとほぼ一致していると考えられることが分かっ32

た。しかし、シミュレーション結果では推進効率はη=49.5%となり、これまでに行われてきた研究結果によるものと比べると小さい値となってしまった。そこで、次回の実験の計画を立てる１つとして、実験の条件を踏まえつつ、パラメータを変えて、シミュレーションによって算出される推進効率が高い体系を調査する。

5.2.1

体系について永久磁石によって作られる磁場は磁石の表面付近では強いが、Fig.2-6 のように、表面から離れると磁場は急速に弱まってしまうため、より強い磁場の中で実験を行うのは困難である。そこで、次回の実験では、磁場をコイルによって作ることにし、永久磁石の代わりにコイルを設置した体系を考える。想定する体系の概略図を Fig.5-1 に示す。
ここでは、レーザーエネルギーやターゲットに関するパラメータは実験のものと同じにする。変化させるパラメータは、コイルとターゲットとの距離 L 、コイル半径 r 、ター
ゲット位置における初期磁場 B 0 である。

レーザー

I

45°
L
r

ターゲット位置に置ける磁場B0

Fig.5-1 想定する体系
ここで、コイル半径の上限について述べる。コイル半径がむやみに大きいと、照射するレーザーの直線上にコイルがあたってしまうという問題が起こる。照射するレーザー
方向と、コイルとターゲットが向きあう方向がなす角度は Fig.5-1 のように 45°である
ので、問題を避ける条件として次の式を満たさなければならない。
(5.1)

r＜L

また、 B 0 は r , L ,コイルに流れる電流 I によって次の式で求められる。

B0 

 0 Ir 22(r 2  L2 ) 3 / 2(5.2)

ただし、  0 は真空透磁率である。これを変形するとI

2 B0 (r 2  L2 ) 3 / 2
0 r 233(5.3)

(5.3)式により、 B 0 を設定することにより I を求められる。 I の値が小さいと、実験においてその電流を流すのが容易になる。

5.2.2

計算パラメータ

共通する計算パラメータは Table 5-1 のようになる。
Table 5-1 計算パラメータ
イオン質量6.5 amu

イオン価数

2.5

時間ステップ54.3 ps

プラズマ質量5.2×10-10 kg

プラズマ半径

1.5mm

プラズマエネルギー

0.11 J

電子温度

0 eV

計算領域

100mm×100mm×100mm

メッシュ数

100×100×100

粒子数100000

Table 5-1 のパラメータを用い、コイルの位置 L 、半径 r 、及びターゲット位置における磁束密度 B 0 を変えた体系でシミュレーションを行う。なお、今回の調査においては、
拡散項の導入や、初期摂動の付加は行っていない。
体系（Ⅰ）～（Ⅳ）における L , r , B 0 を以下に示す。
（Ⅰ） L  1cm, r  0.5cm, B0  0.1T
（Ⅱ） L  1cm, r  0.5cm, B0  0.2T
（Ⅲ） L  1cm, r  0.8cm, B0  0.1T
（Ⅳ） L  2cm, r  1cm, B0  0.1T
(4.7)式、(5.3)式より、 , I をそれぞれの体系において求めた。計算は 40000time stepの 2.17µs+40ns=2.21µs まで行った。また、指標として I の上限を 5kA、推進効率ηの
目標を 60%以上とする。

5.2.3

結果と考察

各体系における計算結果を以下に示す。
体系Ⅰの計算結果を Fig.5-2 に示す。I =8.9kA、t=2.21µs のとき、 =52.1%となった。
体系Ⅱの計算結果を Fig.5-3 に示す。 I =17.8kA、t=2.21µs のとき、η=67.2%となっ34

た。
体系Ⅲの計算結果を Fig.5-4 に示す。 I =5.2kA、t=2.21µs のとき、η=53.1%となった。
体系Ⅳの計算結果を Fig.5-5 に示す。 I =17.8kA、t=2.21µs のとき、η=66.2%となった。
いずれの体系においても電流値の上限を超えてしまった。一方、この４つの体系のう
ち、目標の値を超えた推進効率となったのは体系Ⅱと体系Ⅳであった。しかし、コイルに流すべき電流値が上限を大きく超えており、実験で行うのは難しい。
電流値が最も小さくなったのは体系Ⅲの場合であり、目標の値をわずかに超えるものとなった。体系Ⅲは体系Ⅰと比べてコイル半径が大きいという点しか異なってないものの、電流値は小さくなり、推進効率はほぼ等しい値となった。これより、コイル位置とターゲット位置における磁場を固定した場合、コイル半径はできるだけ大きいほうが電流値が小さくなり、推進効率も変わらないので、実験により適した体系になると考えられる。
また、体系Ⅱは体系Ⅰの B 0 を、体系Ⅳは体系Ⅰの L と r をそれぞれ２倍にし、電流値が同じになるようにしたが、推進効率は 1%ほどしか差がなかった。
以上のことより、電流値が大きいときに推進効率は大きくなるが、そのときのコイル
位置とターゲット位置における磁場を固定し、コイル半径を大きくすることで、電流値を小さくでき、推進効率の値をほとんど変化させずにおける可能性があることが分かった。今後、このことを考慮しつつ、電流値が 5kA 以内、推進効率が 60%以上になるよ
うに、更なるパラメータの検討を行っていく必要がある。

35

Coil

Coil

Coil

Fig.5-2

体系Ⅰの計算結果（80mm×80mm）（左）Z-X 方向、（右）Y-X 方向

上から t= 40ns、1.12µs、2.21µs
36

Fig.5-3

体系Ⅱの計算結果（80mm×80mm）（左）Z-X 方向、（右）Y-X 方向

上から t= 40ns、1.12µs、2.21µs
37

Coil

Coil

Coil

Fig.5-4

体系Ⅲの計算結果（80mm×80mm）（左）Z-X 方向、（右）Y-X 方向

上から t= 40ns、1.12µs、2.21µs
38

Coil

Coil

Coil

Fig.5-5

体系Ⅳの計算結果（80mm×80mm）（左）Z-X 方向、（右）Y-X 方向

上から t= 40ns、1.12µs、2.21µs
39

第６章

結論

本研究では、磁気ノズルの実現可能性を検討する目的で、磁場中のプラズマ粒子の挙動が変化するのかを実験により確かめ、その結果とシミュレーションによる結果を比較した。また、次回の実験についての調査も行った。
実験を行った結果、レーザーによって生成されるプラズマ粒子の挙動は磁場中において大きく変化することが分かった。また、3D hybrid code を用いたシミュレーション結果と比較してみたところ、プラズマ挙動は大きく異なることはなく、ほぼ一致しているということが分かった。しかし、シミュレーションによって計算された推進効率は 49.1%と低い値だった。
この結果より、次回の実験の計画を立てる一環として、コイルを用いたシミュレーシ
ョンを行い、推進効率の向上を試みた。その結果、体系Ⅱが最もよいという結果になったが、電流値が大きくなってしまった。

今後、次回の実験に向けての準備、及び体系を想定していき、その体系におけるシミュレーションを行って、実験可能なパラメータ、高い推進効率であるかの調査を行っていくことが課題となる。

40

参考文献
(1) R. A. Hyde, L. L. Wood, Jr.,J. H. Nuckolls: AIAA Paper, No.72-1063, (1972)(2) R. A. Hyde: UCRL-88857, (1983)(3) Y. Nagamine, H. Nakashima: Fusion Technology, 35 (1999), 62(4) K. V. Vchickov, H. Nakashima, F. Ichikawa, Y. P. Zakharov: Vacuum, 73 (2003),
427(5) K. V. Vchickov, H. Nakashima, Y. P. Zakharov, T. Esaki, T. Kawano, T. Muranaka:
Jpn. J. Appl. Phys., 42 (2003), 6590
(6) N. Sakaguchi, Y. Kajimura and H. Nakashima: Trans. Japan Soc. Aero. Space Sci.,
48 (2005), 180(7) 川渕 亮,’’磁気ノズルを用いたレーザー核融合ロケットの推力ベクトルの制御技術’’, 九州大学工学部エネルギー科学科卒業論文 (2006)(8) 松田 伸夫,’’整形ターゲットを用いたレーザー核融合ロケットの推進効率の向上’’, 九州大学工学部エネルギー科学科卒業論文 (2007)(9) 林田 憲治, ‘’原子力ロケット用磁気ノズルの考案とプラズマ挙動の数値解析’’, 九州大学
工学部エネルギー科学科卒業論文 (2008)(10) D. S. Harned: J. Comp. Phys., 47 (1982), 452(11) E. J. Horowiz, D. E. Schumaker, and D. V. Anderson: J. Comp. Phys., 84 (1989),
279
(12)

A. Mankofsky, R. N. Sudan, and J. Denavit, J. Comp. Phys., 70 (1987), 89(13) C. K. Bridsall and A. B. Langdon: ‘’Plasma Physics via Computer Simulation ’’,
IOP Publishing Ltd, (1997)(14) 村中 崇信,’’実験室プラズマの巨視的構造形成に関する数値解析’’, 九州大学総合理
工学府博士論文 (2000)(15) K. R. Moore, V. A. Thomas, and D. J. McComas, J. Geophys. Res.,96 (1991) ,A5,
7779(16) Y. P.Zakharov, et al. J. Plasma and Fusion Res., SERIES 2 (1999) 398

41

謝辞

本研究において、研究の場を提供して頂き、また御指導頂いた中島秀紀教授に深く感謝致します。研究に関して的確な御指導をして頂いた渡辺幸信准教授に深く感謝致します。研究室生活の基礎的なことから様々なことについて御指導して頂いた山本直嗣助教に深く感謝致します。快適な研究環境を提供して下さった大神めぐみ秘書に深く感謝致します。
今回の実験をするにあたって、実験のできる環境を提供して頂き、さらに忙しい中、実験の打ち合わせ、及び作業をして下さった大阪大学レーザー研究所の藤岡慎介氏、砂原淳氏、城崎知至氏、光産業創成大学院大学の森芳孝氏に深く感謝致します。
同じ核融合ロケットグループでは、計算機やパソコンのことについて多くのアドバイ
スをして下さった松田伸夫氏に深く感謝致します。研究内容について理解ができていな
かった私に色々と教えて下さった前野旭弘氏に深く感謝致します。
先輩として様々なアドバイスをして頂いたり、休憩時における雑談などで充実した研究生活にして頂いた鶴哲平氏、本岡親英氏、永田英隆氏、大宅将史氏、小谷優介氏、新谷將氏、内藤有紀氏、安部晋一郎氏、渡辺健人氏、叶涛氏、蘭長林氏に深く感謝致します。また、同じ学部生として共に頑張ってきた江崎徹氏、西山貴章氏、平山嵩祐氏に深く感謝致します。
最後に、時々会って楽しい一時を与えてくれた友人、そして今まで支えて下さった家族に深く感謝致します。

42