卒

題

目

業

論

文

整形ターゲットを用いたレーザー核融合ロケットの推進効率の向上

氏

名

松田伸夫

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

指 導 教 員

中島秀紀

教授

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

平成 19 年 2 月

目次
第 1 章 序論................................................................................................................................1
1.1 背景 ...................................................................................................................................1
1.2 目的 ...................................................................................................................................2
第 2 章 レーザー核融合推進システム.......................................................................................4
2.1 レーザー核融合ロケット ....................................................................................................4
2.2 推進原理 ...........................................................................................................................5
第 3 章 シミュレーションコード....................................................................................................6
3.1 SPH 法について................................................................................................................6
3.1.1 衝突解析における数値解析法..................................................................................6
3.1.2 SPH 法の概念 ............................................................................................................7
3.1.3 SPH 法の解析手法 ....................................................................................................8
3.1.4 フローチャート ..........................................................................................................17
3.2 3D hybrid code について ................................................................................................18
3.2.1 はじめに....................................................................................................................18
3.2.2 基礎方程式 ..............................................................................................................18
3.2.3 計算方法 ..................................................................................................................20
3.2.4 プラズマ粒子の運動 ................................................................................................20
3.2.5 イオン数密度、電流密度 .........................................................................................21
3.2.6 電磁場の計算方法 ..................................................................................................22
3.2.7 Digital Filter .............................................................................................................24
3.2.8 Hybrid code に課せられる制限 ...............................................................................25
第 4 章 計算結果......................................................................................................................26
4.1 評価方法 .........................................................................................................................26
4.2 球状ターゲット .................................................................................................................26
4.2.1 計算体系 ..................................................................................................................26
4.2.2 計算結果 ..................................................................................................................28
4.3 整形ターゲット .................................................................................................................30
4.3.1 計算体系 ..................................................................................................................30
4.3.2 計算結果 ..................................................................................................................32
第 5 章 結論..............................................................................................................................35
参考文献 .......................................................................................................................................37
謝辞 ...............................................................................................................................................38

第1章
1.1

序論

背景

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

Fig. 1-1

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

レーザー核融合ロケットについては第 2 章で述べる。
磁気ノズルは過去に設計、シミュレーションが行われてきた。Hyde は、超伝導(SCM)コイルを2 つ用いたモデルを提案している(1)。また、2 次元 MHD シミュレーションを用いてプラズマ挙動と推進効率を計算し、1 つのコイルで推進効率 65％という結果を得ている(2)。長峯と中島(3)は 3次元(3D) hybrid code を用いてプラズマ挙動と推進効率を計算し、推進効率 65％という結果を得ている。Vchivkov(4)は同様のコードで推進効率の最適化を行い、推進効率 70％という結果を得ている。また、実験との比較も行い、推進効率の一致を得ている (5)。また、川渕により 3D
hybrid code を用いて核融合プラズマの発生位置をずらしたり、２つのコイルを用いることでロケットの姿勢・推進方向制御についての研究が行われた(6)。
1

1.2

目的

本研究では、さらに推進性能を向上すべく新たな方式を検討する。それは、ターゲットを整形しプラズマの噴出方向を変えることにより、推進性能の向上を目指す。(Fig1-2 参照)
核融合燃料ペレットの周囲に推進剤(Moderator)を配置し、レーザー核融合で発生したプラ
ズマと推進剤の衝突により発生するプラズマの運動量を磁場発生用コイル方向(磁場の強い方向)に集中させることで推進効率の向上を目指す。

膨張過程

Fig. 1-2

以前、川渕がプラズマ粒子の噴出が球状の場合と円錐の場合を 3D hybrid code で計算させ、それぞれの推進効率を求めた(6)。Fig.1-3 は磁気ノズル中での初期配置と膨張後の図である。

Fig.1-3 上段：磁気ノズル中での初期配置

2

下段：膨張後

この計算により、それぞれの推進効率は球状の場合が 69%、円錐の場合が 73%となっており、
運動量を集中させることで推進効率の上昇をシミュレーションできた。
このシミュレーションはターゲットが膨張して円錐の形になったものとして計算させている。
そこで本研究では、Smoothed Particle Hydrodynamics(SPH)法を用いて核融合プラズマと推進剤の衝突を計算させ、衝突後の粒子の位置、速度を求める。そのデータを 3D hybrid codeの粒子の位置、速度データとして入力し、磁気ノズル中でのプラズマの挙動を計算し、推進効率を求める。

この論文の構成は以下の通りである。
第 1 章 本研究の背景と目的を述べた。
第 2 章 レーザー核融合推進システムについて述べる。
第 3 章 本研究に用いた計算コードについて述べる。
第 4 章 解析結果について述べる。
第 5 章 結論を示す。

3

第2章
2.1

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

レーザー核融合ロケット

現在、NASA などの宇宙機関において火星への有人探査の計画が進められている。火星は、人類が地球圏外へ進出する際、比較的早い段階に移住できると考えられる唯一の天体である。
火星の大気には二酸化炭素が豊富に存在し、大量の水が凍った状態で眠っているからである。
また、火星は人類にとって次の居住地となるばかりではなく、火星での経験が、他の惑星の探査・開発にとって重要なものとなり、前哨基地としての役割も担うことになるだろう。しかし、近い将来、火星に恒久的な基地や居住施設を建設することを考えると、従来の宇宙輸送手段である化学推進のロケットでは建設に必要な膨大な量の資材を運搬することは容易ではない。
そこで、その際に使用される宇宙輸送手段として考えられているのが、レーザー核融合ロケットである。1972 年、Hyde によって始めて超伝導(SCM)コイルを 2 つ用いたレーザー核融合ロ
ケットのモデルが提案された(1)。ロケットの推進システムにエネルギー密度が大きい核融合を用いると、化学推進に比べて高い比推力と高い比出力を同時に得ることができる。推進材の排出速度は 107m/sec 程度になると見られ、比推力（排出速度を重力加速度で割った値）は 106sec となる。この値は、化学推進ロケットの 2000 倍程度であり、ペイロード（積載貨物）の割合が多く取れることになる。さらに、ロケットの外部空間をエンジンとして利用するので、その分宇宙船に占めるエンジンの質量が小さいので高比出力が達成されると考えられる。これは、輸送時間の短縮にもつながり、宇宙線による被爆など人体への悪影響も抑えられる。
しかし、レーザー核融合ロケットの実現には多くの課題がある。レーザー核融合を実現するためのレーザードライバ、大型超伝導コイル、大量の核融合燃料生成などの開発が急務であるが、これらは現在活発に研究が行われている地上発電用核融合技術がそのまま流用できる。
さらに、2003 年に Orth が提案したレーザー核融合ロケットの概念設計(7)では、宇宙船全体の質量が 6500ton もあり、2005 年に本研究グループが行った新しい概念設計でも 250ton の巨大な宇宙船になるという計算結果が導き出された。この巨大な構造物を宇宙で組み立て、さらに高速での運転を行うための技術も必要である。

