平成１７年度
九州大学大学院 総合理工学府
先端エネルギー理工学専攻

修 士 論 文論文名
３次元ハイブリッドコードを用いた
磁気プラズマセイルの推力特性解析

氏

名

篠原

大介

指導教員名

中島

秀紀

1

教授

目次
第1章

1

序論

1.1

深宇宙探査のための推進システム

1.2

太陽エネルギーを利用する推進システム･･･････････････････････

2

1.3

太陽風-地球磁気圏相互作用 ･･････････････････････････････････

5

1.4

磁気セイルの原理 ･･･････････････････････････････････････････

6

1.5

磁気プラズマセイル(MPS)の原理 ･････････････････････････････

7

1.6

これまでの研究 ･････････････････････････････････････････････

8

1.7

研究目的 ･･･････････････････････････････････････････････････

9

第2章･････････････････････････････

1

3 次元ハイブリッドコード

10

2.1

はじめに ････････････････････････････････････････････････････

10

2.2

基礎方程式 ･･････････････････････････････････････････････････

10

2.3

計算方法 ････････････････････････････････････････････････････

12

2.3.1

プラズマ粒子の運動 ････････････････････････････････････

12

2.3.2

イオン数密度、電流密度････････････････････････････････

14

2.3.3

電磁場 ････････････････････････････････････････････････

15

2.3.4

Digital Filter ･･････････････････････････････････････････

17

2.4

規格化 ･･････････････････････････････････････････････････････

18

2.5

ハイブリッドコードに課せられる制限

19

第3章

･･････････････････････････

磁気セイルの推力解析

20

3.1

計算体系 ････････････････････････････････････････････････････

20

3.2

推力評価方法 ････････････････････････････････････････････････

21

3.3

結果と考察 ･･････････････････････････････････････････････････

23

第4章

3.3.1

各ケースの結果 ････････････････････････････････････････

23

3.3.2

抗力係数 ･･････････････････････････････････････････････

30

磁気インフレーションの解析

31

4.1

計算体系 ････････････････････････････････････････････････････

31

4.2

結果と考察 ･･････････････････････････････････････････････････

32

統合解析への試み

36

5.1

計算体系 ･･･････････････････････････････････････････････････ ･

36

5.2

解析条件 ･･･････････････････････････････････････････････････ ･

36

5.3

結果と考察 ･････････････････････････････････････････････････ ･

38

結論

42

第5章

第6章
参考文献

45

謝辞

47

2

第1章
1.1

序論

深宇宙探査のための推進システム

宇宙開発における深宇宙や惑星探査のミッションはアメリカやロシア、また
近年では日本でも、のぞみ、LUNAR-A、はやぶさなどで知られる探査機を用いて数多く行われてきている[1]。そして、それらによる調査や研究は地球や太陽系、さらには宇宙の誕生と進化を解明するためには必要不可欠なものとされている。このような探査ミッションには、大きなペイロード（探査機総質量から燃料質量を引いたもの）や短期間での運用が求められており、推進システムの性能向上が重要な課題となっている。
現在、探査機のみならず観測衛星や打ち上げロケットにも用いられている化学推進(Fig.1-1)は短時間で大きな推力を発生させることができる。しかしその反面、燃料の消費が非常に早いため、長期的な探査ミッションにおいてはペイロードを大きくできないという欠点が生じる。この問題から、近年では電気推進(Fig.1-2)の研究も盛んに行われており、イオンエンジンなどが既に探査機として実用化されている。イオンエンジンは推進剤を電離させてプラズマをつくり、電極を用いてこれを高速に噴射し推力を得るものである[2]。イオンエンジンに代表される電気推進はその特性から比推力が大きく、化学推進と比較して探査機に大きなペイロードを持たせることが可能となっている。このため電気推進
は今後の宇宙探査の主力となりつつある。だが、電気推進は、電力供給の難しい宇宙空間において高い電力を必要としながらも推力の絶対量が小さいため、
ミッション時間を大幅に短縮することは難しいとされている。

Fig.1-1 化学推進Fig.1-2 電気推進

(Copyright JAXA)(Copyright JAXA)

1

1.2

太陽エネルギーを利用する推進システム

1.1 節での記述から分かるように、今後の深宇宙探査においてより短期間で低コストのミッションを実現するためには、現存の化学推進や電気推進に代わる新たな推進システムを考案する必要がある。
そのような次世代宇宙推進システムとして近年注目されているのが宇宙帆船型の推進システムである。これは、宇宙空間に大きな帆（セイル）を展開し、太陽から放出されるエネルギーを推力へと変換するシステムである。
このようなセイル型の推進システム中で最も古くから考えられていたのがソーラーセイルである(Fig.1-3)。ソーラーセイルは太陽光圧を鏡によって受けて推進する機構であり、他に燃料を必要としないためその比推力は無限大ということになる。2004 年にはソーラーセイルに用いられる薄膜（直径 10ｍ）の展開実
験が実際に宇宙で行われ、実用化が視野に入るところまで進歩している[3]。しかし単位面積あたりの光圧は小さいため、イオンエンジンほどの推力を得るた
めには直径 100ｍ以上の薄膜が必要であると考えられており、技術的に解決すべき課題が残されている。

Fig.1-3 ソーラーセイル(Copyright JAXA)

2

ソーラーセイルとは別に、20 世紀の終盤に Zubrin によって考案されたのが磁気セイル(Magnetic Sail)である。これは太陽風を利用する推進システムであり、探査機から発生させた磁場と太陽風との相互作用によって推力を発生する[4]。
太陽風は探査機の磁場によって偏向を受け、その運動量変化から探査機は推力を得ることができる。磁気セイルはソーラーセイル同様推進剤を必要としない
ため、こちらも深宇宙探査に対して適した推進システムであると考えられた。
しかし、地球近傍での太陽風動圧は太陽光圧と比較してもさらに小さいため、
数 N オーダーの十分な推力を得るためには直径数十 km にも及ぶ巨大な超伝導コイルが必要であると考えられている[4]。このような巨大な構造物を伴う探査機の運用は、現在の技術では難しく、磁気セイルの実現は不可能と考えられてきた。
その後、機械的に広大な磁場を生成するのではなく、探査機からのプラズマ噴射によって小スケールの磁場を広範囲に展開するという考えが注目を浴びた[5]。この方法はワシントン大学の Winglee 博士によって考案され、M2P2(Mini
Magnetospheric Plasma Propulsion)と名付けられており、日本の研究機関ではこれを磁気プラズマセイル(Magneto Plasma Sail : MPS)と呼んでいる(Fig.1-4)。磁気プ
ラズマセイル(MPS)は磁気セイルと異なり、磁場展開に使用するプラズマ噴射のために燃料を必要とするが、太陽風を推進に利用するということから、高いエ
ネルギー効率が考えられている。Table 1-1 に、Zubrin によって提案された磁気
セイル、Winglee によって提案された M2P2 の各パラメータ[6]を示す。また
Fig.1-5 に比推力を横軸に、推力電力比を縦軸にとり各電気推進と MPS の占める領域をプロットしたグラフを示す[7]。ただしここに示す MPS の性能は、理想的な条件のもと簡単な仮定をおいて概算して求められたものである。これらから、
磁気プラズマセイルは深宇宙探査に有効な推進システムとして期待されている。

