卒

題

目

氏

業

論

文

原子力ロケット用磁気ノズルの考案とプラズマ挙動の数値解析

名

林田憲治

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

指 導 教 員

中島秀紀

教授

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

平成 20 年 2 月

目次
第1章

序論 ................................................................................................................................2

1.1

背景 ...................................................................................................................................2

1.2

目的 ...................................................................................................................................3

第2章

原子力ロケット ................................................................................................................5

2.1

推進原理 ...........................................................................................................................5

2.2

原子力ロケット ...................................................................................................................6

第3章

3D hybrid code ...............................................................................................................7

3.1

はじめに.............................................................................................................................7

3.2

基礎方程式 .......................................................................................................................7

3.3

計算方法 ...........................................................................................................................8

3.4

Hybrid code に課せられる制限 ......................................................................................14

第4章

解析結果 ......................................................................................................................15

4.1

コイル設計 .......................................................................................................................15

4.2

評価方法 .........................................................................................................................16

4.3

計算体系Ⅰ .....................................................................................................................16

4.4

計算体系Ⅱ .....................................................................................................................22

4.5

計算体系Ⅲ .....................................................................................................................29

4.6

計算体系Ⅳ .....................................................................................................................32

4.7

計算体系Ⅴ .....................................................................................................................38

第5章

結論 ..............................................................................................................................41

参考文献
謝辞

1

第1章
1.1

序論
背景

46 億年前に太陽が誕生して以来、太陽は地球に恵みを与え続けてきた。その無尽蔵とも思
われる膨大なエネルギーは、太陽内部における核融合によって作り出されている。今日、核融合は新しいエネルギー源として注目され、盛んに研究が行われている。
この研究の中で、磁気ノズルというプラズマの制御技術が考えられている。磁気ノズルを用いると核融合により爆発的に膨張するプラズマを磁場との相互作用を利用し、固体壁と相互作用すること無く運動方向を変えて推力を得ることが出来る。この技術を用いたロケットとしてレーザ
ー核融合ロケットが考えられており、他の宇宙推進システムと比較して、高い排出速度（即ち高い推進材利用効率）と大きな推力を同時に達成可能である。概念図を Fig.1-1 に示す。

Fig. 1-1

レーザー核融合ロケットの概念図

磁気ノズルは過去に設計、シミュレーションが行われてきた。Hyde は、超伝導(SCM)コイルを2 つ用いたモデルを提案している(1)。また、2 次元 MHD シミュレーションを用いてプラズマ挙動と推進効率を計算し、1 つのコイルで推進効率 65％という結果を得ている(2)。長峯と中島(3)は 3次元(3D) hybrid code を用いてプラズマ挙動と推進効率を計算し、推進効率 65％という結果を得ている。Vchivkov(4)は同様のコードで推進効率の最適化を行い、推進効率 70％という結果を得ている。また、実験との比較も行い、推進効率の一致を得ている(5)。坂口は 2 つのコイルを同軸上に並べた計算体系において最適化を行い、3D hybrid code を用いて推進効率 75%とい2

う結果を得ている(6)。また、川渕により 3D hybrid code を用いて核融合プラズマの発生位置をずらしたり、２つのコイルを用いることでロケットの姿勢・推進方向制御についての研究が行われた(７)。
しかし、未だ核融合推進の実用化の目処は立っていない。そこで、レーザー核融合ロケットの前段階のシステムとして、小型原子炉を電源とした、原子力ロケットが提案されている。原子力ロケットの概念図を Fig.1-2 に示す。

Fig. 1-2

1.2

原子力ロケットの概念図

目的

原子力ロケットでは、レーザー核融合ロケットにおけるこれまでの研究成果を生かして、推進システムに磁気ノズルを用いる。さらに磁気ノズルに冗長性を持たせるために、コイルを複数搭載する。しかし、レーザー核融合ロケットと違い、原子力ロケットでは磁気ノズルの中心軸上に機体があるため、これまで考えられてきた、プラズマが中心軸上を抜けてしまう計算体系では、
機体をプラズマから保護する必要が出てくる。(Fig.1-3 参照) また、原子力ロケットでは核融合反応は生じさせないので、中性子の遮蔽を考える必要はない。そのため、原子力ロケットに適した磁気ノズルを新しく考えることが望まれる。
本研究では、原子力ロケット用の磁気ノズルを考案し、3D hybrid code を用いて磁気ノズル
内でのプラズマ挙動を計算し、磁気ノズルの構成がプラズマ挙動に与える影響を調べる。

3

Coil

Coil

Fig. 1-3