4

2.2

推進原理

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

Fig. 2-1

磁気ノズルでの推進原理

5

第3章
3.1

シミュレーションコード

SPH 法について

ここでは超高速衝突における数値解析法、Smoothed Particle Hydrodynamics(以後 SPH)法
の概念、SPH 法の解析手法等について述べる。
また本研究で使用している SPH 法は G. R. Liu と M. B. Liu のものを参考にしている(8)。

3.1.1

衝突解析における数値解析法

衝突解析を行う際には衝突物体を流体とみなして解析を行うが、その流体の運動を表す方
法として一般的なのは Euler と Lagrange の方法である。Euler 法は固定座標系であるため基本
的にどのような問題でも解けるが、物体がない場所も計算しなければならないため無駄な計算時間がかかるという問題がある。また、Lagrange 法は移動座標系であるため物体のあるところだ
け計算を行えばよい。そのため Euler 法と比べて計算時間が早く、非定常流などの対称性のよい流れを扱う際などには安定性はよいが、大変形を伴う解析時にはメッシュの潰れや歪みが生じるため、その都度リメッシュをかけなければならないという問題がある。これらの解析手法はメッシュによる制限があり、解析にかかる時間や解析結果の精度等がメッシュに依存したものであるということがいえる。
これらの双方の利点を取り入れたもっとも効率のよい数値解析法を見い出すため、さまざまな方法が試されてきた。例えば、異なる２つの物体を取り扱う場合に片方を Euler 座標にとり、もう一方を Lagrange 座標にとるというハイブリッド法などがある。これらの手法は複雑な計算に対して効果を発揮するが、同時にさらなる不正確さも持ち合わせている。
そこで近年盛んに開発が行われているのがメッシュレス解析法である。メッシュレス解析とは名前のとおりメッシュを作らずに解析を行う手法で Element Free、Mesh Free、Meshless、
Gridless 、 Pointless な ど が あ る 。 そ の 中 で も 代 表 的 な も の に Belytschko ら の 開 発 し た
EGF(Element Free Galerkin)法がある。この手法は移動最小自乗近似(MLS : Moving Least
Square Interpolation)が基礎となっており、物体を表すのに節点データのみが必要であり要素の連続性というものは必要としない。これによりさきほど述べたリメッシュをしなくてもよいため、大変形を伴う変形問題を解くような問題の際には非常に有用である。解き方としてはまず接点配置を決定し、その後、移動最小自乗近似により各接点ごとに形状関数を決定し、次に、
Hu-Washizu の原理を用いて剛性マトリクス、荷重ベクトルを計算し、最終的に線形方程式を解き、求めたい点での変位、応力の値を計算するというものである。本研究で用いている SPH 法
も EGF 法の一種である。EGF 法では各節点形状関数を決定する際、移動最小自乗近似が用いられているのに対して、SPH 法では Gauss 分布関数等であらかじめ決定しておいた平滑化関数を用いる。

6

3.1.2 SPH 法の概念
流体力学現象を計算させるための解析手法として SPH 法は 1977 年に Lucy により考案された。SPH 法は、本来、非軸対称な天文物理学の問題を解くために開発された。改良が重ねられ、降着円盤、星形成、超新星爆発、星同士の衝突、銀河の形成、銀河の相互作用を含むさ
まざまな天文物理学の研究において幅広く使われるようになった。その中で中心となっていたのが 1970 年代から 1980 年代前半までの Gingold、Monaghan ら、また 1980 年代の中頃から
1990 年代の Benz、Monaghan らであった。近年になり、SPH 法は天体物理学以外の分野の問題に使用されるようになってきた。これらは主に流体を扱っていたのに対し、1990 年の前半にLibersky らが流体解析だけでなく固体の問題も解析できるように改良を加えたことによりさまざ
まな分野において SPH 法が用いられるようになってきた。最近になって大変形を伴う解析に使用できる利点から超高速衝突解析にも用いられるようになってきた。
SPH 法の最大の利点は偏微分を計算する際に格子を用いないため、格子の絡まりについて考えなくてもよいという点があげられる。この利点により Lagrange 座標を用いて大変形を伴う衝突や爆発を引き起こすような現象を解く際に生じるメッシュの潰れや歪みとさまざまな問題を解消することができる。Lagrange 座標や Euler 座標などの格子をもとに解析を行う手法は偏微分を計算するために節点と節点の連続性を仮定する。しかし、SPH 法ではどの点が近傍粒子かという仮定を持たないランダムに配置された内挿点に基づいた kernel 近似を用いることにより偏微分を計算する。Fig.3-1 に SPH 法の考えをもとにして空間に配置した粒子のイメージを表す。

Fig.3-1

SPH 法を基に空間に配置した粒子のイメージ

7

粒子” i “は 2h の半径内にあるその他の粒子” j “からの影響を受けるが、2h の半径外の粒子からの影響は受けないとい。この h を平滑化距離(Smoothing Length)と呼ぶ。それぞれの粒子は平滑化関数(Smoothing Function)である W (xi − x′j , h ) により重みを掛けられている。つまりSPH 法は連続体と近似した計算対象を素片に分割してその素片の運動を追跡する手法で、各素片は互いに重なり合ってその集合体として、例えば粒子群を記述するといった考え方に基づいている。素片の持つ物理量は素片の中心における値として表され、その点と共に移動する。
そして粒子の影響範囲内で分布関数によって分布され、その重なり合いの総和として着目仮想粒子の物理量が求められる。粒子の存在点によって連続体を表現し、その密度に比例して連続体を緻密に表すように粒子の影響範囲が決められるので、密度の高いところは詳しく、薄いところは大きい粒子によって粗に表現されるので、計算も効率的である。