3

.

1000

Electric Propulsion
Efficiency=100%

Resist Jet

推力電力比︵ mN/kW)

MPS

100

DC arcjetHall thruster
MPD arcjet
Ion Thruster

101000 2000 3000

4000

5000

比推力(s)

Fig.1-5 MPS の予測性能

Fig.1-4 MPS イメージ図(Copyright JAXA)

Table 1-1

磁気セイル、M2P2(MPS)のパラメータ比較

Key Parameters

Magnetic Sail

M2P2

Subsolar Point [km]90

20Coil Diameter [km]31.6

1.E-04

Coil Current [A]

50000

10

Coil Mass [kg]

5000

10

1

1000

Magnetic Moment of Coil [Am ]Coil Surface Magnetic Field [T]

2.E+14

315

8

0.06

Magnetic Field at SpaceCraft [T]

1.E-06

0.06

Injected Plasma Density [m-3]05.E+19

0

0.5

0.02

0.01

Number of Turns in the Coil
2

-1

Fuel Consumption [kg day ]-2

Acceleration [ms ]※Subsolar Point は磁気圏半径で、コイル中心から太陽風と磁場が釣り合う位置
までの距離

4

1.3

太陽風-地球磁気圏相互作用

磁気セイルや磁気プラズマセイルのコンセプトは太陽風と地球磁気圏の相互作用(Fig.1-6)に基づいて考えられている。
太陽コロナのプラズマは 100 万度という高温になるため、静的な釣り合い状
態をとることができず、絶えず外向きに加速されて超音速の流れとなり、惑星間空間に向かって吹き出される[8]。これが太陽風と呼ばれる。地球近傍での太陽風の速度はおよそ 400km/sec、密度はおよそ 5×106m-3 と観測されている[9]。
コロナや太陽風のプラズマの電気伝導度は非常に高く、多くの場合、電気伝導度は無限大として考えられる。そのため、コロナ上空まで延びた磁力線は、太陽風に運ばれて惑星間空間へと広がる。これより、地球近傍にもこのような惑星間磁場が存在するがその磁場強度は平均して数 nT ほどと極めて小さい。コロ
ナの磁力線が開いたところから吹き出してくる太陽風は定常的な惑星間空間の構造をつくるが、磁力線が閉じた領域では、太陽フレアやコロナガス噴出と呼ばれる爆発的な活動現象が発生する。この場合、太陽風は通常よりも高い密度と速度を持つこともある。しかし、この点に関しては理論的にシミュレーショ
ンで予測することは困難なため、本研究では太陽風を定常的なプラズマ流れとして取り扱う。
地球の磁場は、上記のような太陽風に対してあたかも壁のごとく立ちはだか
り、太陽に面した側の地球磁場は圧縮され、流れの圧力とつり合う。地球磁場
は、ほぼ球形の空間に閉じ込められて磁気圏を形成する。太陽風は進路を曲げられ、上下や脇にそれて夜側へと流れていく。このように地球磁気圏は太陽風から粒子の運動量およびエネルギーの注入を受けることになり、これが探査機への推力に利用できるとして考えられるようになった。

Fig.1-6 太陽風-磁気圏相互作用(Copyright JAXA)

5

1.4

磁気セイルの原理