磁気ノズル内でのプラズマ挙動（上から 0、11.08 s）

この論文の構成は以下の通りである。
第 1 章 本研究の背景と目的を述べた。
第 2 章 原子力ロケットについて述べる。
第 3 章 本研究に用いた計算コードについて述べる。
第 4 章 解析結果について述べる。
第 5 章 結論を示す。
4

第2章
2.1

原子力ロケット
推進原理

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

[a]

[b]

[c]

Fig. 2-1

磁気ノズルの推進原理

5

2.2

原子力ロケット

現在の探査機には、電源として主に太陽電池が使われているが、将来、木星以遠の深宇宙探査を目的とした場合、太陽から距離があるため、太陽電池では十分な電力を得ることが困難である。従って、太陽からの距離に依存しない電源が求められている。
そこで、将来の宇宙推進システムとして核融合ロケットが有望視されている。このシステムでは太陽からの距離に依存しないだけでなく、高い比推力と比出力も同時に得ることができる。
本研究グループではレーザー核融合ロケットについて研究を行ってきた。しかし、レーザー核
融合ロケットには、レーザー核融合を実現するためのレーザードライバや、大型超伝導コイル、
さらに 250 ton もの巨大な建造物を宇宙で組み立てるための技術などが必要であり、実用化の目処は立っていない。
このような状況下で、その中間の段階である原子炉を搭載した推進システムは米国、ロシア
等で、1950 年代から研究が行われてきた(８)。さらに、2004 年には原子炉を電源とした原子力電
気推進システム計画、JIMO(Jupiter Icy Moon Orbiter)が計画された。一方、日本では、まだ宇宙推進用原子炉を本格的に検討した例はない。
本研究グループが提案する原子力ロケットでは、原子炉は高温ガス炉を採用し、電気への変換は MHD 方式を、推力発生には磁気ノズルを用いる。また、本研究グループで行った概念設計では、機体総質量は 6000 kg 程度になるという計算結果が導き出されている。

6

第3章
3.1

3D hybrid code
はじめに

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

3.2

基礎方程式

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

mi

dv i
 ZeE  v i  B 
dt

(3.1)

dx i
 vi
dt

(3.2)

ここで v i はイオンの速度、 Z はイオンの価数、 e は電気素量、 x i はイオンの位置、 mi はイ
オンの質量、また、 E 、 B はそれぞれ電場、磁場を表している。
流体として扱う電子の運動方程式を電場の解法に用いる。電子の運動方程式は以下の式で表される。

ne medv e ene E  v e  B   Pe
dt

(3.3)

ここで v e は電子の速度、 me は電子の質量、 n e は電子数密度、 Pe は電子圧力である。電子は慣性無しの流体としてみなすので、慣性による効果は無視し、(3.3)式の左辺は 0 にできる
( me  0 )。さらに電子圧力 Pe に関しては電子の状態方程式より以下のようになる。
7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  Zeni を考慮し次式を得る。

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

(3.11)

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

B   E
t

3.3(3.12)

計算方法(3.1)式から(3.12)式を数値計算によって解くためには、基礎式を差分式に直し、数値的に安定で、物理的に妥当な空間幅、時間幅で解かなければならない。
計算の大まかな手順を以下に示す。
8

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

3.3.1

プラズマ粒子の運動

本コードでは、いくつかのイオンを超粒子と呼ばれる集団として扱い、この超粒子 10 万個を用いた。つまり、超粒子を個別に運動させることでプラズマ粒子の運動を追う。(3.2)式に、
leap-frog 法(10)による時間中心差分を施すと次式となる。

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

 、  、  、  を次のように定義する。

9(3.15)



qt n
Bx
2 mi



qt n
By
2 mi

 

qt n
Bz
2 mi

(3.16)

  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.17)(3.18)

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

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

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

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

(3.20)

イオン数密度、電流密度

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

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

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

10

イオン電流密度 J i は、数密度 ni を用いて次式により求めることが出来る。

J i  Zeni v i(3.22)

Fig. 3-1 PIC 法

3.3.3

電磁場の計算方法(1) プラズマ領域での計算
電磁場を求めるためには、(3.11)式、(3.12)式を用いる。(3.12)式を片側差分して、 B n 1 / 2 だ
けを左辺に残すと次式となる。

B n1 / 2  B n t  En2(3.23)

また、(3.11)式は以下の式になる。

En 

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

(3.24)

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

11(3.25)

1
B npred B n1 / 2 t
1  E npred2(3.26)

1
1
ここで、下付き添え字 pred は予測値であることを表している。 E npred、 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.27)

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

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

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