3.1.3 SPH 法の解析手法
a) 基本式
SPH 法では、位置 x の物理量 f (x) を表す時に、物理量 f (x) そのものではなく平滑化関数
W ( x − x′, h) によって平滑化された物理量を表す式(3.1-1)を用いる。
f ( x) = ∫ f ( x′)W ( x − x′, h)dx′(3.1-1)同様に１階微分 ∇f (x) は W ( x − x′, h) の１階微分 ∇W ( x − x′, h) を用いて式(3.1-2)のように表せる。

∇f ( x) = ∫ f ( x′)∇W ( x − x′, h)dx′(3.1-2)

ここで、 W ( x − x′, h) は規格化のため式(3.1-3)の関係を満たさなければならない。

∫ W ( x − x′, h)dx′ = 1(3.1-3)

また、 W ( x − x′, h) の広がりは h の大きさによって決められ、もし h➝0 ならば式(3.1-4)のようになる。
f ( x) ⎯h⎯⎯→ f ( x)→0(3.1-4)

一般的に平滑化関数は x − x′ = 0 の部分でピークをとる関数であるので、h が 0 に近づくと Diracのδ関数と同等となるため f (x) は f (x) と等しいとみなすことができる。これにより、 f (x) は n個の離散点によって表せると考えることができるため式(3.1-5)のような数密度分布により空間的に分配した形で表せる。N は粒子数である。
n( xi ) = ∑ δ (xi − x j )
N

j =1

8

(3.1-5)

これを粒子 j の有限の体積 ∆V j と結びつけて考えると式(3.1-6)が得られる。

∆V j =

mj(3.1-6)ρj

ここで ρ j は粒子 j の密度である。式(3.1-1)を式(3.1-6)を用いて書き換えると式(3-7)が得られる。
f ( x) = ∫ f ( x ′)W ( x − x ′, h)dx ′N≅ ∑ f ( x j )W ( x − x j , h)∆V j
j =1

N

mj

j =1

ρj

f ( xi ) = ∑f ( x j )W ( x − x j , h)

(3.1-7)

このような近似を用いることにより Lagrange 方程式において加速度や歪み速度などを見積もることができるということが SPH 法の基礎となっている。
例として流体の密度を表現することを試みる。実際の粒子の集合を考えると、式(3.1-7)より密度は式(3.1-8)となる。
N

ρ ( xi ) = ∑ m jW ( x − x j , h)

(3.1-8)

j =1

ここで m j は流体素片の質量であり、 x j は j 番目の粒子の重心を表す。つまり、N 個の分布関数の重ね合わせとして、言い換えると、ある広がりを持った粒子の重ね合わせとして流体の密度を決定する。(概念図 Fig.3-2)

Fig.3-2 流体の密度の概念図

9b) Navier-Stokes 方程式
動的流体の流れについての支配方程式は一組の偏微分方程式で表すことができる。その一組の偏微分方程式は、よく知られている Navier-Stokes 方程式であり、その式は質量保存、運動量保存、エネルギー保存である。
1) 連続の式

2)

3)

Dρ
∂v β
= −ρ βDt∂x(3.1-9)Dv α
1 ∂σ αβ=ρ ∂x βDt(3.1-10)De σ αβ ∂v α=ρ ∂x βDt(3.1-11)

運動方程式

エネルギー方程式

上の式で σ は全応力テンソルである。これは二つの部分からなり、一つは等方性の圧力 p で、もう一つは粘性応力 τ である。
σ αβ = − pδ αβ + τ αβ
(3.1-12)ニュートン流体において粘性せん断応力は動的粘性係数μを介してεで表されるせん断歪率に比例しなければならない。
τ αβ = µε αβ(3.1-13)
ここで∂v β ∂v α 2+ β − (∇ ⋅ v)δ αβ(3.1-14)α3∂x∂x
等方性の圧力と粘性応力が分割し考える場合、エネルギー方程式は次のように書き換えられる。

ε αβ =µ αβ αβDep ∂v β
+
ε ε
=−
β
ρ ∂x
2ρ
Dt

10

(3.1-15)c) 構成式i) 連続の式
前述のとおり、SPH 法において密度の式は式(3.1-8)だが、簡略化すると式(3.1-16)となる。
N

ρ i = ∑ m jWij

(3.1-16)

j =1

ここで、 Wij = W ( xi − x j , h) = W ( xi − x j , h) = W ( Rij , h) である。
Rij =

rij

=

xi − x jは粒子 i と j の相対距離であり、 rij は二粒子間の距離である。
hh
式(3.1-16)は粒子の密度が重みをかけられた粒子の密度の平均で近似できることを端的に表している。
粒子を密度について近似するもう一つの方法は密度の連続性を考慮することである。これはSPH 法の概念を用いた流体力学の連続の式に従って近似するというものである。
流体力学の連続の式は式(3.1-17)で表される。これは式(3.1-9)と同じものである。

Dρ∂v β= −ρ βDt∂x
これに SPH 法の概念を用いると式(3.1-18)となる。
N m
∂Wij
Dρ i
j
v βj ⋅ β= −ρi ∑
Dt∂xi
j =1 ρ j(3.1-17)(3.1-18)

SPH 法のカーネルと 1 の勾配での粒子の近似を考慮すると式(3.1-19)を得る。
∇1 = ∫ 1 ⋅ ∇W ( x − x ′, h)dx ′m j ∂Wij⋅ β =0∂x j
j =1 ρ j(3.1-19)⎛ N m j ∂Wij ⎞⎟= ρ i viβ ⎜ ∑⋅⎜ j =1 ρ ∂x β ⎟
∂xi
j
i ⎠
⎝

(3.1-20)N=∑

これは次のように書き換えられる。
N

mj

j =1

ρj

ρi ∑

viβ ⋅

∂Wij

式(3.1-20)を式(3.1-18)に代入すると式(3.1-21)を得る。
N m
∂Wij
Dρ i
j
= ρi ∑
vijβ ⋅ βDt∂xi
j =1 ρ j

