卒 業 論 文

題

目

氏

磁気ノズルを用いたレーザー核融合ロケットの推力ベクトルの制御技術

名川渕亮

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

指 導 教 員中島秀紀 教授

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

平成 18 年 2 月

目次
第1章

序論····························································································································· 2

1.1

背景 ································································································································ 2

1.2

目的 ································································································································ 2

第2章

レーザー核融合推進システム ···················································································· 42.1

推進原理 ························································································································ 42.2

レーザー核融合ロケット·································································································· 5

第3章

3D hybrid code ············································································································ 6

3.1

はじめに·························································································································· 63.2

基礎方程式 ···················································································································· 6

3.3

計算方法 ························································································································ 7

3.3.1

プラズマ粒子の運動································································································ 83.3.2

イオン数密度、電流密度 ························································································· 93.3.3

電磁場の計算方法 ································································································ 10

3.3.4

Digital Filter··········································································································· 12

3.4

Hybrid code に課せられる制限 ···················································································· 13

第4章

解析結果··················································································································· 14

4.1

評価方法 ······················································································································ 144.2制御方法 1 ··················································································································· 14

4.2.1

計算体系 ··············································································································· 14

4.2.2

結果と考察············································································································· 164.3制御方法 2 ··················································································································· 23

4.3.1

計算体系 ··············································································································· 23

4.3.2

結果と考察············································································································· 244.4制御方法 3 ··················································································································· 32

4.4.1

計算体系 ··············································································································· 32

4.4.2

結果と考察············································································································· 35

第5章

結論··························································································································· 36

参考文献 ···································································································································· 38
謝辞············································································································································ 39

1

第1章
1.1

序論
背景

46 億年前に太陽系が誕生して以来、地球に恵みを与え続けてきた太陽。その無尽蔵とも思われる膨大なエネルギーは、太陽内部における核融合によって作り出されてきた。今日、核融合は新しいエネルギー源として注目され、盛んに研究が行われている。
この研究の中で、磁気ノズルというプラズマの制御技術が考えられている。磁気ノズルを用いると核融合により爆発的に膨張するプラズマを磁場との相互作用を利用し、固体壁との相互作用無く運動方向を変えて推力を得ることが出来る。従って、他の宇宙推進システムと比較して、高い排出速度（即ち高い推進材利用効率）と大きな推力を同時に達成可能である。この
磁気ノズルを用いたロケットとしてレーザー核融合ロケットが考えられている。概念図を Fig.
1-1 に示す。
磁気ノズルは過去に設計、シミュレーションが行われてきた。Hyde は、超伝導(SCM)コイル
を 2 つ用いたモデルを提案している(1)。また、2 次元 MHD シミュレーションを用いてプラズマ挙動と推進効率を計算し、1 つのコイルで推進効率 65％という結果を得ている(2)。長峯と中島(3)は 3 次元(3D) hybrid code を用いてプラズマ挙動と推進効率を計算し、推進効率 65％という結果を得ている。Vchivkov(4)は同様のコードで推進効率の最適化を行い、推進効率 70％という結果を得ている。また、実験との比較も行い、推進効率の一致を得ている(5)。

1.2

目的

上に挙げた研究は推進効率を求めることが目的のため、推力方向は一方向に固定されていた。しかし、ロケットを実際に運用するためには推進方向を制御する必要がある。
現在、宇宙空間で行われている姿勢・推進方向制御はおもに主推進装置とは独立した装置によって行われている。一方、レーザー核融合ロケットは巨大な構造物であるため、この方法では制御装置も大きくなってしまうと思われる。そこで、主推進装置を用いた姿勢・推進方向の制御ができれば、ロケットの効率的な運用につながる。
そこで、本研究では 3D hybrid code を用い、3 通りの制御方法で磁気ノズルにおける推力ベクトル制御の可能性と課題を探った。
一つ目は Orth が提案した制御方法(6)について 3D hybrid code を用い、初期プラズマ発生
位置と推力ベクトル方向の関係を求めることで姿勢・推進方向制御の可能性を探った。この方法を『制御方法 1』とした。
二つ目は坂口によって最適化された(7)、2 つのコイルを使った磁気ノズルについて、一方のコイルを傾けることで推力ベクトルの制御を目指した。この方法を『制御方法 2』とした。
2