(3.29)

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.29)式により求められた B n1 を次の予測値
1
B npredとし、(3.25)式から(3.29)式までを次式で定義する条件 (11)を満足するまで、繰り返し計算を行う。
max i , j ,k Bin, j ,1k( q 1)  Bin, j ,1k( q )max i , j ,k Bin, j ,1k( q )(3.30)

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

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

21  sin  NGMX 
12(3.31)

ここで NGMX は x、y または z 方向の格子数である。  の値は通常 1.8 から 1.9 付近である。
また、   1 の場合には、SOR 法はガウス－ザイデル法に等しい。
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.32)

 (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.33)

 6 E p (i, j , k ) , ( p  x, y, z )3.3.4

Digital Filter

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

newi , j ,k

cornerM i , j ,k  S iside, j , k  K i , j , k(3.34)
M  6S  12 K
 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.35)

 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.36)

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

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

1
4

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

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

数密度、イオン電流密度および電場を計算した後に、それぞれの値にこの Digital Filter を13作用させる。

3.4

Hybrid code に課せられる制限

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

t 

x
VA

(3.38)

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

VA 

B0
 0 ni mi(3.39)

また、付加的な数値的安定条件は以下である。

 ci t  0.2(3.40)

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

(3.41)

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

14

第4章
4.1

解析結果
コイル設計

まず初めに、磁気ノズルを構成する超伝導コイルの個数が増えることによる質量増加の影響を考えるために、超伝導コイルの質量を計算した (13)。超電導電磁石用線材としては、NbTi と
Nb3Sn が実用的である。ここでは、NbTi を採用する。
安定化材には、熱伝導率が高くて軽い Al を用い、複合構造とする。NbTi と Al の断面積比を 1:2 とする。NbTi の臨界温度は 9 K なので、冷却材として液体 He を用いる。構造は、バンド
ル導体型とし、線材は有効冷却面積を大きくするために、より線にし、容器におさめ、冷却材は素線の隙間を通す。液体 He 温度（4.2 K）では、磁場強度を 2 T とすれば、NbTi の臨界電流密度は、105～106 A/cm2 となる。コイル素線 1 本 1 本には、線材が張り付けられ、素線間の影響が出ないようにしている。しかし、コイルの堆積の 40％は冷却材や絶縁材で占められ、残り 60％を
コイル素線と安定化材が占められるものとする。コイル電流 I c 、コイル素線部の断面積 S csc 、及びコイル質量 mcoil は、次式のようになる。
Ic 

2 B0 Rc

0

(4.1)

Ic
Jc

(4.2)

mcoil  2Rc S csc (  Nb3Sn  2  Al )(4.3)

S csc ここで、 J c はコイル素線部の臨界電流密度、 Rc はコイル半径、  NbTi は NbTi の密度、  Al は Alの密度であり、 Rc  0.5 m、  NbTi  8.57 g/cm3、  Al  2.7 g/cm3 を用いる。また、冷却材と絶縁材の重量は無視する。Jc は臨界電流密度以下である必要があり、ここでは安全係数を 10 として、
Jc = 50000 A/cm2 を用いる。
次に、コイルカン（容器）の質量を求める。コイルは、それ自身の作り出す磁場によって、外側へのフープ力を受ける。その力は、次式により求められる。

Fhoop  8R
13
 0 I c 2  ln c  
2
4
 ac

(4.4)

ここで、 Fhoop はコイルの受ける力を表し、 I c はコイル電流、 Rc はコイル半径、 a c はコイル断面円形とした場合の断面半径を表す。このフープ力により、コイルはフープ応力を受けるが、コイルカンによって、そのフープ応力を吸収させる。コイルカンには、グラスファイバーを用いる。グラスファイバーの密度は  gf  1.89 g/cm3 であり、また使用応力として  u  330 MPa 用いる。
15

コイルカン断面積 S cc 、重量 mcc は次式で求まる。

S cc 

Fhoop
2 u

(4.5)

mcc  2Rc S cc  gf

(4.6)

その結果、コイル質量とコイルカンの質量を合わせると、26 kg となった。原子力ロケットの機体総質量はおよそ 6000 kg なので、コイルを増やすことによる機体質量増加の影響はかなり小さいと言える。

4.2

評価方法

今回、磁気ノズルにおけるプラズマ粒子の挙動をシミュレーション結果を評価する指標として、
推進効率ηとαを用いた。
推進効率ηとαは以下の式で求めた。

 mv  100  全Z方向運動量   100 [ % ]全初期運動量 
mv

(4.7)

機体方向－Z方向に抜けたプラズマ粒子 数
100 [ % ]全プラズマ粒子数(4.8)