ここで、 vij = vi − v j である。

11

(3.1-21)

より一般的な粒子の連続性の形は勾配演算子の中に密度を置いた次の恒等式を用いる。
−ρ

(

)

⎛ ∂ ρv β
∂v β
∂ρ ⎞
⎜⎜
=
−
− v β ⋅ β ⎟⎟
β
β
∂x
∂x ⎠
⎝ ∂x

(3.1-22)

同様に、SPH 法による近似をあらゆる勾配に適用し、二番目の勾配における速度が、その勾配が評価される粒子で計算されるならば、最もよく使用される式(3.1-23)のような密度の連続の式
を得る。
N
∂Wij
Dρ i
= ∑ m j vij ⋅(3.1-23)
Dt
∂xi
j =1
この式は、粒子の密度の時間変化の割合が密接に領域におけるある一つの粒子と他のすべての粒子の間の相対的な速度に関連するのを明確に示している。平滑化関数の勾配はこれらの相対的な速度の寄与を決定する。

ii) 運動方程式
運動量の式[式(3.1-10)]の勾配に SPH 法の概念を適用すると式(3.1-24)を得る。

∂Wij
σ αβ
Dviα
1 N
j
=
m
∑
j
Dt
ρ i j =1
ρ j ∂xiβ(3.1-24)

上式に次の恒等式(3.1-25)を加えると、式(3.1-26)を得る。

σ iαβ ∂Wij σ iαβ ⎛⎜ N m j ∂Wij ⎞⎟
=
=0
mj
∑
∑
ρ i ρ j ∂xiβ
ρ i ⎜⎝ j =1 ρ j ∂xiβ ⎟⎠
j =1
N

N
∂Wij
σ iαβ + σ αβ
Dviα
j
= ∑mj
Dt
ρi ρ j∂xiβ
j =1(3.1-25)(3.1-26)次の恒等式(3.1-27)を考慮し SPH 法による近似を加えると式(3.1-26)は式(3.1-28)のようになる。
1 ∂σ αβ∂ ⎛ σ αβ ⎞ σ αβ ∂ρ⎟+⎜=ρ ∂x β∂x β ⎜⎝ ρ ⎟⎠ ρ 2 ∂x β

(3.1-27)αβN m∂Wij σ iαβ N m j
∂Wij
Dviα
j σi
=∑
+
ρ
∑
j
β
Dt
ρ 2j j =1 ρ j∂xiβj =1 ρ j ρ j ∂xi(3.1-28)

簡単な再配列の後、式(3.1-28)は次のように書き換えられる。
N⎛ σ iαβ σ αβ⎞ ∂WijDviαj= ∑mj⎜ 2 + 2 ⎟ β⎜ ρDtρ j ⎟⎠ ∂xi
j =1
⎝ i

12

(3.1-29)

式(3.1-12)を用いると、式(3.1-26)と式(3.1-29)はそれぞれ次のように書き換えられる。
NN
pi + p j ∂Wij
∂Wij
µ i ε iαβ + µ j ε αβDviαj= −∑ m j
+
m
∑
j
α
Dt
ρ i ρ j ∂xi
ρi ρ j
∂xiβ
j =1
j =1

(3.1-30)

NN⎛ pi⎛ µ i ε iαβ µ j ε αβ⎞ ∂Wijp j ⎞ ∂Wij
Dviα
j
⎜
⎟
⎜
⎟
= −∑ m j 2 + 2
+
m
+
∑
j
α
2
2
⎜
⎟
⎜
⎟ ∂x α
Dt
ρ
ρ
∂
x
ρ
ρ
j =1
j =1
j ⎠
i
j
⎝ i
⎝ i
⎠ i

(3.1-31)式(3.1-30)と式(3.1-31)の第一項は圧力についての SPH 法による近似であり、第二項は粘性力についての近似である。また式(3.1-14)において、粒子 i について ε αβ の SPH 法による近似を考慮すると次のように書き換えられる。
N

mj

j =1

ρj

ε iαβ = ∑

v βj

∂Wij

N

mj

∂xi

j =1

ρj

+∑
α

⎛ 2 N mj
⎞ αβ
⎜ ∑
⎟δ
−
⋅
∇
v
W
j
i
ij
⎟
∂xiβ ⎜⎝ 3 j =1 ρ j
⎠

∂Wij

v αj

(3.1-32)式(3.1-19)を考慮すると、次の恒等式を得る。

⎛ N m j ∂Wij ⎞

∂Wij

mj

N

∑ ρ v β ∂x α = v β ⎜⎜ ∑ ρ ∂x α ⎟⎟ = 0(3.1-33)⎛ N m j ∂Wij ⎞
α⎜
⎟=0
=
v
i ⎜∑
β ⎟
ρ
∂xiβ
∂
x
j
=
1
j
i
⎝
⎠

(3.1-34)

j =1

j

N

mj

j =1

j

∑ρ
N

i

i

i

viα

⎝ j =1

j

i

⎠

∂Wij

⎛ N mjmj⎞∑ ρ v ⋅ ∇ W = v ⋅ ⎜⎜ ∑ ρ ∇ W ⎟⎟ = 0
j =1

i

i

ij

i

⎝ j =1

j

i

ij

j

⎠

(3.1-35)

粘性のひずみ率に速度差を含めるために式(3.1-32)からこれらの恒等式を引くと、粒子 i における、 ε αβ の最終的な SPH 法による近似は式(3.1-36)である。
N

mj

j =1

ρj

ε iαβ = ∑

v βji

∂Wij

N

mj

∂xi

j =1

ρj

+∑
α

v αji

⎛ 2 N mj⎞
⎜ ∑−v ji ⋅ ∇ iWij ⎟δ αββ⎜ 3 j =1 ρ
⎟
∂xi
j
⎝
⎠

∂Wij

(3.1-36)

粒子 j についても同様にして求めることができる。粒子 i と j についての ε αβ を計算した後に、式(3.1-30)か式(3.1-31)により加速度を計算することがでる。

iii) エネルギー方程式
式(3.1-15)での内部のエネルギーe の発展に関しては、式(3.1-36)を使用することでひずみ率にかかわる部分を近似することができる。圧力仕事にかかわる第一項には対応する SPH 法の公式を導くいくつかの異なった方法がある。
−