1.2 節で述べたように磁気セイルは太陽風と磁場の相互作用により推力を発生させる推進システムである。1-3 節において、磁場は太陽風から運動量やエネルギーを注入されると説明したが、その具体的なメカニズムにおいては電磁気学的な考え方が存在する。固体壁を帆とする場合、流れの作る圧力分布が直接固体壁に作用し、推力が発生すると考えられる。一方、磁気セイルのように磁場を帆とする場合、磁気圏境界に太陽風プラズマ流れの圧力分布が存在することで、境界面に相応する誘導電流が流れる(Fig.1-7）。この誘導電流によってコイル周辺での磁場(B1,B2)が誘起され、コイル電流(J1,J2)との間で J×B のローレンツ力(F1,F2)が発生する(Fig.1-8）。図においてコイル中心から見て左側（昼側）と右側（夜側）ではコイル電流の向きが反対となるため、互いに逆向きのローレンツ力が発生することが分かる。ただし、左側のほうが磁気圏境界電流に近いことから誘起される磁場が右側よりも大きくなることにより、ローレンツ力は左側のほうが大きくなると考えられる。このローレンツ力の差によってコイルは太陽風方向に力を受けることとなる。これが磁気セイルの推力発生メカニズムである。

Fig.1-7

Fig.1-8

磁気圏境界の誘導電流6コイルに働くローレンツ力

1.5

磁気プラズマセイル（MPS）の原理

MPS は磁気セイルのような太陽風-磁気圏相互作用のメカニズムの他に、探査機からの噴射プラズマによる磁気インフレーションのメカニズムが重要である。
インフレーションには、磁場のプラズマへの凍結という物理現象が利用される。
これは、プラズマが磁力線と共に、あるいは磁力線がプラズマと一緒に動く現象のことを意味する。この現象について以下に説明する。
一般化されたオームの法則とマクスウェル方程式より、次のような磁場の拡散方程式が得られる。
∂B
η 2= ∇ × ( v × B) +
∇ B
∂t
µ0

(1.1)

ここで、 η は電気抵抗率を表し、 µ 0 は真空透磁率を表す。右辺第一項が磁場の凍結に寄与する項である。第二項は磁場の拡散減衰を表す項であり、磁場の凍結を破る働きを持つ。プラズマの抵抗率が小さい場合には、もしプラズマと磁場の動きにずれがあると、プラズマ中の磁束の変化を妨げようとして、大きな渦電流が誘起され、それによって作られる磁場がもとの磁場を変形し、結果として磁力線がプラズマに凍結する。逆に抵抗率が大きい場合、この誘導電流が十分でないために、プラズマと磁力線との間に運動のずれが生じる。(1.1)式の第一項と第二項の比Rm =∇ × ( v × B)η∇ B / µ 02≈µ 0 vLη

(1.2)

を磁気レイノルズ数と呼ぶ。v は流体の代表速度、L は系の特性長を表す。Rm は、
磁力線凍結の程度を見積もるための無次元量として用いられ、 Rm が大きいときほど磁場の凍結が成立する。太陽風プラズマなどは特性長が極めて大きいため、

Rm も地球近辺で約 1015 と非常に大きく、これはいかなる磁場も太陽風に凍結することを表している [10]。
MPS はこのような磁場の凍結を利用したインフレーションにより磁気圏を広げた後、太陽風との相互作用によって推力を発生する推進システムである。相互作用に関しては、磁気セイルと同様のメカニズムである。

7

1.6

これまでの研究2000 年に Winglee によって M2P2 が考案されて以来、アメリカや日本において、様々な研究が行われている。MPS は実験的なアプローチが難しいため、現在のところは数値解析による研究結果が多く発表されている。
MPS のアイデアを提唱した Winglee は電磁流体 (MHD)モデルを用い、磁場インフレーション、太陽風 -磁場相互作用双方を含めた全系での解析を行い、その推力評価によって MPS の実現性を示した [5]。 Nishida は MHD シミュレーショ
ンを行い、磁気セイルの太陽風 -磁場相互作用による推力発生メカニズムを明らかにし、更に推力ベクトル制御の可能性を示した [11]。 Asahi は MHD 解析によって、 MPS の磁気インフレーションの解析や性能評価を行った [12]。最近の研究としては、 Nagata が MHD を用いて、イオンスラスタと比較した性能評価を行っている [13]。このように磁気セイルや MPS の研究は、電磁流体解析により、実現性を議論する踏み込んだ内容まで進められている。
ここで一つ考慮しなければならないことは、MHD が電磁流体近似を仮定した解析であり、rL / L <1 の範囲の現象でなければ成立しないということである。ここで rL はイオンのラーマ半径、 L は系の代表長を表す（ここで言う代表長は磁気圏の半径と定義する）。太陽風のラーマ半径 rLs は、磁気セイルや MPS の磁気圏境界付近において約 100km であることが計算により分かる。この値から、地球磁気圏のような広大な代表長を持つ系ならば流体近似は十分に成立することが分かる。だがそれに対し、MPS では、直径数メートルほどのコイルで生成可能な小規模の磁場を太陽風イオンラーマ半径である 100km 以上に展開することは難しいと考えられ、 MPS の太陽風 - 磁気圏相互作用メカニズムにおいては

rLs / L >1 の粒子的効果が支配的となると考えられる。また、磁気インフレーションメカニズムにおいても、噴射プラズマのラーマ半径 rLi はコイルから離れた
弱磁場領域では系よりも大きくなってしまい粒子性が強まることが考えられる。
そのような有限ラーマ半径効果を考慮するため、粒子コードによる研究も必要である。近年、Fujita は太陽風と磁場の干渉を、イオン、電子ともに粒子として扱う Full-PIC コードによって解析し、荷電分離が MPS の推力発生において重要な役割を果たすことを明らかにした [14]。また同氏はイオンを粒子、電子を流体として扱うハイブリッドコードによって磁気セイルの推力特性評価も行っている [15]。Khazanov は磁気セイルについてハイブリッドコードにより解析を行い、粒子運動が支配的な領域では大きな推力が見込めないことを指摘している [16]。
Omidi は相互作用に加え、インフレーションを 2 次元のハイブリッドコードによって解析し、噴射プラズマによる磁場の展開が難しいことを示した [17]。
8

1.7

研究目的

前節で示したように磁気セイルや MPS の研究において、粒子運動を考慮した
数値解析は実現可能性を左右する重要なものだと考えられる。そのため本研究では、 3 次元のハイブリッドコードを用いて粒子性を考慮した磁気セイル（太陽風 -磁気圏相互作用のみ）や磁気インフレーションの解析を行い、最終的にそれら双方のメカニズムを同時に解析して、MPS における推力特性解析を実施す
ることを目的としている。ハイブリッドコードは、粒子コードの中でも巨視的な系で有効とされる数値解析コードである。前節で紹介したように、ハイブリッドコードを用いた解析の事例は過去にもある。しかし、本論文で、最終的な目的としている太陽風 - 磁気圏相互作用と磁気インフレーション双方の同時解
析は Omidi が 2 次元のハイブリッドコードで行っているのみであり、プラズマ噴射の手法についても十分な検討がされていない。これらの背景を考慮し、本研究では以下のような解析を行う。
相互作用メカニズムのみが考慮される磁気セイルにおいては、ラーマ半径代
表長比 rLs / L（ rLs：太陽風ラーマ半径）に対する推力依存性の調査を目的とする。
本研究では、代表長を変化させた 5 ケースでの抗力係数を算出して評価を行う。
また、Fujita はこれと同様の解析を行っているため、その結果と比較することで検証結果の妥当性を確認する。
磁気インフレーションでは、プラズマ噴射により磁場を遠方へ展開できることを示す必要がある。通常、ダイポール磁場は B ∝ r −3 で減衰するが、本研究ではこの減衰勾配をプラズマ噴射によりどの程度減らせるかを確認し、Asahi の結果と比較する。
MPS 解析の足がかりとして、本コードにおいて現実的な時間で計算可能なパラメータを用いて相互作用と磁気インフレーション双方のメカニズムを同時に解析する。本研究では、MPS で重要パラメータとなる噴射プラズマ β 値、 rLi / L（ rLi ：噴射プラズマラーマ半径）の依存性を調査することを目的とする。
本論文の構成を以下に示す。
第 2章

解析に用いる 3 次元ハイブリッドコードに関して述べる。

第 3章

磁気セイルの解析方法、解析結果と考察を述べる。

第 4章

磁気インフレーションの解析方法、解析結果と考察を述べる。

第 5章

太陽風 - 磁気圏相互作用と磁気インフレーションの統合解析について
解析方法、解析結果と考察を述べる。

第 6章

総括を述べる。
9

第2章
2.1

3 次元ハイブリッドコード

はじめに

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

2.2

基礎方程式

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

(2.1)

dx i
= vi
dt

(2.2)

である。ここで v i はイオンの速度、 Z はイオンの価数、 e は電気素量、 x i はイ
オンの位置、 mi はイオンの質量、また、 E 、 B はそれぞれ電場、磁場を表している。

10

流体として扱う電子の運動方程式を電場の解法に用いる。電子の運動方程式は以下の式で表される。
ne me

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

(2.3)

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

である。ここで Te は電子温度である。(2.4)式を (2.3)式に代入して Pe を消去する。
0 = −ene (E + v e × B ) − ∇(neTe )(2.5)
荷電準中性条件および、ダーウィン近似を用いたアンペアの法則により∇ × B p = µ 0 (J e + J i )(2.6)

が得られる。但し B p は、プラズマにより生成される磁場で、以後特に断りのな
い限り、磁場は外部磁場とプラズマにより生成される磁場とのベクトル和であるとする。また J i はイオン電流密度、 J e は電子電流密度、 µ 0 は真空中の透磁率である。 J i はイオンの数密度と電流密度から計算することができる。 J e は電子電流密度の定義より、
J e = −ene v e(2.7)

である。 (2.6)、 (2.7)式を用いて J e を消去し v e について求めると、
ve = −⎞1 ⎛ 1⎜⎜ ∇ × B p − J i ⎟⎟ene ⎝ µ 0⎠(2.8)

となる。これを (2.5) 式に代入して電場 E について整理する。荷電準中性条件より ene = Zeni を考慮し、
E=

⎫
1 ⎧1
⎨ (∇ × B p ) × B − J i × B − ZTe ∇ni ⎬Zeni ⎩ µ 0
⎭

(2.9)

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

11

(2.10)

計算方法

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

外部磁場および、イオンの粒子位置、速度の初期設定を行う。

②

格子点上でのイオンの数密度、電流密度を計算する。

③(2.9)式により、電場を求める。

④①で設定した外部磁場と、③で求めたプラズマによる電磁場から、(2.1)式、(2.2)式を用いて新たなイオン粒子の位置、速度を求める。

⑤(2.10)式により磁場を求める。

②〜⑤を繰り返し、プラズマの時間的変化を追う。
3 次元ハイブリッドコードに対する各方程式の差分式を以下の小節で示す。

2.3.1

プラズマ粒子の運動

本コードでは、いくつかのイオンを超粒子と呼ばれる集団として扱っている。
つまり、超粒子を個別に運動させることでプラズマ粒子の運動を追う。 (2.2) 式に、 leap-frog 法 [19]による時間中心差分を施すと、
x in +1 − x in
= v in +1 / 2∆t(2.11)となる。ここで各変数の上付き添え字は時間を表し、下付き添え字は i 番目のイオンを表す。また、 ∆t は時間幅を表す。
(2.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
⎭

(

)

(2.12)となる。n − 1 / 2 はある時刻 n から ∆t / 2 だけ以前の値であることを示す。同様に、
n + 1 / 2 は ∆t / 2 だけ進んだ値であることを示す。上付き添え字が n 、 n − 1 / 2 である変数は既知量である。
プラズマ粒子の運動を求めるためには、まず (2.12) 式において、 v in +1 / 2 を左辺にまとめて、未知量である v in +1 / 2 を解く。次に、求めた v in +1 / 2 を (2.11)式の右辺に代入し、新たな粒子の位置 x in +1 を求めることが可能となる。

12

(2.12)式を変形すると、Av n +1 / 2 = S(2.13)となる。 α 、 β 、 γ 、 δ を次のように定義する。

α=

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

(2.14)(2.15)

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

v n +1 / 2 = A −1S
ここで、 A −1 は以下のように得られる。
⎛ 1 + α 2 αβ + γ αγ − β ⎞⎟1⎜
A −1 = ⎜ αβ − γ 1 + β 2 βγ + α ⎟δ⎜2 ⎟⎝αγ + β βγ − α 1 + γ ⎠13(2.16)(2.17)

2.3.2

イオン数密度、電流密度

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

δni (i, j , k ) =V2
NV(2.18)

で与えられる。すべての粒子についての和を δni についてとれば、ある格子点上の ni を定義することが出来る。
イオン電流密度 J i は、数密度 ni を用いて、
J i = Zeni v iより求めることが出来る。

Fig.2-1 PIC 法

14(2.19)

2.3.3(1)

電磁場

プラズマ領域での計算

電磁場を求めるためには、 (2.9) 式、 (2.10) 式を用いる。 (2.10) 式を片側差分して、 B n +1 / 2 だけを左辺に残すと、

B n +1 / 2 = B n −∆t∇ × En2(2.20)となる。また、 (2.9)式は、
En =

⎫
T
1⎧ 1
1 n
∇ × B np × B n −J i × B n − e ∇nin ⎬⎨ni ⎩ µ 0 Ze
Ze
e
⎭

(

)

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

+1
E npred= −E n + 2E n +1 / 2+1
B npred= B n +1 / 2 −∆t+1∇ × E npred2(2.22)(2.23)+1+1
ここで、下付き添え字 pred は予測値であることを表している。 E npred、 B npred
を求
3/ 2
3/ 2
3/ 2
めた後、それらによるイオンの運動を解き、 v in +pred、 nin +pred、 J in +predを求める。次+3 / 2に、 B npred
を次式により求める。

+3 / 2+1= B npred−
B npred∆t+1∇ × E npred2(2.24)+3 / 2
そして、(2.21)式で E n を求めた時と同様にして、 E npredを求める。最後に次式に示す線形内挿から次のタイムステップの値、 E n +1 、 B n +1 を求める。

E n +1 =1 n +1 / 2 1 n +3 / 2E+ E pred2
2

B n +1 = B n +1 / 2 −∆t∇ × E n +1
2

(2.25)(2.26)

これらの値 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 タイムステップの値を再度予測値として用いる。具15

+1
体的には、 (2.26) 式により求められた B n +1 を次の予測値 B npredとし、 (2.22) 式から(2.26)式までを次式で定義する条件を満足するまで、繰り返し計算を行う [20]。

max i , j ,k Bin, +j ,1k( q +1) − Bin, +j ,1k( q )max i , j ,k Bin, +j ,1k( q )≤ε(2.27)

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

(2)

真空領域での電場の計算

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

α=

2

1 + sin (π NGMX )(2.28)

ここで 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(2.29)+ (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)

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

16

2.3.4

Digital Filter

PIC 法で格子点上の量を決める場合、例えば粒子の位置を格子上に割り振る

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

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

M + 6 S + 12 K

(2.31)

φ 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(2.32)

Φ 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(2.33)+ φ i −1, j +1,k + φ i +1, j +1,k + φ i +1, j −1,kここで、 M 、 S 、 K の値をそれぞれ 12、 4、 1 とすると (2.31)式は次式となる。

14

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

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

17

2.4

規格化

コード内で扱われる物理量は全て、規格化された値となっている。規格化のために用いられるパラメータは、初期磁場 B0 [T]、サイクロトロン周波数 ω ci [1/s]、
光速 c [m/s]、タイムステップとなる微少時間幅 ∆t [s]である。ここで、サイクロ
トロン周波数、

ω ci =

ZeB0
mi(2.35)

は、イオン質量 mi [kg]、初期磁場 B0 [T]、価数 Z 、電荷素量 e [C]からなる。
規格化する物理量として距離、速度、電磁場、電流密度、および粒子密度が挙げられる。これらは、以下のように規格化される。
∆tω ci
∆tZe
x
~
= E⋅
⇒
E⋅
⇒
x
mi ⋅ cB0 ⋅ cc∆t
c 2 ∆t 3
v
⇒
⇒
J⋅
v~
Ze
c
3
∆tω ci∆tZe~ni (c∆t )
⇒
= B⋅
⇒
B
B⋅
mi
B0

~
E
~
J
n~i

ここで、文字の頭についている波線 (ティルディ )は規格化された量であることを
意味する。

18

2.5

ハイブリッドコードに課せられる制限

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

(2.36)

これは、 ∆t , ∆x → 0 で解が収束する条件である。 V A はアルフベン速度であり、B0(2.37)VA ≡µ 0 ni miで定義される。また、付加的な数値的安定条件は
ω ci ∆t < 0.2(2.38)である。ここで、 ω ci はイオンサイクロトロン周波数である。
メッシュ幅の最小値を求める条件として∆x > c ω pi

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

ωD が存在し、時間刻み ∆t はこの逆数より小さい値としなければ、解は発散するとされている [22]。

ωD =(π / ∆x) 2 Bµ 0 ne19(2.40)

第3章
3.1

磁気セイル解析

計算体系

まずは MPS の太陽風 -磁気圏相互作用のみに着目して推力特性解析を行うた
め、磁気インフレーションを伴わない磁気セイルの解析を行った。それに関する計算方法、計算結果をこの章で解説する。磁気セイル解析の計算体系を Fig.3-1に示す。計算体系は右手系直行座標系として、系の中心にダイポール磁場を配置し、コイルは YZ 平面にあるものとした。実際の地球付近の宇宙空間には、
惑星間磁場 (IMF)が存在するが数 nT と極めて小さいために、本解析では 0 と置いた。太陽風を構成する超粒子（代表粒子）は初期状態で全体系内にランダムに配置され、それら全てに対し Z 軸正の方向への初期速度を与えている。超粒子は Z 軸負の境界から継続的に入射され、境界外側に出ると計算対象外となる。
境界での電磁場はその勾配を一つ内側のメッシュに対し 0 とすることで与えられる。
Table 3-1 に計算で用いたパラメータを示す。超粒子として用いるイオンには、

太陽風の中で 96%を占める水素イオンを採用した。太陽風の密度、速度、電子
温度は地球近傍の値を用いて代表値とした。ω ci はイオンサイクロトロン周波数
で、ここに示しているのは磁場強度が最も大きい位置における値である。

L

Fig.3-1 磁気セイル計算体系

20

Table 3-1 磁気セイル計算パラメータ
Ion
Ion Number Density [m-3]Ion Initial Velocity [km/s]

5×106
400

Number of Super Particle

3000000

Weight of Super Particle

1017〜1021
10Electron Temperature [eV]ΔtTotal Time Step [step]

3.2

Hydrogen＜0.01ωci-1
50000step

推力評価方法

第 1 章で示したように、本研究の目的は有限ラーマ半径効果を考慮した推力評価を行うことである。そのため、ラーマ半径代表長比 rLs / L（ rLs は太陽風のラーマ半径）を変化させた場合の推力への依存性を調査することとした。本解析においては、代表長を 800km から 40km まで、つまり rLs / L が１より小さい領域
から１より大きい領域での 5 ケースを例にとって計算を行った。各ケースのパラメータを Table 3-2 に示す。ここで代表長 L は初期のダイポール磁気圏の半径を指す。磁気圏半径は、太陽風の動圧と磁気圧が電磁流体的に釣り合う位置からコイル中心までの距離と定義した。釣り合いを考えるにあたって以下のようなプラズマ β 値が用いられる。

1 2
ρv
β= 2 2
B
2µ 0

(3.1)

Table 3-1 に示した太陽風のパラメータを用い、 β を 1 として B を求めるとおよそ 41nT となる。つまり、磁気圏半径 L はコイル中心から磁場強度が 41nT になる位置までの距離ということになる。本解析では、コイル磁場の強度を変える
ことによって磁気圏代表長を変化させ、 L をパラメトリックに変更している。
ラーマ半径 rL はこの磁場強度と太陽風代表速度 400km/s から、磁気圏周辺でおよそ 100km と計算することができ、今回はこの値を参考に rLs / L を求めた。
21

Table 3-2 磁気セイル解析ケース
rLs/L

System Size(km)

21

0.13

4000

250

10.5

0.25

2000

200

125

5.25

0.50

1000

4

80

50

2.1

1.25

400

5

40

25

1.05

2.50

200

Case

L (km)

Coil Radius(km) Coil Current[kA]

1

800

500

2

400

3

推力特性を調査するための値として、次のような抗力係数（無次元推力値）C d
を用いた。

Cd =

F
1 2
ρv S
2

(3.2)

分母は太陽風の動圧を表す。 S は磁気圏断面積を示し S = πL2 である。抗力係数は、推力を規格化（無次元化）することにより太陽風 -磁気圏の相互作用の強さそのものを計ることができる値である。ここで仮に、 C d =0.5 とすると F =1Nを得る場合、その代表長 L は 30km となる。
推力 F の算出方法について説明する。1.4 節で示したように、本来、磁気セイ
ルや MPS は電磁気学的相互作用によってコイルがローレンツ力を受けて加速する。しかし、最近の研究によりそのローレンツ力と磁気圏周辺の太陽風の運動量変化は近似できるということが示されている [10]。本解析においては、コイ
ル半径に対して十分なメッシュ数を用いていないため、上記のような電磁気学的な力を直接算出することは困難である。そこで、計算体系内の太陽風の運動量変化を求めることで、推力を見積もることとしている。

22

3.3

結果と考察

3.3.1

各ケースの結果

各ケースにおける計算結果を Fig.3-2〜 3-6 に示す。各図はそれぞれ X-Z 面における、(a)太陽風イオンの粒子分布と磁場ベクトル、(b)密度分布、(c)磁場強度分布であり、推力が定常状態に落ち着いた時間における図を表示した。 Fig.3-2における矢印は太陽風粒子の初期速度方向である。これらの図に関して以下に説明する。
Fig.3-2 は、Table 3-2 で示した Case1 における結果である。このケースでは rLs / L

が 0.13 と小さく MHD 近似が十分に成り立つ領域であるといえる。
Fig.3-2(a)、(b)の結果から分かるように、太陽風粒子が体系中心のダイポール磁場に強く影響を受けて、軌道が大きく変化していることを確認できた。イオンが反射される領域では、密度が高くなっており、太陽風密度の約 5 倍という結果が得られた。これは、理論値 [23]ともよい一致を得られる結果である。また、(c)の磁場も後方に大きく引き伸ばされており、総合的に強い相互作用が起きていることを見て取ることができた。 Case2(Fig.3-3)においても MHD 条件が成り立ち、各分布図も Case1 とほぼ変化がないことが分かった。 Case3(Fig.3-4)はラーマ半径と磁気圏半径の比が 0.5 の結果である。 Fig.3-4(a)からは、粒子が磁気圏内に入り込んでしまい、Case1,2 と比較して相互作用が弱くなっている様子を見て取れた。Case4(Fig.3-5)では rLs / L が１より大きくなり、粒子の磁場に対する偏向は Case3 よりもさらに小さくなった。 Case5(Fig.3-6)ではそれが更に顕著に表れており、粒子の偏向はこの図からは視認しにくいほどまで弱くなり、磁場に対する影響も小さく、太陽風と電磁流体的に釣り合うとされる 41nT の位置の変化がほとんど無いことが確認できた。

23

(a)

粒子分布

(b)

数密度分布

(c)

磁場強度分布Fig.3-2

Case1 結果 (t=12.0s)

24

(a)

粒子分布

(b)

密度分布

(c)

磁場強度分布Fig.3-3

Case2 結果 (t=6.0s)

25

(a)

粒子分布

(b)

密度分布

(a)

磁場強度分布Fig.3-4

Case3 結果 (t=3.0s)

26

(a)

粒子分布

(b)

密度分布

(c)
Fig.3-5

磁場強度分布
Case4 結果 (t=1.2s)

27

(a)

粒子分布

(b)

密度分布

(c)

磁場分布

Fig.3-6

Case5 結果 (t=0.6s)

28

Fig.3-7 に体系中心付近でのイオンの粒子軌道を示す。(a)が Case1( rLs / L =0.13)の場合、 (b)が Case5( rLs / L =2.5)の場合である。 (a)に示すように、ラーマ半径が磁気圏代表長と比較して小さい場合、プラズマが磁場にトラップされ、磁気圏境界での粒子の偏向角が大きく相互作用が強い様子を確認できる。それに対し、
磁気圏代表長が小さい (b)の場合、相互作用が弱くなり偏向角も小さくなっていることを見て取ることができる。これは、粒子の有限ラーマ半径効果が支配的となることで、電磁流体的に釣り合うとされる磁気圏境界 (41nT) の位置でトラップされなくなるからである。この場合、粒子は磁気圏境界のさらに内側の強磁場領域で初めて相互作用する。この作用が推力発生につながると考えられる。

(a)

Case1

(b)

Case5

Fig.3-7粒子軌道

29

3.3.2

抗力係数

Fig.3-8 に各ケースの抗力係数を示す。横軸は磁気圏代表長、縦軸は無次元抗力係数である。曲線は、 Fujita の結果である [15]。鎖線は rLs / L =1 を表し、その右側の領域では電磁流体近似が成り立ち、左側では粒子的効果が支配的となる。
これから分かるように、代表長が 400km 以上では抗力係数は高く一定となり、代表長が小さくなるほど太陽風 -磁気圏相互作用は弱くなり、抗力係数は小さくなった。このように、粒子性がより強くなる領域ほど、太陽風の運動量輸送を効率的に受け取れなくなる。しかし、Case5(L=40km, C d =0.89)での推力は 3N と求められた。実際にイオンエンジンを搭載した探査機「はやぶさ」の推力が 30mNほどであるため [3]、この推力値は探査機に用いるシステムとして十分な値であると言える。以上から、磁気インフレーションによってこの程度のサイズの磁気圏を実現できるかが課題となる。

rLs / L = 1

MHD 領域

Drag Coefficient

4
3.5
3

Case1

Case2

2.5

Case3

2

Case4

1.5
1

Present Research

Case5

0.5

Fujita

0
0

200

400

600

800

Magnetsphere Size : L [km]Fig.3-8

抗力係数

30

1000

第4章

磁気インフレーション解析

4.1 計算体系
この章では、磁気インフレーションのみを考慮した解析に関して、計算方法、計算結果を述べる。計算体系を Fig.4-1 に示す。本解析では、探査機をコイル半径と同様の半径をもつ球体と仮定する。プラズマの噴射方法としては、その球体の表面から、数密度 1.0×1016m-3 のプラズマを 1 タイムステップごとに、初期
速度 × 時間刻み ∆t の幅に定義し、球面に対し外向き法線方向に速度を与え、噴射した。本計算では、初期ダイポール磁場の 41nT（太陽風と磁場が電磁流体的に釣り合う）の位置が、初期と比較してどの程度展開されるのかを調査し、そのときの磁場の分布を確認する。また、一般的にダイポール磁場は B ∝ r −3 で減衰するが、この減衰勾配をプラズマ噴射によりどの程度減らせるかを確認する。
本計算で用いた解析パラメータを Table 4-1 に示す。
本来、磁気インフレーションは磁場の凍結現象を応用したものであるため、(1-2)式で示した磁気レイノルズ数 Re を考慮する必要がある。だが、本解析に用いたパラメータでは Re がおよそ 105 と大きい。よってプラズマと磁場は完全に凍結すると仮定している。具体的には、コード上で用いる運動方程式に電気抵抗の項を考慮せずに計算している。

Fig.4-1 磁気インフレーション計算体系

31

Table 4-1 磁気インフレーション解析パラメータ
Ion

Argon

Atomic Mass

40

Atomic Valence1-3

Ion Number Density [m ]Ion Initial Velocity [km/s]

1×1016

Number of Super Particle

1450000

Weight of Super Particle

1022

Electron Temperature [eV]10

System Size [km]160

Coil Radius [m]2Coil Current [A]

7.1×104

Magnetic Field around coil [µT]
Δt [s]

Total Time Step [step]40

5.7
9.0×10-43630

4.2 結果と考察
解析結果を Fig.4-2、Fig.4-3、Fig4-4 に示す。Fig.4-2 は t =0s、Fig.4-3 は t =0.32s、
Fig.4-4 は t =0.64s の時の磁場強度分布である。 (a)は y=0 での ZX 平面の磁場強度分布であり、 (b)は (a) における点線部の磁場強度分布と B =41nT（直線）を表
す。 (a)の線は内側から B =100nT,41nT,10nT,1nT の磁場の等高線である。
各図の (a)を見ると、磁場の等高線が時間の経過ともに全体として広がっていることが分かる。これより、プラズマ噴射による磁場展開現象が確認されたことが分かる。
各図の (b)を見ると t =0s では r =15km 付近にあった B =41nT の磁場が、t = 0.64 s
では r = 50km 付近に展開されている。これから、 B =41nT の磁場を遠方まで展開するという今回の目的は達成された。
しかし、コイル半径は現実的な MPS のシステムで想定されている長さ (〜数十 m)よりも非常に大きい値である。これにはコード上の制限からの理由があり、詳細は 6 章で説明する。

32

(a)

(b)
Fig.4-2

磁場強度分布 ( t =0s)

(a)

(b)
Fig.4-3

磁場強度分布 ( t =0.32s)

33

(a)

(b)
Fig.4-4

磁場強度分布 ( t =0.64s)34

t =2.7s での磁束密度分布と B ∝ r −1.51 (Asahi の MHD による解析結果 )のラインを比較したものを、Fig.4-5 に示す。t =2.7s ではコイル中心から 100km 離れた位置でさえも B =41nT 以上の値となっていた。
Fig.4-5 では r =40km 以上の磁場がスムーズに展開されているのに対して、

r =40km 内の磁場は大きく変動していた。この原因として、電磁波の励起によるもの、空間の精度が粗いこと、超粒子の重みが大きいことから一個の粒子の動きによる磁場への影響が大きいことが挙げられるが、今後検討し、明らかにする必要がある。

Fig.4-5

磁場強度分布 ( t =2.7s)

35

第5章
5.1

統合解析への試み

計算体系

この章では、太陽風 -磁気圏相互作用と磁気インフレーションのメカニズム双方の同時解析について、研究方法、研究結果を紹介する。 4 章でも触れたが、本ハイブリッドコードの制限上、実際の MPS モデルでの計算が現実的な時間では、困難であり、太陽風のような低密度領域では実スケールと比較して大きな体系しか計算することができない（詳細は 6 章）。しかし、 MPS の性質を議論する中で重要となる無次元パラメータの噴射プラズマ β 値や rL / L の値を現実的なパラメータに合わせることで、MPS のリアルスケールでの性質も予測できると考え、本解析を実施した。 Fig.5-1 に計算体系を示す。今回の解析では磁気セイルの計算と比較するため、磁気セイル解析の体系にプラズマ噴射を加えるという条件で行った。この解析では太陽風プラズマと噴射プラズマに関して、それぞれの超粒子（代表粒子）を用いた。噴射プラズマ以外の解析パラメータは Table 3-1 に示すものと同様である。噴射方法は磁気インフレーション解析と同様で、噴射プラズマ代表粒子にはアルゴンイオン（原子量 40、価数 1）を用い、初期速度は 40km/s とした。

5.2

計算条件

太陽風 -磁気圏相互作用においては、磁気圏代表長が小さくなるにつれ粒子性が支配的となることが明らかになったが、磁気インフレーションにおいても粒子性、流体性の両方を考慮することが重要である。例えば、コイル半径を 10m、コイル表面磁場 0.01T を仮定し、一般的なイオンエンジン（粒子：アルゴン）の排気速度 4km/s を用いて考えると、コイル周辺の強磁場領域 (0.01T)ではラー
マ半径はおよそ 0.2m となり流体近似が成り立つが、噴射プラズマがコイル中心から 100m の位置（磁場が距離の 3 乗で減衰すれば 10-5T）に達するとラーマ半
径はおよそ 160m となり粒子性が強くなってしまう。よって流体条件、粒子条
件のいずれに対してもインフレーションが実現できることが望ましいとされる。
以上の理由より、本解析においてはラーマ半径代表長比 rLi / L を 0.1〜 10 に変えることにより、その依存性を調査した。ただし、ここで言う rLi は噴射プラズマの代表ラーマ半径であり、太陽風と磁気圏の境界付近 (41nT) での値である。また、 L は磁気セイル同様、計算初期状態の磁気圏代表長を定義している。
MPS に関しては、上記の rLi / L 以外にも重要な技術課題として推力電力比を検36

討しなければならない。推力電力比は単位電力に対する推力の比のことである。
つまり、推進システムとして実用的であるためには、噴射プラズマのエネルギー（に要する電力）を極力小さくして高推力を得なければならないのである。
よって本研究では、噴射プラズマの β i 値を変化させた場合の解析も行った。 β i
値の定義式は、 (3.1) 式と同じであり、分子が噴射プラズマの動圧で、分母が噴
射プラズマ初期位置における磁気圧となる。尚、太陽風 -磁気圏間の β 値と区別するため、今後 β i とする。
今回の解析では、 rLi / L を 0.1,10 とし、 β i の値を 0.1,1,10,100 とした計 8 ケースの計算を行った。rLi / L に関しては磁気圏代表長を変化させ、 β i に関しては噴射プラズマ密度を変化させた。 Table 5-1 に解析パラメータを示す。

Fig.5-1

MPS 計算体系

Table 5-1

MPS 解析ケース

Case1

Case2

0.1

1

rLi/L

Case3

Case4

Case5

Case6

10

100

0.1

1

0.1

βi

Case7

Case8

10

100

10

Initial Magnetosphere Radius (km)

4000

40

Coil Radius (km)

2500

25

Coil Current (kA)

10.5

0.105

Δt (s)

1.27E-03

1.27E-05
50000

Total Time Step (step)-3

Injected Plasma Density (m )5.0E+06 5.0E+07 5.0E+08 5.0E+09 5.0E+06 5.0E+07 5.0E+08 5.0E+09

37

5.3

結果と考察

Fig.5-2 に rLi / L =0.1、 Fig.5-3 に rLi / L =10 における粒子分布ならびに磁場強度分布を示し、噴射開始からの時刻を示す。それぞれについて β i 値が 0.1 と 100の際の結果を示した。粒子分布図に関しては赤い点が太陽風を表し、青い点が噴射プラズマを表す。Fig.5-4 に各ケースの推力を示す。ここで、 F0 は同じ代表長における磁気セイル（プラズマ噴射無しの場合）の推力である。
本解析では、β i <1 のケースでも磁場が広がり、推力増分がある結果となった。
これは第４章で述べたように、磁気レイノルズ数が無限大という仮定のもとで解析を行っているため、磁場がプラズマに対し完全に凍結するからである。ただし、β i 値によってその影響力に変化があることが Fig.5-4 から分かる。これは、

β i が小さいと粒子の運動エネルギーが小さく、磁場に与える影響が弱くなるためである。 Fig.5-2(a),5-3(a) の粒子分布から、 β i が小さいケースでは噴射プラズ
マ粒子の軌道が磁場によって影響を受けている様子が見て取れる。
Fig.5-2 は rLi / L =0.1 であり、 rLs / L <1 であるため、すべてのプラズマは流体的にふるまうパラメータにおける計算である。 rLs / L <1 であるため、 β i を変えて
噴射プラズマのエネルギーを大きくし磁場拡大を促進したとしても推力の増分は小さい。このことは Fig.5-4 での rLi / L <1 の結果からも判断できる。一方、
Fig.5-3 は rLi / L =10 で rLs / L >1 であり、粒子効果が顕著に表れるパラメータによる 計 算 で あ る 。 こ の 場 合 、 βi が 強 い と き 磁 気 圏 が 拡 大 さ れ て い る 様 子 がFig.5-3(b)から分かる。 rLs / L >1 という条件下の計算では、わずかな磁気圏拡大でも推力の増分は大きいということが Fig.3-8 から分かる。そして Fig.5-4 でのrLi / L >1 の結果もそのことを示している。ただし、本解析では rLi / L を変化させる際に L を変化させて計算を行っている。これは体系自体の大きさを変化させていることに相当し、今後、同じ体系長での rLi / L の依存性を調査するために磁場の値を変化させて解析する（ rLi を変化） ことも必要である。

38

(a) Case1 ( β i =0.1)(b) Case4 ( β i =100)
Fig.5-2

rLi / L =0.1 での粒子と磁場強度の分布 (t=60s)
39

(a) Case5 ( β i =0.1)(b) Case8 ( β i =100)
Fig.5-3

rLi / L =10 での粒子と磁場強度の分布 (t=0.6s)
40

2
rLi/L=0.1
rLi/L=10

1.8

F/F0

1.6
1.4
1.2
1
0.1

1

10

Injected Plasma β i
Fig.5-4

rLi / L と β i に対する推力の変化

41

100

第6章

結論

本研究では、太陽風と磁場の相互作用を利用する推進システムである磁気セイル解析や磁気インフレーションの解析、また双方のメカニズムの統合解析を
3 次元ハイブリッドコードを用いて行った。それぞれの結果は以下に示す通りである。

・ 磁気セイル解析においては、太陽風 - 磁気圏間の相互作用における有限ラー
マ半径効果を調査するために、磁気圏境界における太陽風ラーマ半径と磁気圏半径の比 rLs / L を変化させた場合の抗力係数を算出した。その結果として、流体近似が成り立つ領域 ( rLs < L )では抗力係数 C d は高く一定であり、粒子性が支配的となる領域 ( rLs > L )では磁気圏が小さくなるほど太陽風との相互作
用が弱くなり C d 値は小さくなった。磁気圏半径 40km で得られる推力はおよそ 3N である。
・ 磁気インフレーションの解析では、太陽風が吹かない条件下でのプラズマ噴
射による磁場の展開を確認することができた。結果として通常 B ∝ r −3 のダイ
ポール磁場の減衰勾配を B ∝ r −1.51 程度になることが確認できた。またコイル中心から 15km の位置にある 41nT の磁場を、コイルから 100km 以上の位置まで展開することができた。
・ 統合解析においては、噴射プラズマのラーマ半径と磁気圏半径の比 rLi / L 、
噴射プラズマ β i などの無次元パラメータを変化させて結果を比較した。これ
より、磁気レイノルズ数が無限大の条件のもとでは、現実的なパラメータである β i <1, rLi / L >1 においても磁場は展開され、推力の増分が得られることが分かった。

42

以上のような結果が得られたが、磁気インフレーション解析や統合解析においては本ハイブリッドコードの制限から、依然として次のような問題を抱えている。
本来の MPS のアイデアでは初期磁場は小さく、コイル半径は大きくても 10m程度と考えられている。しかし正確なダイポール形状を描くためには、コイル
半径に 1 メッシュ以上使用することが必要であり、数十 km 程度まで 41nT の磁場を広げるような計算を行うためには、1 メッシュ 1m とすると１軸に対し数万メッシュ必要であり、我々の計算資源ではおよそ 107day 以上の計算時間が求められてしまい非現実的である。これに対する解決方法としては、不均一のメッシュや球座標を用いたコードに改良することで解決できると考えられる。

また、本ハイブリッドコードは、磁場を求める際に使用するアンペアの式に
ダーウィン近似を用いていることにより、 (2-9)で求められる電場の値が、低密度、高磁場、メッシュ幅が小さいなどの条件で発散してしまう。これにより生じる問題を以下に挙げる。
① 太陽風のような低密度領域において体系が 10km 以下小スケールの計算を行えない。
② 本来はコイル付近の磁場は 0.01T ほどと考えられているが [5]、本解析コード
では、太陽風などの低密度 (5×106m-3)粒子を扱う際に、磁場の上限値がおよそ 300nT と制限されてしまう。
③ コイル付近の磁場強度が小さい場合、β i 値を基準にした際の噴射プラズマの数密度は太陽風と同程度の小さな値となってしまう。最近の MHD 解析によれば、 β i <10-6 で噴射しなければイオンエンジンと比較して推力電力比を高くすることはできないとされている [13]。本解析コードで、そのような β i 値を用いて計算しようとすれば、噴射プラズマの数密度は太陽風以下となってしまう。
(2-40) 式で求められる磁場の分散周波数 ω D はこの発散現象の指標となっており、この値の逆数以下の時間刻み ∆t で計算することが望ましいとされている[22]。MPS 解析で必要とされるパラメータ ∆x =1m, N =5×106m-3, B =0.01T を用い

ると ∆t <10-17s が要求されてしまい、膨大な計算時間が必要である。
43

前項の電場発散現象に対して現段階で考えられる対策案の一つは、イオン、電子ともに粒子として扱う Full-PIC コードを用いることである。同じ粒子コードでありながら、こちらにはダーウィン近似を用いたことによる問題が生じな
いため、 MPS 解析に必要な低密度、高磁場、小スケール（磁気セイルと比較して）の体系を容易に実現できるメリットがある。このコードには荷電分離の観点からデバイ長の制限があるが、低 β （低密度）の噴射プラズマを用いた解析ならばそれも緩和できると考えられる。
ハイブリッドコードを用いた中での対策法としては、並列計算を行うことで処理速度を高め計算領域を拡大する方法や、インフレーションを MHD で解析し、その結果をハイブリッドコードの初期条件に引き渡し、太陽風との相互作用を解析するような方法が考えられる。
これらをふまえ、より現実のモデルに近い体系での評価を実施する必要があ
る。

44

参考文献
[1] Internet Homepage，http://www.jaxa.jp．
[2] 栗本恭一，荒川義博，「電気推進ロケット入門」，東京大学出版会，2003.
[3] アストロアーツ「Space Guide 宇宙年鑑 2005」[4] R. M. Zubrin, D. G. Andrews, “Magnetic Sails and Interplanetary
Travel,”Journal of Spacecraft and Rockets, Vol.28, No.2, p.3118, 1991.

[5] R. M. Winglee, J.Slough, T. Ziemba, A. Goodson, “Mini-Magnetospheric
Plasma Propulsion: Tapping the energy of the solar wind for spacecraft
propulsion”, Journal of Geophysical Research, Vol.105, No.21, pp.67-77,
2000.
[6] P. K. Toivanen, P. Janhunen, H. E. J. Koskinen, “Magnetospheric
Propulsion”, ESTEC contract, No. 16361, 2004.
[7] 山川宏，小川博之，藤田和央，野中聡，“磁気プラズマセイルによる外惑星探査”,日本航空宇宙学会論文集, Vol.52, No.603, 2004.
[8] 恩藤忠典，丸橋克英「宇宙環境化学」，オーム社，2000.
[9] Internet Homepage, http://www2.nict.go.jp．
[10] 高村秀一「プラズマ理工学入門」，森北出版株式会社，1997.
[11] 西田浩之，稲谷芳文，小川博之，船木一幸，“磁気セイルの推進特性解析”,
ISAS Research Note, ISAS RN 791, 2005.
[12] 朝日龍介，船木一幸，山川宏，藤田和央，“磁気プラズマセイルの電磁流体解析および性能評価”，ISAS Research Note, ISAS RN 789, 2005.
[13] 永田靖典，大津広敬，小川博之，船木一幸，山川宏，藤田和央，“磁気プラズマセイルの 3 次元 MHD 解析”，第 49 回宇宙科学技術連合講演会予稿集，
2005.
[14] K.

Fujita,I.

Funaki,H.

Ogawa,H.

Yamakawa,“NumericalInvestigations of Acceleration Process in Magneto Plasma Sail”, ISTS
2004-o-3-07v, June 2004.
[15] K. Fujita, “Particle Simulation of Moderately-Sized Magnetic Sails”
Journal of Space Technology and Science, Vol.20, No.2, pp.26-31, 2005.
[16] N. Omidi, H. Karimabadi, “Kinetic Simulation / Modeling of Plasma
Sail”, AIAA paper 2003-5256, pp.20-23, 2003.
[17] G. Khazanov, P. Delamere, K. Kabin, T. J. Linde “Fundamentals of the
Plasma Sail Concept: Magnetohydrodynamic and Kinetic Studies”,
45

Journal of Propulsion and Power, Vol.21, No.5, 2005.
[18] D. S. Horowiz, J. Comput. Phys., Vol.47, p.452, 1982.
[19] E. J. Horowiz, D. E. Schumaker, D. V. Anderson, J. Comp. Phys.,Vol.

84, p.279, 1989.
[20] A. Mankofsky, R. N. Sudan, J. Denavit, J. Comp. Phys., Vol.70,
p.89,1987.
[21] C. K. Bridsall, A.B. Langdon, “Plasma Physics via Computer Simulation”,
IOP Publishing Ltd, 1997.
[22] A. P. Matthews, “Current Advance Method and Cyclic Leapfrog for 2D
Multispecies Hybrid Plasma Simulations”, Journal of Computational
Physics 112, pp.102-116, 1994.
[23] 大林辰蔵，「宇宙空間物理学」，裳華房，2003.

46