推進効率  z0α

つまり、はどれだけプラズマの運動方向を＋Z 方向に変化させられているか、はどれだけ機
体方向にプラズマが抜けてしまっているかを表す。は高いほど、は小さいほど性能が良いと
言える。

4.3

計算体系Ⅰ

4.3.1

計算体系Ⅰ－A

計算体系を Fig.4-1 に示す。この計算体系では、コイルを縦に置くことにより、－Z 軸方向に抜けるプラズマ粒子をなくすことを目的としている。
コイルは半径 1 m、 Z-Y 平面上にある。初期プラズマは原点を中心に 0.3 m の球体内に一様に分布しており、プラズマ粒子は原点を中心にして放射状に初速を持って広がる。計算領域
は Z 軸方向に-7 m～7 m、X,Y 軸方向にそれぞれ-6 m～6 m の直方体とした。
計算に用いたパラメータを Table 4-1 に示す。主な値は長峯(3)と同じものを用いた。計算は、が収束したとみなせる 40000 time step の 11.08 s まで計算した。
16

計算体系

Fig. 4-1
Table 4-1

計算に用いるパラメータ

Coil radius [m]1.0Coil current [MA]3.57

Coil position along Z [m]-1.0

Plasma coordinates (x,y,z) [m](0, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]12.0×12.0×14.0

Mesh number

120×120×140

Number of particles

4.3.2

100000

結果と考察

計算結果を Fig.4-2 に示す。11.08 s では、＝49.1%、＝0.0%だった。
－Z 軸方向に広がるプラズマはコイル付近の磁場が強い領域で跳ね返されることで、－Z 方
向に抜けるプラズマ粒子はなくなっているが、X、Y 方向に大きくプラズマ粒子が広がっている
ため、は低い値になっていると思われる。
17

C1

C1

C1

Fig. 4-2

計算結果（左は Z-X 平面での粒子位置と磁場、右は X-Y 平面での粒子位置）
上から順に、0、5.54、11.08 s

18

4.3.3計算体系Ⅰ－B

計算体系を Fig.4-3 に示す。この計算体系では、2 つのコイルを縦に並べ、コイルの間に磁場の壁を作ることによって、－Z 軸方向に抜けるプラズマ粒子を減尐させることを目的としている。
コイルはそれぞれ半径 1 m、 Z-Y 平面に平行である。初期プラズマは原点を中心に 0.3 mの球体内に一様に分布しており、プラズマ粒子は原点を中心にして放射状に初速を持って広がる。計算領域は Z 軸方向に-7 m～7 m、X,Y 軸方向にそれぞれ-6 m～6 m の直方体とした。
計算に用いたパラメータを Table 4-2 に示す。主な値は長峯(3)と同じものを用いた。計算は、が収束したとみなせる 40000 time step の 11.08 s まで計算した。

Fig. 4-3

計算体系

19

Table 4-2計算に用いるパラメータ
C1 : 1.0

Coil radius [m]C2 : 1.0
C1 : 3.57

Coil current [MA]C2 : 3.57
C1 : (1.0, 0.0, -2.0)

Coil position (x,y,z) [m]C2 : (-1.0, 0.0, -2.0)

Plasma coordinates (x,y,z) [m](0, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]12.0×12.0×14.0

Mesh number

120×120×140

Number of particles

4.3.4

100000

結果と考察

計算結果を Fig.4-4 に示す。11.08 s では、＝47.1%、＝0.4%だった。
－Z 方向に広がるプラズマ粒子が、コイルの間にできた磁場の壁によって跳ね返されている様子が確認できた。しかし、X、Y 方向に大きくプラズマ粒子が広がってしまっているため、ηは低い値となっている。
この計算体系では、片方のコイルが壊れた場合、コイルを平行移動させることにより、計算体
系Ⅰ－A と同じ体系とすることができる。

20

C1
C2

C1
C2

C1
C2

Fig. 4-4

計算結果（左は Z-X 平面での粒子位置と磁場、右は X-Y 平面での粒子位置）
上から順に、0、5.54、11.08 s

21

4.4

計算体系Ⅱ

4.4.1

計算体系Ⅱ－A