p ∂v βp ⎛p Dρ
∂v β ⎞
⎟= 2
⎜
ρ
=
−
β
2 ⎜
β ⎟
ρ ∂x
ρ ⎝
∂x ⎠ ρ Dt

13(3.1-37)

式 (3-37) の関係を考慮すると、関係を持つ粒子でpρ2の値を求めながら、圧力仕事は式(3.1-18) や式 (3.1-21) 、式 (3.1-23) の密度の連続性を用いることにより近似できる。例えば、式(3.1-23)の密度の連続性を用いると、圧力仕事は次式のように表される。
∂Wij
pi N
p ∂viββmv−=⋅(3.1-38)∑jij
ρ ∂xiβ ρ i2 j =1∂xiβ
別の方法は次の恒等式を用いる。
p ∂v β∂ ⎛ pv β ⎞ β ∂ ⎛ p ⎞
⎟+v
⎜
⎜ ⎟
=
−
ρ ∂x β∂x β ⎜⎝ ρ ⎟⎠∂x β ⎜⎝ ρ ⎟⎠−(3.1-39)

すると次の式を得る。
−

N
∂Wij
pj
p ∂viβ=m j 2 vijβ ⋅ β
∑
β
ρ ∂xi
ρj
∂xi
j =1

(3.1-40)式(3.1-38)と式(3.1-40)を加え一般的な形に変形すると、−関係を持つ粒子で求めた

⎛ p
pj ⎞
∂Wij
p ∂viβ 1 N= ∑ m j ⎜ 2i + 2 ⎟vijβ ⋅ ββρ ∂xi
2 j =1 ⎜⎝ ρ i
ρ j ⎟⎠∂xip(3.1-41)の値と式(3.1-21) を用いることで、式 (3.1-37)の時間の増減率にρ2ついて近似すると次式を得る。

pi N m j β ∂Wijp ∂viβ=∑ vij ⋅ ∂x β
ρ ∂xiβ ρ i j =1 ρ ji(3.1-42)p ∂v β1⎡ ∂∂p ⎤−= − ⎢ β ( ρv β ) − v β β ⎥βρ ∂xρ ⎣ ∂x∂x ⎦(3.1-43)−

また式(3.1-43)の変形を用いることで圧力仕事についての SPH 法の近似の別の形が求まる。

−

p j β ∂Wijp ∂viβ1 N
mjvij ⋅ β=∑βρ ∂xi
ρ i j =1
ρj∂xi

(3.1-44)式(3.1-42)と式(3.1-44) を加え一般的な形に変形すると式(3.1-45)を得る。

−

pi + p j β ∂Wijp ∂viβ 1 N
vij ⋅ β
= ∑mj
β
ρ ∂xi
2 j =1ρi ρ j∂xi(3.1-45)

圧力仕事を近似するには多くの形式があるので、内部エネルギーの発生のための SPH 法の公式は自然と多くの代替手段を持つことができる。

14

もっともよく使われるのは次の二式である。

⎛ p
pj ⎞
∂Wij
µ
∂ei 1 N= ∑ m j ⎜ 2i + 2 ⎟vijβ ⋅ β + i ε iαβ ε iαβ∂t2 j =1 ⎜⎝ ρ i
2ρ i
ρ j ⎟⎠∂xi(3.1-46)

pi + p j β ∂Wij∂ei 1 Nµvij ⋅ β + i ε iαβ ε iαβ= ∑mj
ρi ρ j∂t
2 j =12ρ i∂xi(3.1-47)

また、密度、運動量、およびエネルギーを導くためのオイラー方程式に対応する SPH 法のそれ
ぞれの方程式は以下のとおりになる。
連続の式
N
∂Wij
Dρ i
= ∑ m j vijβ ⋅ β(3.1-23)
Dt
∂xi
j =1
運動方程式
N
pi + p j ∂WijDviα= −∑ m jDtρ i ρ j ∂xiα
j =1

(3.1-48)N⎛ p
p j ⎞ ∂WijDviα= −∑ m j ⎜ 2i + 2 ⎟ α⎜ρDtρ j ⎟⎠ ∂xi
j =1
⎝ i

(3.1-49)

⎛ p
pj ⎞
∂Wij
∂ei 1 N= ∑ m j ⎜ 2i + 2 ⎟vijβ ⋅ β∂t
2 j =1 ⎜⎝ ρ i
ρ j ⎟⎠∂xi(3.1-50)pi + p j β ∂Wij∂ei 1 Nvij ⋅ β= ∑mj∂t
2 j =1ρi ρ j∂xi(3.1-51)

エネルギー方程式d) 人工粘性
SPH では人工粘性を考慮している。レーザー核融合で生じたプラズマと推進剤との衝突を
SPH 法で記述しようとすると Fig.3-3 のように相対速度が音速を超える粒子群のぶつかり合いとして表される。

Fig.3-3

SPH 法による衝撃波面近傍の模式図

15しかし、音速を超えている場合には粒子が粒子をすり抜けるという非物理的現象が起こる。
それを回避するために衝撃波面にのみ働く人工粘性を考慮している。
人工粘性は次式で表される。

⎧ − α Π cij φij + β Πφij2⎪
ρ ij
Π ij = ⎨⎪0⎩vij ⋅ xij < 0(3.1-52)vij ⋅ xij ≥ 0

ここで、

φij =

hij vij ⋅ xij
xij

2

+ϕ2

(3.1-53)1(c i + c j )2(3.1-54)

ρ ij = ( ρ i + ρ j )1
2(3.1-55)

1(hi + h j )2(3.1-56)

cij =

hij =

vij = vi − v j , xij = xi − x j(3.1-57)