三つ目は、広がるプラズマの運動量を集中化することで推力ベクトル制御を目指した。この方法を『制御方法 3』とした。このシミュレーションは、江藤、大村、加藤、久保田、松田との共同研究として行った。
この論文の構成は以下の通りである。
第 1 章 本研究の背景と目的を述べた。
第 2 章 レーザー核融合推進システムについて述べる。
第 3 章 本研究に用いた計算コードについて述べる。
第 4 章 解析結果について述べる。
第 5 章 結論を示す。

Fig. 1-1 レーザー核融合ロケットの概念図

3

第2章
2.1

レーザー核融合推進システム
推進原理

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

Fig. 2-1 磁気ノズルでの推進原理

4

2.2

レーザー核融合ロケット

将来の宇宙輸送手段として、レーザー核融合ロケットは魅力的なロケットである。1972 年、
Hyde によって始めて超伝導(SCM)コイルを 2 つ用いたレーザー核融合ロケットのモデルが提案された(1)。ロケットの推進システムにエネルギー密度が大きい核融合を用いると、化学推進に比べて高い比推力と高い比出力を同時に得ることができる。推進材の排出速度は 107m/sec
程度になると見られ、比推力（排出速度を重力加速度で割った値）は 106sec となる。この値は、
化学推進ロケットの 2000 倍程度であり、ペイロード（積載貨物）の割合が多く取れることになる。
さらに、ロケットの外部空間をエンジンとして利用するので、その分宇宙船に占めるエンジンの質量が小さいので高比出力が達成されると考えられる。これは、輸送時間の短縮にもつながり、
宇宙線による被爆など人体への悪影響も抑えられる。
しかし、レーザー核融合ロケットの実現には多くの課題がある。レーザー核融合を実現するためのレーザードライバ、大型超伝導コイル、大量の核融合燃料生成などの開発が急務であるが、これらは現在活発に研究が行われている地上発電用核融合技術がそのまま流用できる。
さらに、2003 年に Orth が提案したレーザー核融合ロケットの概念設計(6)では、宇宙船全体の質量が 6500ton もあり、2005 年に本研究グループが行った新しい概念設計でも 250ton の巨
大な宇宙船になるという計算結果が導き出された。この巨大な構造物を宇宙で組み立て、さらに高速での運転を行うための技術も必要である。

5

3D hybrid code

第3章
3.1

はじめに

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

3.2

基礎方程式

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

dv e= −ene (E + v e × B ) − ∇Pe
dt

(3.3)

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

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)式を数値計算によって解くためには、基礎式を差分式に直し、数値的に安定で、物理的に妥当な空間幅、時間幅で解かなければならない。
計算の大まかな手順を以下に示す。
7

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

3.3.1

プラズマ粒子の運動

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

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
mi ⎩
2
∆t
⎭