計算体系を Fig.4-5 に示す。この計算体系では、2 つのコイル（C1 と C2）を横に並べ、その中心軸上にプラズマを置くことにより、－Z 方向の運動量を持つプラズマ粒子がコイル付近の、磁場が強い空間へ広がる。その結果、の減尐と、の向上が見込めるのではないかと考えた。
コイルはそれぞれ半径 1 m で、X-Y 平面に平行におき、図中の矢印の向きに電流を流す。
初期プラズマは原点を中心に 0.3 m の球体内に一様に分布しており、プラズマ粒子は原点を中心にして放射状に初速を持って広がる。計算領域は Z 軸方向に-7 m～7 m、X,Y 軸方向にそれぞれ-6 m～6 m の直方体とした。
計算に用いたパラメータを Table 4-3 に示す。主な値は長峯(3)と同じものを用いた。コイル間の距離を L とし、L＝30, 40, 50 cm と変化させ、計算を行った。また、L＝30 cm の時に、C1 の電流の向きを逆に流した場合も計算を行った。計算は、が収束したとみなせる 40000 time
step の 11.08 s まで計算した。

L

Fig. 4-5

計算体系

22Table 4-3計算に用いるパラメータ
C1 : 1.0

Coil radius [m]C2 : 1.0
C1 : 3.57

Coil current [MA]C2 : 3.57
C1 : -1.0

Coil position along Z [m]C2 : -1.0

Plasma coordinates (x,y,z) [m](0, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]12.0×12.0×14.0

Mesh number

120×120×140

Number of particles

4.4.2

100000

結果と考察

L＝30 cm の時の計算結果を Fig.4-6 に示す。11.08 s では＝73.8%、＝3.0%であった。
L＝40 cm の時、11.08 s では＝71.0%、＝3.3％であった。
L＝50 cm の時、11.08 s では＝69.7%、＝4.6%であった。
L が大きくなるにつれ、は減尐し、は増加している。これは、コイル間の距離が大きくなる
につれ、コイルの内側を抜けるプラズマ粒子は減尐するが、それ以上にコイルの間を抜けるプ
ラズマ粒子が増加するためだと考えられる。また、X-Y 平面の図より、プラズマ粒子が Y 方向に
大きく広がっているのが分かる。そのため、この Y 方向への広がりを抑えることで運動量を＋Z方向に集中させ、推進効率をより向上させることができると思われる。
L＝30 cm として、C1 の電流の向きを逆にした時の計算結果を Fig.4-7 に示す。11.08 s では＝69.7%、＝4.3%であった。
C1 の電流の向きを逆にしたことにより、コイルの間に磁場の壁ができている。その結果、計算体系Ⅰ－A と比べて、コイルの間を抜けるプラズマ粒子は減尐しているが、両コイルの内側を抜けるプラズマ粒子は増加している。
これらの計算体系では、片方のコイルが壊れた場合、コイルの位置を平行移動することで、
Fig.1-3 と同じ体系にすることができる。このとき、＝、＝となる。
23

C1
C2

C1
C2

C1
C2

Fig. 4-6

L＝30 cm での計算結果（左は Z-X 平面での粒子位置と磁場、右は X-Y 平面での粒子位置）
上から順に、0、5.54、11.08 s
24

C1
C2

C1
C2

C1
C2

Fig. 4-7

計算結果（Z-X 平面での粒子位置と磁場）上から順に上から順に、0、5.54、11.08 s

25

4.4.3

計算体系Ⅱ－B

計算体系を Fig.4-8 に示す。この計算体系では、計算体系Ⅱ－A において L＝30 cm の状
態から両コイルをそれぞれ内側に 30°傾けている。こうすることで、X 方向に広がるプラズマ粒子を抑え、の向上を目指す。
コイルはそれぞれ半径 1 m で、 電流は図中の矢印の向きに流す。初期プラズマは原点を中心に 0.3 m の球体内に一様に分布しており、プラズマ粒子は原点を中心にして放射状に初速を持って広がる。計算領域は Z 軸方向に-7 m～7 m、X,Y 軸方向にそれぞれ-6 m～6 m の直方体とした。
計算に用いたパラメータを Table 4-4 に示す。主な値は長峯(3)と同じものを用いた。計算は、が収束したとみなせる 40000 time step の 11.08 s まで計算した。

Fig. 4-8

計算体系

26

Table 4-4計算に用いるパラメータ
C1 : 1.0

Coil radius [m]C2 : 1.0
C1 : 3.57

Coil current [MA]C2 : 3.57
C1 : (1.15, 0.0, -1.0)

Coil position (x,y,z) [m]C2 : (-1.15, 0.0, -1.0)

Plasma coordinates (x,y,z) [m](0, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]12.0×12.0×14.0

Mesh number

120×120×140

Number of particles

4.4.4

100000

結果と考察