上式の α Π と β Π は定数である。因数 ϕ = 0.1hij は二つの粒子が近づいているときに数値発散を防ぐために挿入されている。c と v はそれぞれ音速と粒子の速度ベクトルを表している。
この人工粘性を SPH 法の公式に導入すると SPH 法で計算に使用される式は以下の一組の式のようになる。
N⎧ Dρ i=⎪ Dt ∑ m j (vi − v j ) ⋅ ∇Wij
j =1⎪
⎪ Dv αN⎛ p⎞p⎪ i = −∑ m j ⎜ 2i + 2j + Π ij ⎟∇Wij
⎜ρ
⎟
ρj
⎪ Dt
j =1
⎝ i
⎠
⎨
⎛ p
⎞
⎪ Dviα 1 Npj= ∑ m j ⎜ 2i + 2 + Π ij ⎟(vi − v j ) ⋅ ∇Wij⎪⎟
2 j =1 ⎜⎝ ρ i
ρj
⎪ Dt
⎠
⎪ Dx
⎪ i = vi
⎩ Dt

16

(3.1-58)

3.1.4 フローチャート
SPH 法を用いたシミュレーションコードのフローチャートを示す。
粒子の初期化(INPUT)

近傍粒子の検索

平滑化計算

相互作用計算人工粘性項計算

エネルギー、運動量の変化を計算

粒子の位置、速度、エネルギー、密度を計算

データ書き出し(OUTPUT)

3D hybrid code へ
Fig.3-4 フローチャート

17

3.2

3D hybrid code について

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

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

dx i
= vi
dt

(3.2-2)

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

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

ここで v e は電子の速度、 me は電子の質量、 ne は電子数密度、 Pe は電子圧力である。電子は慣性無しの流体としてみなすので、慣性による効果は無視し、(3.2-3)式の左辺は 0 にできる

( me → 0 )。さらに電子圧力 Pe に関しては電子の状態方程式より以下のようになる。

18

Pe = neTe

(3.2-4)

ここで Te は電子温度である。(3.2-4)式を(3.2-3)式に代入して Pe を消去する。

0 = −ene (E + v e × B ) − ∇(neTe )(3.2-5)

荷電準中性条件および、ダーウィン近似を用いたアンペアの法則により以下の式が得られる。

∇ × B p = µ 0 (J e + J i )(3.2-6)

ここで µ 0 は真空中の誘電率であり、 J i 、 J e はそれぞれイオン電流密度、電子電流密度である。また、 B p はプラズマにより生成される磁場で以下のようになる。

B p = B − B ext − B pick

(3.2-7)

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

∇ × B pick = µ 0 J pick

(3.2-8)

但し、 J pick はピックアップコイルに流れる電流密度である。

J i はイオンの数密度と電流密度から計算することが出来る。 J e は電子電流密度の定義によ
り以下のようになる。

J e = −ene v e(3.2-9)(3.2-6)、(3.2-9)式を用いて J e を消去し v e について求めると以下のようになる。
ve = −⎞1 ⎛ 1⎜⎜ ∇ × B p − J i ⎟⎟ene ⎝ µ 0⎠(3.2-10)

これを(3.2-5)式に代入して電場 E について整理する。荷電準中性条件より ene = Zeni を考慮し次式を得る。
E=⎫
1 ⎧1⎨ (∇ × B p )× B − J i × B − ZTe ∇ni ⎬Zeni ⎩ µ 0⎭(3.2-11)

但し、(3.2-11)式から得られる電場はプラズマが存在する領域内でのみ使用され、プラズマが存在しない真空領域での電場は後節で述べるラプラス方程式から得られる。

19

プラズマによる磁場を記述する方程式として、以下のファラデーの電磁誘導の法則を用いる。

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

3.2.4

プラズマ粒子の運動

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

x in +1 − x in
= v in +1 / 2(3.2-13)∆t
ここで各変数の上付き添え字は時間を表し、下付き添え字は i 番目のイオンを表す。また、∆t は時間幅を表す。
(3.2-1)式に時間中心差分を施すと以下の式となる。
v in +1 / 2 − v in −1 / 2 Ze ⎧ n 1 n +1 / 2⎫=+ v in −1 / 2 × B n ⎬(3.2-14)⎨E + v i
∆t
mi ⎩
2
⎭

(

)

n − 1 / 2 はある時刻 n から ∆t / 2 だけ以前の値であることを示す。同様に、 n + 1 / 2 は ∆t / 2 だ

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

20

(3.2-14)式を変形すると次式となる。

Av n +1 / 2 = S(3.2-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.2-16)(3.2-17)

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

v n +1 / 2 = A −1S(3.2-18)−1

ここで、 A は以下のように得られる。

⎛ 1+α 21⎜
A −1 = ⎜ αβ − γδ⎜⎝αγ + β3.2.5

αβ + γ αγ − β ⎞⎟1 + β 2 βγ + α ⎟βγ − α 1 + γ 2 ⎟⎠(3.2-20)

イオン数密度、電流密度

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

21

δni (i, j , k ) =V2
NV(3.2-21)

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

J i = Zeni v i(3.2-22)

Fig.3-5 PIC 法

3.2.6

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

B n +1 / 2 = B n −∆t∇ × En2(3.2-23)

また、(3.2-11)式は以下の式になる。
En =

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

(

)

22

(3.2-24)

ここで、 B n 、 E n は既知であるとすると、 (3.2-23) 式から B n +1 / 2 を求めることが出来る。次に、

B n +1 / 2 、 J n +1 / 2 、 v n +1 / 2 、 nin +1 / 2 を(3.2-24)式に代入して E n +1 / 2 を求め、線形外挿より電磁場の予+1+1、 B npredを次式より求める。
測値 E npred+1
E npred= −E n + 2E n +1 / 2+1= B n +1 / 2 −B npred(3.2-25)∆t+1∇ × E npred2(3.2-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.2-27)+3 / 2
そして、(3.2-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.2-28)∆t∇ × E n +1
2

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

ここで q は収束計算の回数を表し、 i 、 j 、 k はそれぞれ x、y および z 方向の格子番号を表している。 ε は収束計算を打ち切るための値である。
電場を計算するために用いる(3.2-24)式は、分母に数密度を含んでいるため、 ni がゼロとな

23

る領域ではこの式は発散してしまう。さらに、数値計算の特性上 ni があまり小さい値になると、発散もしくは計算誤差が生じる。このことを避けるために、数密度がある特定の値 nv 以上の場
合においてのみこのような電場の計算を行う。