(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

α 、 β 、 γ 、 δ を次のように定義する。

8(3.15)

α=

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 である。さらに、 (2.15)式を次のように変形して、速度を求める式が得られる。

v n +1 / 2 = A −1S(3.18)

ここで、 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 を定義することが出来る。
9イオン電流密度 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

10(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)式までを次式で定義する条件(10)を満足するまで、繰り返し計算を行う。
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 があまり小さい値になると、発散もしくは計算誤差が生じる。このことを避けるために、数密度がある特定の値 nv 以上の場合においてのみこのような電場の計算を行う。

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

α=21 + sin (π NGMX )11(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 は、格子上に定義される値、例えば電場、数密度等を計算する際に発生する高周波の雑音を除去する目的にも用いられる (11)。Digital

Filter の式は以下のように表される。

φnewi , j ,k←cornerMφ i , j ,k + SΦ iside, j , k + KΦ i , j , k

M + 6 S + 12 K

(3.34)

φ 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)式は次式となる。

14

φ inew
φ i , j ,k +, j ,k ←1 side1
Φ i , j ,k + Φ icorner, j ,k12
48(3.37)

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

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 はイオンプラズマ周波数である。

13

第4章
4.1

解析結果
評価方法

今回、3 つの制御方法についてシミュレーションにより検討する際、推進効率 η と Steering

angle β を磁気ノズルの評価の指標とした。
推進効率 η は以下の式で求めた。

推進効率η =∑ mv (全Z方向運動量)∑ m v (全初期運動量)
z

0

(4.1)

つまり、広がるプラズマ粒子全体の+Z 方向の運動量を増加させると推進効率は高くなる。さらに、広がってゆくプラズマ粒子全ての運動量ベクトルを合成したものを推力ベクトルとし、 Z軸方向とのなす角度を Steering angle β とした。
それぞれの制御方法に対するシミュレーション結果を、推進効率 η と Steering angle β の関係としてまとめ、それについて考察した。

4.2

制御方法 1

4.2.1

計算体系

計算体系を Fig. 4-1 に示す。この制御方法は Orth によって提案され(2)、コイル 1 つのシン
プルな構造の磁気ノズルにおいて Steering angle を得ることを目的としている。
コイルの中心は Z 軸上の原点から-1m に有り、半径は 1m、X-Y 平面に平行である。初期プラズマは初期位置を中心に半径 0.3m の球体内に一様分布している。まず、プラズマ粒子は(tanα, 0, 0)を中心に放射状に初速を持って広がる。計算領域は Z 軸方向に-7m〜7m、X,Y 軸
方向にそれぞれ-6m〜6m の直方体とした。
計算に用いたパラメータを Table 4-1 に示す。主な値は長峯(3)と同じものを用いた。プラズマの初期位置を決める α は 0, 10, 20, 30, 45[deg]の 5 パターンで計算を行った。計算は η、β が収束したとみなせる 40000 time step の 11.08[µ s]まで実行した。

14

Fig. 4-1 制御方法 1 の計算体系

Table 4-1 計算に用いるパラメータ

1.0
Coil radius [m]Coil current [A]
3.57×106-1
Coil position along Z [m](tanα, 0, 0)
Plasma coordinates [m]Plasma radius [m]0.3
4
Plasma energy [MJ]
110
Plasma mass [mg]197
Atomic mass [AMU]79
Atomic numberTime step ∆t [ns]0.277
12.0×12.0×14.0
Calculation region [m]120×120×140
Mesh number
Number of particles
10000015

4.2.2

結果と考察

α = 0[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-2 に示す。0[µ s]における初期プラズマは原点を中心に分布している。その後、プラズマ粒子は Z軸に対称に広がるため Steering angle β はほとんど角度を持たず推移し、11.08[µ s] では、η =0.690、β = -0.1[deg]であった。
α = 10[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-3 に示す。0[µ s]における初期プラズマは X 軸上の x = 0.176[m]を中心に分布している。11.08[µ s]では、η = 0.685、β = 7.2[deg]であった。

α = 20[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-4 に示す。0[µ s]における初期プラズマは X 軸上の x = 0.364[m]を中心に分布している。11.08[µ s]では、η = 0.659、β = 15.6[deg]であった。

α = 30[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-5 に示す。0[µ s]における初期プラズマは X 軸上の x = 0.577[m]を中心に分布している。11.08[µ s]では、η = 0.619、β = 23.8[deg]であった。

α = 45[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-6 に示す。0[µ s]における初期プラズマは X 軸上の x = 1.0[m]を中心に分布している。11.08[µ s]では、η = 0.510、β = 37.2[deg]であった。
以上の結果より、11.08[µ s]における α と β、η の関係を Fig. 4-7 に示す。
この中で、Steering angle β は α に比例して増加した。近似式を以下に示す。

β = 0.8092α
2

(4-1)
2

このときの決定係数 R は 0.9977 であった。決定係数 R は関数式の近似度を見るための尺
度で０と１の間を取り、１に近いほど近似度が高いことを意味する。
この関係は、初期プラズマの位置が Z 軸上から離れるにつれてプラズマに対する磁場の偏りが大きくなり、プラズマが+X 方向により偏向された結果、合成した推力ベクトルの X 方向成分が増加し、さらにプラズマが磁力線に沿った方向に広がっていったためだと考えられる。
また、推進効率 η は α が大きくなると減少している。これは、推力ベクトルの X 成分が大きくなり、その分 Z 成分が小さくなったためである。
また、Fig. 4-8 に示すように α = 45[deg]ではプラズマが宇宙船構造物(Fig. 4-8 の点線で示す)に衝突すると予想されることがわかった。よって、α の限界値は 30[deg]〜45[deg]の間にあると考えられる。
この制御方法では核融合を起こす位置を変えるため、燃料ペレットをレーザー照射地点へ
正確に射出し、さらにレーザーを誘導するミラーを精密に制御する技術が新たに必要になる。

16

Fig. 4-2 α = 0[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]17Fig. 4-3 α = 10[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]18

Fig. 4-4 α = 20[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]19

Fig. 4-5 α = 30[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]20

Fig. 4-6 α = 45[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]21

Thrust efficiency ηSteering angle β [deg]40

Steering angle β [deg]Thrust efficiency η

0.70

0.65

30
0.60

20
0.55

0.50

10

0

10

20

30

40

0

Angle of initial plasma position α [deg]Fig. 4-7 11.08[µ s] における α と β、η の関係

Fig. 4-8 α = 45[deg], 11.08[µ s] におけるプラズマ粒子位置と推力ベクトル
22

4.3

制御方法 2

4.3.1

計算体系

計算体系を Fig. 4-9 に示す。この制御方法では坂口によって最適化された(7)磁気ノズルを基にすることで、高い推進効率と Steering angle を同時に得ることを目的としている。
メインコイルの中心は Z 軸上の、原点から-1m に有り、半径は 1m、X-Y 平面に平行である。
リアコイルの中心は Z 軸上の、原点から-0.5m に有り、半径は 2m、中心軸を極角 γ で傾ける。
初期プラズマは原点を中心に半径 0.3m の球体内に一様分布している。まず、プラズマ粒子
は原点を中心に放射状に初速を持って広がる。計算領域は Z 軸方向に-7m〜7m、X,Y 軸方
向にそれぞれ-6m〜6m の直方体とした。
計算に用いたパラメータを Table 4-2 に示す。リアコイルが傾きを持つ以外は坂口によって最適化されたパラメータ(7)を用いた。今回、リアコイルの傾き角度 γ = 0, 15, 30, 45, 60, 75[deg]の 6 パターンで計算した。
制御方法 1 と同様に、η、β が収束したとみなせる 11.08[µ s]まで計算した。

Fig. 4-9 計算体系

23

Table 4-2 計算に用いたパラメータ1.0
Main coil radius [m]Main coil current [A]
3.57×106-1
Main coil position along Z [m2.0Rear coil radius [m]Rear coil current [A]
1.19*106-0.5
Rear coil position along Z [m](0,0,0)
Plasma coordinates [m]0.3
Plasma radius [m]4Plasma energy [MJ]
110
Plasma mass [mg]0.277Time step ∆t [ns]12.0×12.0×14.0
Calculation region [m]120×120×140
Mesh number
100000
Number of particles
4.3.2

結果と考察

γ = 0[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-10 に示す。プラズマ粒子は Z 軸に対称に広がるため Steering angle β はほとんど角度を持たず推移
し、11.08[µ s]では、η = 0.747、β = -0.2[deg]であった。

γ = 15[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-11に示す。11.08[µ s]では、η = 0.742、β = 4.3[deg]であった。

γ = 30[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-12に示す。11.08[µ s]では、η = 0.731、β = 8.2[deg]であった。

γ = 45[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-13に示す。11.08[µ s]では、η = 0.688、β = 12.3[deg]であった。

γ = 60[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-14に示す。11.08[µ s]では、η = 0.519、β = 17.8[deg]であった。

γ = 75[deg]に対する 0、5.54、11.08[µ s] でのプラズマ粒子位置と推力ベクトルを Fig. 4-15に示す。11.08[µ s]では、η = 0.381、β = 5.3[deg]であった。
以上の結果より、 11.08[µ s]における γ と β 、η の関係を Fig. 4-16 に示す。この中で、

0 ≤ γ ≤ 60 [deg]において Steering angle β は α に比例して増加した。近似式を以下に示す。
24

β = 0.2868α

(4-2)

2

このときの決定係数 R は 0.9954 であった。この関係は、-X 方向に広がるプラズマがリアコ
イルによって遮られ、+X 方向に偏向されたためである。しかし、75[deg]での β はこの式には当てはまらなかった。Fig. 4-15 でわかるようにリアコイルが水平に近く、 γ ≤ 60 [deg] に比べて

-X 方向へ広がるプラズマ粒子が多くなったことで、+X 方向へのプラズマ粒子運動量の偏りが小さくなったためだと考えられる。
また、推進効率 η は α が大きくなると減少している。これは、推力ベクトルの X 成分が大きくなり、その分 Z 成分が小さくなったためである。さらに、γ = 60, 75[deg]ではリアコイルがプラズ
マ粒子の+Z 方向運動を遮るので、急激に減少している。
この制御方法 2 は制御方法 1 に比べ高い推進効率を得られたが、2 つのコイルを用いているため機械構造的に複雑化し、質量も大きくなってしまう。さらに、核融合で発生する中性子の衝突割合が大きくなるので、その分発生する熱量も大きくなり冷却がさらに必要になるなど、マイナスの面も存在する。

25

Fig. 4-10 γ = 0[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]26

Fig. 4-11 γ = 15[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]27

Fig. 4-12 γ = 30[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]28

Fig. 4-13 γ = 45[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]29Fig. 4-14 γ = 60[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]30Fig. 4-15 γ = 75[deg] でのプラズマ粒子位置と推力ベクトル
上から順に 0[µ s]、5.54[µ s]、11.08[µ s]

31

25

Thrust efficiency ηSteering angle β [deg]Steering angle β [deg]Thrust efficiency η

0.8

0.7

20

0.6

15

0.5

10

0.4

5

0.3

0
0

15

30

45

60

Rear coil slope γ [deg]75Fig. 4-16 11.08[µ s] における γ と β、η の関係

4.4

制御方法 3

4.4.1

計算体系

計算体系を Fig. 4-17 に示す。この制御方法は、広がるプラズマの運動量を集中化することで、推進効率を上げながら推力ベクトル制御を目指している。プラズマ運動量集中化の原理
を Fig. 4-18 に示す。核融合燃料ペレットの周りに整形した Moderator（減速材）を配置すると、レーザーを照射して核融合を起こした際にプラズマの運動量が一部に集中し、球形ではない形で広がっていくと考えられる。このペレットを用いると、球形の燃料ペレットに比べて推進効
率が上がり、燃料ペレットの向きを変えて核融合を起こすだけで Steering angle が得られると思われる。
メインコイルの中心は Z 軸上の、原点から-1m に有り、半径は 1m、X-Y 平面に平行である。
初期プラズマは、円錐形の領域内に分布させることによってプラズマの運動量集中化を模擬した。初期プラズマ分布領域の Y = 0 平面における断面を Fig. 4-19 に示す。他の制御方法と比較するため、円錐形の初期プラズマ分布領域の体積を半径 0.3m の球体の体積と同じになるように決定した。初期プラズマ分布領域は、プラズマ粒子に初速度を原点から放射状に与え、初期運動量の総和が 0 となる位置に配置し、推進効率が高くなると予想される形を設定した。そして、この円錐の中心軸を極角 θ で傾けた。計算領域は Z 軸方向に-7m〜7m、X,Y 軸
方向にそれぞれ-6m〜6m の直方体とした。

32

計算に用いたパラメータを Table 4-3 に示す。主な値は長峯(3)と同じものを用いた。円錐中
心軸の極角 θ は 0, 30, 60, 90, 120, 150, 180[deg]の 7 パターンで計算を行った。今回、η、β が収束したとみなせる 11.08[µ s]まで計算した。

Fig. 4-17 計算体系

Fig. 4-18 プラズマ運動量集中化の原理

33

Fig. 4-19 初期プラズマ分布領域

Table 4-3 計算に用いたパラメータ

1.0
Coil radius [m]Coil current [A]
3.57×106-1
Coil position along Z [m]4
Plasma energy [MJ]
110
Plasma mass [mg]197
Atomic mass [AMU]79
Atomic number
0.277
Time step ∆t [ns]12.0×12.0×14.0
Calculation region [m]120×120×140
Mesh number
100000
Number of particles

34

4.4.2

結果と考察

11.08[µ s]における θ と β、η の関係を Fig. 4-20 に示す。θ =0[deg]において、推進効率は0.73 であった。これは、球形の初期プラズマ領域を用いたときの推進効率 0.69 より高い値であり、プラズマの運動量集中化による推進効率の上昇をシミュレートできた。Steering angle β の
大きさは θ = 150[deg]での β = -10.6 が一番大きかったが、このときの推進効率は 0.55 であり、この結果だけでは制御方法 1、制御方法 2 に比べて有効であるとはいえなかった。今回はプラ
ズマ運動量の集中化を初期分布領域の整形によって模擬したが、より実際の現象に近づけるためにプラズマ運動量の集中化を含めたシミュレーションを行わなければならない。その上で、様々なパターンの燃料ペレットについてシミュレーションし、性能向上を目指す必要があると考える。

Thrust efficiency ηSteering angle β [deg]

0.7

10

5

0.65
0
0.6
-5

0.55

-10

0.5

0.45

Steering angle β [deg]Thrust efficiency η

0.75

0

30

60

90

120

150

Plasma cone angle θ [deg]Fig. 4-20 11.08[µ s] における θ と β、η の関係

35

-15
180

第5章

結論

3D hybrid code を用いたシミュレーションにより、3 種類の推進方向制御の可能性を検証した。

X-Y 平面に平行に置いた 1 つのコイルで構成された磁気ノズルを用いて、初期プラズマを(tan α, 0, 0)を中心に分布させた制御方法１においてシミュレーションした。その結果、Steering
angle β とプラズマ初期位置角度 α との関係を以下の式として表せることが分かった。

β = 0.8092α(5-1)

上の式は 0 ≤ γ ≤ 45 [deg]において成り立つ。
また、α = 45[deg]ではプラズマが宇宙船構造物に接触すると予想されるため、α の限界値
は 30[deg]〜45[deg]の間にあると考えられる。
しかし、この制御方法では核融合を起こす位置を変えるため、燃料ペレットをレーザー照射地点へ正確に射出し、さらにレーザーを誘導するミラーを精密に制御する技術が新たに必要になる。
さらに、2 つのコイルで構成された磁気ノズルのリアコイルを角度 γ で傾ける制御方法 2 においてシミュレーションした。その結果、Steering angle β とリアコイル傾き角 γ との関係を以下の式として表せることが分かった。

β = 0.2868α(5-2)

上の式は 0 ≤ γ ≤ 60 [deg]において成り立つ。
しかし、 60 ≤ γ ではリアコイルによってプラズマの+Z 方向運動が妨げられ推進効率が大幅に低下し、有効な推力方向制御ができないことがわかった。
この制御方法 2 は制御方法 1 に比べ高い推進効率を得られたが、2 つのコイルを用いているため機械構造的に複雑化し、質量も大きくなってしまう。さらに、核融合で発生する中性子の衝突断面積が大きくなるので、その分発生する熱量も大きくなり冷却がさらに必要になるな
ど、マイナスの面も存在する。
制御方法 3 では、プラズマの運動量集中化を模擬し、推進効率を上げながら推力ベクトル制御を目指した。今回の計算結果では、制御方法 1、制御方法 2 に比べて有効性は見られなかった。今回はプラズマの運動量集中化を大まかに模擬して計算したので、今後より実際の現象に近づけたシミュレーションを行う必要がある。
今後、Fig. 5-1 のような多数のコイルを繋げた磁気ノズルなど、様々な磁気ノズルで計算を行い、推進効率の向上や Steering angle の増加を目指す。

36

Fig. 5-1 新しい磁気ノズルの概念

37

参考文献
(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) C. D. Orth: VISTA - A Vehicle for Interplanetary Space Transport Application Powered by
Inertial Confinement Fusion, UCRL-LR-110500 (2003)(7) N. Sakaguchi, Y. Kajimura and H. Nakashima: Trans. Japan Soc. Aero. Space Sci., 48 (2005),
180(8) D. S. Harned: J. Comp. Phys., 47 (1982), 452(9) E. J. Horowiz, D. E. Schumaker, and D. V. Anderson: J. Comp. Phys., 84 (1989), 279(10) A. Mankofsky, R. N. Sudan, and J. Denavit, J. Comp. Phys., 70 (1987), 89(11) C. K. Bridsall and A. B. Langdon: "Plasma Physics via Computer Simulation", IOP
Publishing Ltd, (1997)38謝辞
本研究だけでなく 3 年の課題集約演習の時から、直接御指導頂いた中島秀紀教授に深く感謝致します。また卒論の書き方をはじめ、的確な御指導及び貴重な助言をして頂いた渡辺
幸信助教授に深く感謝致します。計算コードをはじめ、様々なことを御指導して頂いた梶村好
宏講師に深く感謝致します。研究室生活の基礎から御指導して頂いた山本直嗣助手に深く感謝致します。快適な研究環境を提供して下さった岡安悦子秘書に深く感謝致します。
多くのアドバイスをして頂いた宮本尚使氏に感謝致します。計算機のことについて膨大な知識を与えてくださった増井博一氏に感謝致します。パソコンをはじめ貴重な研究環境を与えてくださった林真照氏に感謝致します。
スポーツなどの楽しい研究生活をも与えて下さった和泉敬介氏、児玉章裕氏、篠原大介氏、
新屋敷佳祐氏、谷所正彦氏、小川真司氏、加来大輔氏、笠岳幸氏、近岡貴行氏、野田賢治
氏、宗貞伸孝氏、金川隆保氏、近藤慎哉氏、西島康太氏、三好誠氏に感謝致します。特に
篠原大介氏、野田賢治氏には、研究に不可欠なコードや計算機のことなど多大なる知識やア
ドバイスを頂き、深く感謝致します。また、席が隣になったために御迷惑をかけたばかりでなく、忙しい時にも相談にのって下さった児玉章裕氏に謝罪させて頂くと共に、深く感謝致します。
半年の短い間でしたが、課題集約演習で来てくれた江藤修三氏、大村高輝氏、加藤邦男
氏、久保田浩之氏、松田伸夫氏には、頼りない指導に最後まで付きあっていただいて深く感謝致します。
最後に、これまで支えてくださった家族に深く感謝致します。

39