計算結果を Fig.4-9 に示す。11.08 s では＝68.6%、a＝1.6%であった。
計算体系Ⅱ－A と比べて、X 方向のプラズマ粒子の広がりにはほとんど変化がなかった。また、は減尐しているがも下がっている。そこで、Z が－1.0～1.0 m の領域のプラズマ粒子数を数えたところ、計算体系Ⅱ－A と比べて増加していた。つまり、その分＋Z 方向に大きな速度を持つプラズマ粒子の数が尐ないということである。そのため、は下がっているのだと考えられ
る。

27

C1
C2

C1
C2

C1
C2

Fig. 4-9

計算結果（左は Z-X 平面での粒子位置と磁場、右は X-Y 平面での粒子位置）
上から順に、0、5.54、11.08 s

28

4.5

計算体系Ⅲ

4.5.1

計算体系

計算体系を Fig.4-10 に示す。この計算体系では、計算体系Ⅰで分かった、Y 方向にプラズマが広がってしまうという問題を、3 つのコイル（C1、C2、C3）をコイルの中心が正三角形の頂点の位置にくるように配置することにより防ぎ、結果としての向上との減尐を目的としている。
コイルはそれぞれ半径 1 m、 X-Y 平面に平行である。初期プラズマは原点を中心に 0.3 mの球体内に一様に分布しており、プラズマ粒子は原点を中心にして放射状に初速を持って広がる。計算領域は Z 軸方向に-7 m～7 m、X,Y 軸方向にそれぞれ-6 m～6 m の直方体とした。
計算に用いたパラメータを Table 4-5 に示す。主な値は長峯(3)と同じものを用いた。コイル間の距離 L を L=30, 40, 50 cm と変化させ、計算を行った。計算は、が収束したとみなせる
40000 time step の 11.08 s まで計算した。

L

計算体系

Fig. 4-10

29

Table 4-5計算に用いるパラメータ
C1 : 1.0

Coil radius [m]C2 : 1.0
C3 : 1.0
C1 : 3.57

Coil current [MA]C2 : 3.57
C3 : 3.57
C1 : -1.0

Coil position along Z [m]C2 : -1.0
C3 : -1.0

Plasma coordinates (x,y,z) [m](0, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]12.0×12.0×14.0

Mesh number

120×120×140

Number of particles

4.5.2

100000

結果と考察

L＝30 cm の時、11.08 s では＝77.5%、＝2.4%であった。
L＝40 cm の時の計算結果を Fig.4-11 に示す。11.08 [μs]では＝78.0%、＝2.7％であった。
L＝50 cm の時、11.08 s では＝76.3%、＝3.3%であった。
計算体系Ⅰと比べ、の増加が確認できた。これは、コイルを 3 つにしたことにより、Y 方向へのプラズマ粒子の広がりを抑えることができたためと考えられる。しかし、に関しては、計算体系Ⅰと比べてそれほど減尐させることはできなかった。このことから、コイルを増やしても Z 軸上にできた隙間から抜けるプラズマ粒子を減らす効果はそれほど得られないと考えられる。
この体系では、コイルが 1 つ壊れた場合、コイルを平行移動することで、計算体系Ⅱを同じにすることができる。

30

Fig. 4-11

L＝40 cm での計算結果（左は Z-X 平面での粒子位置と磁場、右は X-Y 平面での粒子位置）
上から順に、0、5.54、11.08 s
31

4.6

計算体系Ⅳ

4.6.1

計算体系

計算体系を Fig.4-12 に示す。この計算体系では、C1 の中心軸上にプラズマをおき、C1 の中心を抜けたプラズマ粒子を C2 が作る磁場によって曲げることでを減尐させ、さらに C3 の作る
磁場によってプラズマ粒子の運動を＋Z 方向に変えることにより、の向上を目的としている。
コイルの半径はそれぞれ 1 m で、図中の矢印の向きに電流を流す。C1 と C2 は X-Y 平面に平行で、C3 は Z-Y 平面上にある。プラズマは X 軸と C1 の中心軸の交点にを中心に半径 0.3 mの球体内に一様に分布している。プラズマ粒子は交点を中心にして放射状に初速を持って広がる。計算領域は Z 軸方向に-7 m～7 m、X,Y 軸方向にそれぞれ-6 m～6 m の直方体とした。
計算に用いたパラメータを Table 4-6 に示す。主な値は長峯(3)と同じものを用いた。計算は、が収束したとみなせる 100000 time step の 27.70 s まで計算した。

計算体系

Fig. 4-12

32

Table 4-6計算に用いるパラメータ
C1 : 1.0

Coil radius [m]C2 : 1.0
C3 : 1.0
C1 : 3.57

Coil current [MA]C2 : 3.57
C3 : 3.57
C1 : （1.5, 0.0, -1.0）)
C2 : （0.0, 0.0, -2.5）