(2) 真空領域での電場の計算
真空領域では、ラプラス方程式 ∇ 2 E = 0 を SOR（Successive Over Relaxation）法を用いて解く。
SOR における加速係数 α は正方形領域の場合は次のように表される。
2(3.2-31)α=
1 + sin (π NGMX )ここで 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.2-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.2-33)− 6 E p (i, j , k ) , ( p = x, y, z )3.2.7

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

M + 6 S + 12 K

24

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

14

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

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

3.2.8 Hybrid code に課せられる制限
ハイブリッドコードにおいて、差分化の時間幅 ∆t 、空間幅 ∆x は次の条件を満たす必要がある。
CFL（Courant-Friedrichs-Lewy）条件∆x∆t <(3.2-38)
VA
これは、 ∆t , ∆x → 0 で解が収束する条件である。 V A はアルフベン速度であり次式で定義する。
B0
VA ≡(3.2-39)µ 0 ni mi
また、付加的な数値的安定条件は以下である。
ω ci ∆t < 0.2(3.2-40)

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

∆x > c ω pi

(3.2-41)

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

25

第4章
4.1

計算結果

評価方法

今回、ターゲットである燃料ペレットを整形した場合とそうでない場合のシミュレーション結果を評価する指標として、SPH 法を用いたシミュレーションコード(以下 SPH code)の計算結果による全プラズマ粒子の運動量の角度分布と、3D hybrid code の計算結果による推進効率ηを用いた。
推進効率 η は以下の式で求めた。

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

(4.1)

0

つまり、広がるプラズマ粒子全体の+Z 方向の運動量を増加させると推進効率は高くなる。
それぞれのターゲットに対するシミュレーション結果を運動量の角度分布、推進効率について比較し、それについて考察した。

4.2

球状ターゲット

4.2.1 計算体系
SPH code の計算体系を Fig.4-1 に、3D hybrid code の計算体系を Fig.4-2 に示す。また、それぞれの計算に用いたパラメータを Table4-1、Table4-2 に示す。
Fig.4-1 において、中心の核融合ペレットの半径は 0.5cm である。そして 0.5cm の空隙を設け、そのまわりに厚さ 1.0cm の推進剤(水素)を配置している。核融合ペレットをなす粒子数は 1 万個であり、推進剤をなす粒子の個数は十万個である。どちらも一様に分布させてある。
また Fig4-2 においてコイルの中心は Z 軸上の原点から-1m に有り、半径は 1m、X-Y 平面に平行である。プラズマ粒子は SPH code の計算で得られた位置、速度を持っている。計算領域
は Z 軸方向に-7m〜7m、X,Y 軸方向にそれぞれ-6m〜6m の直方体とした。

Fig.4-1 球状ターゲット(概念図)26

Fig.4-2 磁気ノズル内のプラズマ配置(概念図)

Table 4-1 SPH code の各種パラメータ

Pellet mass[mg]0.141Moderator mass[mg]100Pellet density[kg/m3]

0.2696Moderator density[kg/ m3]

4.548

Time step ∆t [nsec.]

0.1

Number of particles

110000

Table 4-2 3D hybrid code の各種パラメータ

Coil radius[m]1.0Coil current[A]3.57 × 10 6

Coil position along Z[m]-1.0

Plasma coordinate[m](0,0,0)

Plasma mass[mg]100Time step ∆t [nsec.]

0.277

Number of particles

10000027

4.2.2 計算結果
SPH code の計算結果である。右の図は XZ 平面での断面図である。

Fig.4-3 球状ターゲットの膨張過程(上から 0.1,180[nsec.])

球状ターゲットの運動量の角度分布である。横軸は XY 平面と Z 軸からなる角度(概念図をFig.4-4 に示す)、縦軸はその角度に入る粒子の運動量の総和である。
0.08
0.07

運動量 [kg ・m/s]

0.06
0.05
0.04
0.03
0.02
0.01

96
10
4
11
2
12
0
12
8
13
6
14
4
15
2
16
0
16
8
17
6

88

72

80

64

48

56

40

24

32

16

0

8

0
XY平面とZ軸がなす角度 [度]Fig. 4-4

概念図

Fig.4-5 運動量の角度分布28

3D hybrid code の計算結果である。

Fig.4-6 磁気ノズル内での膨張過程

29

4.3

整形ターゲット

4.3.1 計算体系
SPH code の計算体系を Fig.4-7 に、3D hybrid code の計算体系を Fig.4-8 に示す。また、それぞれの計算に用いたパラメータを Table4-3、Table4-4 に示す。
Fig.4-7 において、中心の核融合ペレットは球状ターゲットの場合と同じで半径 0.5cm、空隙
も 0.5cm で同じである。しかし推進剤を図のように整形し、部分的に推進剤の厚さを薄くした。これは核融合ペレットが核融合反応を起こした際に生じる衝撃波が推進剤の薄いところから先に宇宙空間に抜けていき、最終的に Fig.4-8 のように発生したプラズマの形状が円錐のようになり、より多くのプラズマ粒子をコイル方向へ膨張することを狙ったものである。

Fig.4-7 整形ターゲット(概念図)Fig.4-8 磁気ノズル内のプラズマ配置(概念図)30

Table 4-3 SPH code の各種パラメータ

Pellet mass[mg]0.141Moderator mass[mg]100Pellet density[kg/m3]

0.2696Moderator density[kg/ m3]

3.411

Time step ∆t [nsec.]

0.1

Number of particles

110000

Table 4-2 3D hybrid code の各種パラメータ

Coil radius[m]1.0Coil current[A]3.57 × 10 6

Coil position along Z[m]-1.0

Plasma coordinate[m](0,0,0)

Plasma mass[mg]100Time step ∆t [nsec.]

0.277

Number of particles

10000031

4.3.2 計算結果
SPH code の計算結果である。

Fig.4-9 整形ターゲットの膨張過程 右は XZ 平面での断面図(上から 0.1,180,300[nsec.])

32

整形ターゲットの運動量の角度分布である。横軸は XY 平面と Z 軸からなる角度、縦軸はその角度に入る粒子の運動量の総和である。
0.1
0.09
0.08

運動量 [kg ・m/s]

0.07
0.06
0.05
0.04
0.03
0.02
0.01

96
10
4
11
2
12
0
12
8
13
6
14
4
15
2
16
0
16
8
17
6

88

72

80

64

48

56

32

40

16

24

0

8

0
XY平面とZ軸がなす角度 [度]運動量の角度分布

Fig.4-10

一定角度内にある粒子の速度の平均値である。横軸は XY 平面と Z 軸からなる角度、縦軸は粒子の平均速度。
60000

50000

30000

20000

10000

10
4
11
2
12
0
12
8
13
6
14
4
15
2
16
0
16
8
17
6

88

96

80

72

64

56

48

40

32

16

24

0

0

8

平均速度 [m/s]40000

XY平面とZ軸からなる角度 [m/s]Fig.4-11粒子の平均速度33

3D hybrid code の計算結果である。

Fig.4-12

磁気ノズル内での膨張過程

34

第5章

結論

今回、ターゲットが球状の場合と整形させた場合の SPH code によるプラズマと推進剤の衝突を解析し、それで得られた衝突後の粒子の位置と速度データを 3D hybrid code に用い、それぞれのターゲットでの推進効率を求めた。
SPH code の解析の結果、球状ターゲットの場合は、Fig.4-3 から等方的に膨張している様子がわかる。Fig.4-5 の運動量の角度分布図でも角度によって大きな差は見られない。また 3D
hybrid code で計算させた結果、推進効率は 58.4%であった。
整形ターゲットの場合は、推進剤を薄くした部分から衝撃波が抜け、プラズマが噴出している。
本研究の目的は、プラズマの噴出方向を変えることで性能の向上を目指すというものだった。しかし、Fig.4-10 の運動量の角度分布図を見ると角度による運動量の差が見られない。これは、Fig.4-11 のようにプラズマが噴出している部分の粒子の速度は速く、その他の粒子の速度は遅いが、速度が速い部分の粒子数は少なく、速度が遅い部分の粒子数が多くなっている。そのため運動量の角度分布に差が生じなかったと考えられる。また速度が速いプラズマが膨張している部分は Fig.4-12 からコイル方向へ膨張し Z 軸正の方向に偏向しているのがわかる。
以前、中島らにより整形ターゲットの設計がなされた(13)。Fig.5-1 に示す。今回使用した整形ターゲットはこれをモデルにしている。この時、膨張のシミュレーションも行われている。Fig.5-2にその図を示す。

Fig.5-1 整形ターゲット
Fig.5-2 膨張後の粒子位置(66.2[nsec.])
プラズマ粒子が推進剤の薄いところから噴出している様子がわかる。

Fig.5-3 運動量の角度分布35Fig.5-3 は運動量の角度分布図である。実線が推進剤が球状の場合、破線が推進剤を整形した場合だが、推進剤が球状の場合よりも整形させた場合のほうが部分的に運動量を集中できている。本研究で同様の結果を得ることができなかったのは、SPH code のパラメータを正しく設定できていなかった可能性がある。
また、推進効率は 71.1%であった。運動量の角度による差はできなかったが、速度が速いプ
ラズマ粒子の噴出方向をコイル方向に向けることができた。その粒子が磁場により Z 軸正の方向に偏向され、球状ターゲットの場合より同方向に膨張した粒子数が増加したため推進効率が向上したと考えられる。
今後の計画として SPH code の改良を考えている。運動量を集中させることができれば、より正確な推進効率を求めることができる。また、今回のモデルとは違う形状の整形ターゲットを考えている。例として Fig.5-4 に概念図を示す。

Fig.5-4 概念図

点火点で炸薬が 爆轟を始め、爆轟波(衝撃波)が金属ライナに到達すると材料が崩壊し、高速金属微粒子のジェットが生成し始める。これを応用し、金属ライナの部分に推進剤を配置し、点火点で核融合反応を起こし、解析を行う。

36

参考文献
(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)川渕 亮,“磁気ノズルを用いたレーザー核融合ロケットの推力ベクトルの制御技術” H18.2(7) C. D. Orth: VISTA - A Vehicle for Interplanetary Space Transport Application Powered by
Inertial Confinement Fusion, UCRL-LR-110500 (2003)(8) G. R. Liu , M. B. Liu, “ Smoothed Particle Hydrodynamics a meshfree particle method ”
2003 by World Scientific Publishing Co. Pte. Ltd.
(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) H. Nakashima, et al., “Some Design Aspects of Laser Fusion Rocket”, PROCEEDING OF
THE SIXTEENTH INTERNATIONAL SYMPOSIUM ON SPACE THECHNOLOGYAND SCIENS SAPPORO(1988)37

謝辞
本年度の卒業研究だけでなく三年生の課題集約からご指導頂き中島秀紀教授に深く感謝致します。卒論の書き方等、様々な助言をして頂いた渡辺幸信助教授に深く感謝致します。計算に関することをはじめ、様々なことをご指導頂いた梶村好宏講師に深く感謝致します。研究室
における些細なことから研究に関することまで様々な質問に答えて下さった山本直嗣助手に深く感謝致します。楽しく快適な研究室生活を与えて下さった岡安悦子秘書に深く感謝致します。
先輩としての貴重な意見を下さった桑野寛久氏に深く感謝致します。パソコンに関する様々なことを教えて下さった林真照氏に深く感謝致します。
楽しい研究室生活にして頂いた小川真司氏、加来大輔氏、笠岳幸氏、近岡貴行氏、野田賢
治氏、宗貞伸孝氏、金川隆保氏、川渕亮氏、近藤慎哉氏、永田英隆氏、西嶋康太氏、三好誠
氏、大宅将史氏、本岡親英氏に深く感謝致します。
特に川渕亮氏には同じ核融合ロケットグループとして計算に関する様々なことを自分の研究
時間を削ってまで教えて頂いた上に、時折、生意気なことを言う私にあきれることなく最後まで付き合って頂き深く感謝致します。また野田賢治氏も 3D hybrid code など研究に関する助言を与えて頂き深く感謝致します。笠岳幸氏と西嶋康太氏には席が隣であったため何度も何度も質問をしてしまいましたが、その度に親身に答えて下さり深く感謝致します。
最後に、数々の悩みを聞いてくれた友人達、そして今まで私を支えてくれた家族に深く感謝するとともに、お礼の言葉を述べて謝辞とさせて頂きます。
「ありがとう。」
平成 19 年 2 月 15 日

38