Coil position (x, y, z) [m]C3 : （-1.5, 0.0, -1.0）
Plasma coordinates (x,y,z) [m](1.5, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]12.0×12.0×14.0

Mesh number

120×120×140

Number of particles

4.6.2

100000

結果と考察

計算結果を Fig.4-13、Fig.4-14 に示す。27.70 s において=73.0%、=0.0%という結果が得られた。これは、－Z 方向に抜けていくプラズマ粒子が磁場によって運動方向を変えられ、＋Z方向の運動量が増加したためと考えられる。問題としては、－X 側に広がったプラズマ粒子がC3 の作った磁場によって引き戻されてしまっている。そのため、引き戻されるプラズマ粒子をなくすことでさらに推進効率の向上が見込める。そこで、引き戻されるプラズマ粒子を減らすために、整形ターゲットを用いて計算を行った。整形ターゲットとは、ターゲットを整形することにより、プラズマの噴出方向を変える方法である[14]。計算結果を Fig.4-15 に示す。－X 側に広がるプラ
ズマ粒子を減らすことはできなかった。また、＋X 方向への広がりが大きくなってしまったため、推進効率は下がった。今後、整形ターゲットの形を変えて、計算を行う必要がある。
この計算体系では、C3 が壊れた場合も推進効率はわずかに下がるが、＝0.0%という結果が得られる。その場合の計算結果を Fig.4-16 に示す。
33

C1
C2
C3

C1
C2
C3

C1
C2
C3

Fig. 4-13

計算結果（Z-X 平面での粒子位置と磁場）上から順に、0、5.54、11.08 s

34

C1
C2
C3

C1
C2
C3

C1
C2
C3

Fig. 4-14

計算結果（Z-X 平面での粒子位置と磁場）上から順に、16.62、22.16、27.70 s

35

C1
C2
C3

C1
C2
C3

C1
C2
C3

Fig. 4-15

整形ターゲットを用いた計算結果（Z-X 平面での粒子位置と磁場）上から順に、0、5.54、11.08 s

36

Fig. 4-16

計算結果（Z-X 平面での粒子位置と磁場）上から順に、0、5.54、11.08 s

37

4.7

計算体系Ⅴ

4.7.1

計算体系

計算体系を Fig.4-17 に示す。これは、過去に坂口により最適化された計算体系(6)を基に、2つのコイル（C1 と C2）の中心軸上後方に小さなコイル（C3）を置くことで、を減尐させることを目的としている。
コイルはそれぞれ X-Y 平面に平行である。初期プラズマは原点を中心に 0.3 m の球体内に一様に分布しており、プラズマ粒子は原点を中心にして放射状に初速を持って広がる。計算領
域は Z 軸方向に-5 m～5 m、X,Y 軸方向にそれぞれ-4 m～4 m の直方体とした。
計算に用いたパラメータを Table 4-7 に示す。主な値は坂口(6)と同じものを用いた。計算は、が収束したとみなせる 40000 time step の 11.08 s まで計算した。

計算体系

Fig. 4-17

38

Table 4-7計算に用いるパラメータ
C1 : 2.0

Coil radius [m]C2 : 1.0
C3 : 0.2
C1 : 1.19

Coil current [MA]C2 : 3.57
C3 : 1.19
C1 : -0.5

Coil position along Z [m]C2 : -1.0
C3 : -1.75

Plasma coordinates (x,y,z) [m](0, 0, 0)

Plasma radius [m]0.3

Plasma energy [MJ]

4

Plasma mass [mg]110

Atomic mass [AMU]197

Atomic number

79Time step Δt [ns]0.277Calculation region [m]8.0×8.0×10.0

Mesh number

74×74×93

Number of particles

4.7.2

100000

結果と考察

計算結果を Fig.4-18 に示す。11.08 s では＝81.4%、＝3.2%であった。
また、比較のために Table 4-7 のパラメータを用いて、C3 がない体系で計算すると、＝80.8%、＝3.4%であった。
後方に半径の小さなコイルを置くことによるの減尐はわずかしか得られなかった。コイル半径を小さくしても、コイルの中心を抜けるプラズマ粒子を防ぐことは難しいと考えられる。

39

C1
C2
C3

C1
C2
C3

C1
C2
C3

Fig. 4-18

計算結果（Z-X 平面での粒子位置と磁場）上から順に、0、5.54、11.08 s

40

第5章

結論

本研究では、原子力ロケット用の磁気ノズルを考案し、3D hybrid code を用いて磁気ノズル
内でのプラズマ挙動を計算し、磁気ノズルの構成がプラズマ挙動に与える影響を調べた。また、システムの冗長性についても考察を行った。
計算体系Ⅰでは、－Z 方向に抜けるプラズマ粒子を減らすために、コイルを Z-Y 平面に平行に 1 つ置いた場合と 2 つ置いた場合において計算を行った。いずれも－Z 方向に抜ける粒子を減らすことはできたが、X、Y 方向の運動量が大きくなってしまい、推進効率は低い値となった。
計算体系Ⅱでは、システムに冗長性を持たせると同時に、－Z 方向に抜けるプラズマ粒子を減らすために、2 つのコイルの間にプラズマを置いて、計算をおこなった。また、X 方向に広がるプラズマ粒子を減らすために両コイルを内側に 30°傾けた場合の計算も行った。その結果、－Z 方向に広がるプラズマ粒子は、コイル付近の磁場が強い領域で跳ね返され、－Z 方向に広がるプラズマ粒子の減尐と推進効率の向上を同時に達成できた。しかし、コイルを傾けた体系では、X 方向に広がるプラズマ粒子を抑える効果はあまり得られず、逆に＋Z 方向に大きな速度を持つ粒子が減ってしまい、推進効率は下がった。いずれも、プラズマ粒子が Y 方向に大きく広がってしまうという問題が分かった。
計算体系Ⅲでは、計算体系Ⅱのプラズマ粒子が Y 方向に大きく広がるという問題をなくすた
めに、コイルを 3 つに増やして計算を行った。その結果、プラズマ粒子の Y 方向への広がりを抑えることができ、推進効率の向上と、若干ではあるがαの減尐が達成できた。
計算体系Ⅳでは、磁場によってプラズマ粒子の運動方向を曲げることで、－Z 方向に抜ける
プラズマ粒子の減尐と推進効率の向上を同時に達成することを目的として計算を行った。その結果、磁場によってプラズマ粒子の運動方向が曲げられ様子が確認でき、ηの向上とαの減尐が確認できた。また、－X 方向に広がったプラズマ粒子の一部が、磁場によって引き戻されているため、これを防ぐために整形ターゲットを用いて計算を行ったが、引き戻されるプラズマ
粒子を減らすことはできなかった。今後、ターゲットの形を変えて計算を行う必要がある。
計算体系Ⅴでは、坂口により最適化された計算体系(6)を元に、2 つのコイルの後方に半径の小さなコイルを置くことで、－Z 方向に抜けるプラズマ粒子の減尐を目的として計算を行った。
その結果、半径の小さなコイルを置くことによるプラズマ挙動の変化はほとんど見られなかった。
41

以上の結果より、原子力ロケットに用いる磁気ノズルの候補としては、冗長性を重視すると、計算体系ⅢかⅣが適していると考えられる。そのうち、機体に衝突するプラズマ粒子の影響が大きいとすると、計算体系四の方がが小さいため適していると思われる。一方、機体に衝突するプラズマ粒子の影響が小さいならば、推進効率の高い計算体系Ⅲの方が適していると思われる。しかし、いずれの場合も、1 つのコイルが壊れた場合、性能の低下はまぬがれない、そこで、システムに冗長性を持たせる手段として、計算体系を 1 台として、同じものをもう 1 台搭載するという考えもある。
今回は、大まかなプラズマ挙動を調べ、これまでの計算結果と比較するために、レーザー核融合のパラメータを用いて計算を行っている。そのため、今後、より実際の現象に近づけるために、原子力ロケットのパラメータを用いて計算を行う必要がある。

42

参考文献
(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) 川渕 亮,“磁気ノズルを用いたレーザー核融合ロケットの推力ベクトルの制御技術” H18.2(8) M. S. EI-Genk, “A Critical Review of Space Nuclear Power and Propulsion, 1984-1993”
American Institute of Physics Press, (1994)(9) D. S. Harned: J. Comp. Phys., 47 (1982), 452(10) E. J. Horowiz, D. E. Schumaker, and D. V. Anderson: J. Comp. Phys., 84 (1989), 279(11) A. Mankofsky, R. N. Sudan, and J. Denavit, J. Comp. Phys., 70 (1987), 89(12) C. K. Bridsall and A. B. Langdon: "Plasma Physics via Computer Simulation", IOP
Publishing Ltd, (1997)(13) 中島 秀紀 他 “レーザー核融合ロケット” (核融合研究第 66 巻, 第 3 号, 1991)(14) 松田 伸夫,“整形ターゲットを用いたレーザー核融合ロケットの推進効率の向上” H19.2

43