目次
第 1 章 序論 ...................................................................................................................... 1
1.1 背景 ......................................................................................................................... 1
1.2 目的 ......................................................................................................................... 3
第 2 章 レーザー核融合推進システム ............................................................................ 5
2.1 推進原理 ................................................................................................................. 5
2.2 レーザー核融合ロケット .......................................................................................... 6
第 3 章 シミュレーションコード .......................................................................................... 7
3.1 SPH 法について...................................................................................................... 7
3.1.1 衝突解析における数値解析法........................................................................ 7
3.1.2 SPH 法の概念 .................................................................................................. 8
3.1.3 SPH 法の解析手法 .......................................................................................... 9
3.1.4 フローチャート ................................................................................................ 19
3.2 3D hybrid code について ...................................................................................... 20
3.2.1 はじめに.......................................................................................................... 20
3.2.2 基礎方程式 .................................................................................................... 20
3.2.3 計算方法 ........................................................................................................ 22
3.2.4 プラズマ粒子の運動 ...................................................................................... 22
3.2.5 イオン数密度、電流密度 ............................................................................... 24
3.2.6 電磁場の計算方法 ........................................................................................ 25
3.2.7 Digital Filter ................................................................................................... 27
3.2.8 Hybrid code に課せられる制限 ..................................................................... 28
第 4 章 計算結果 ............................................................................................................ 29
4.1 評価方法 ............................................................................................................... 29
4.2 整形ターゲットの概念 ........................................................................................... 29
4.3 球状ターゲットと整形ターゲットの比較 ................................................................ 30
4.4 方向転換 ............................................................................................................... 34
4.4.1 ターゲットの回転 ............................................................................................ 35
4.4.2 Cut モデル ...................................................................................................... 38
4.5 Egg モデル ............................................................................................................ 40
4.6 衝撃点火 ............................................................................................................... 45
4.6.1 衝撃点火ターゲット ........................................................................................ 46
4.6.2 衝撃点火ターゲットの変形 ............................................................................ 52
第 5 章 結論 .................................................................................................................... 54

第1章
1.1

序論

背景

1961 年，人類は初めて宇宙へと進出し，1969 年には月へも降り立った．今年はアポロ
11 号が月面着陸を果たしてからちょうど 40 年になる．
現在も様々な形で宇宙空間の利用，宇宙開発は進められている．気象衛星，GPS 衛星など我々の生活に身近なものや，着々と完成に近付いている国際宇宙ステーションな
どでは宇宙空間，無重力下における様々な実験を行っている．また，SOHO やハッブル
宇宙望遠鏡など，太陽や惑星，外宇宙の観測など太陽系や宇宙全体の謎の解明に向けた ミ ッ シ ョ ン も 行 わ れ て い る ． 日 本 に お い て は ， 「 か ぐ や （ SELENE ） 」 ， 「 は や ぶ さ（MUSES-C）」といった探査機や衛星を打ち上げ，つい先日も温室効果ガスの観測を目的とした「いぶき（GOSAT）」の打ち上げに成功している．また最近では，アメリカ航空宇宙局（NASA）のスペースシャトル退役に伴い，国際宇宙ステーションの物資補給等のミッシ
ョンを民間のロケットに委託したり，日本でも，大阪の中小企業が中心となって作った衛星「まいど 1 号（SOHLA-1）」が打ち上がるなど，宇宙開発は国家主導のものだけではなく，
民間レベルに広く浸透してきたと言える．
またこれらのロケット，探査機に使用されている推進装置は化学推進と電気推進である．
化学推進は大きな推力を出すことができるが，すぐに燃料を使い果たしてしまう．電気推進は燃料利用効率が高く，長時間使用することができるが推力が小さい．ロケットの打ち上げ時には短時間で大きな推力が得られる化学ロケットを使用し，宇宙空間において人工衛星の推進，姿勢制御のように小さな推力で済むものには電気推進を使用し，それぞれの役割で使い分けられている．
しかし，宇宙開発の次なるステップとして挙げられる有人惑星飛行の場合は，無重力による人体への影響，宇宙放射線による被曝を考慮してできるだけ短い時間でミッションを終わらせる必要があり，食料や探査に必要な物資が大量になってしまう．そのため有人
惑星間飛行に使用されるロケットには大きな推力を得られ，燃費がいいエンジンが必要になると考えられ，熱核ロケットや，核パルス推進，磁気セイルなど様々なタイプの推進方法が提案，研究されている．
この研究の中で，磁気ノズルというプラズマの制御技術が考えられている．磁気ノズルについては 2 章で詳しく述べる．磁気ノズルを用いると，膨張するプラズマを磁場との相互作用を利用し，固体壁との相互作用無く運動方向を変えて推力を得ることが出来る．この技術を用いたロケットとして NASA で開発が進められている VASIMR（Fig.1-1）や，レー
ザー核融合ロケットが考えられており，他の宇宙推進システムと比較して高い排出速度（即ち高い推進材利用効率）と大きな推力を同時に達成可能である．レーザー核融合ロケットの概念図を Fig.1-2 に示す．

1

Fig. 1-1 VASIMR

Fig. 1-2

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

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

1.2

目的

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

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

下段：膨張後この計算により，それぞれの推進効率は球状の場合が 69 %，円錐の場合が 73 %となっており，運動量を集中させることで推進効率の上昇をシミュレーションできた．
このシミュレーションはターゲットが膨張して円錐の形になったものとして計算させている．そこで本研究では，Smoothed Particle Hydrodynamics(SPH)法を用いて核融合プラズ
マと推進剤の衝突を計算させ，衝突後の粒子の位置、速度を求める．そのデータを 3D
hybrid code の粒子の位置，速度データとして入力し，磁気ノズル内でのプラズマの挙動を計算し，推進効率を求める．
本研究の目的は，推進性能を向上すべく新たな方式を検討する．それは，整形ターゲットを使用することである．核融合燃料ペレットの周囲に推進剤(Moderator)を配置し，レーザー核融合で発生したプラズマと推進剤の衝突により発生するプラズマの運動量を磁場発生用コイル方向(磁場の強い方向)に集中させることで推進効率の向上を目指す．
また，様々な形状を使用して磁気ノズル内でのプラズマの挙動を解析し，整形ターゲッ
トを用いることによるレーザー核融合ロケットの推進性能の向上の可能性を探る．
3

膨張過程

Fig. 1-4

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

4

第2章
2.1

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

推進原理

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

Fig. 2-1

磁気ノズルでの推進原理
52.2

レーザー核融合ロケット

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

6

第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 分布関数等であらかじめ決定しておいた平滑化関数を用いる．

7

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 法を基に空間に配置した粒子のイメージ

8

粒子” 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 は粒子数である．

9n xi     xi  x j N(3.1-5)j 1

これを粒子 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 流体の密度の概念図

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

2)

3)

D
v 
  
Dt
x(3.1-9)Dv 
1  Dt x 

(3.1-10)

運動方程式

エネルギー方程式
De   v Dt x 

(3.1-11)上の式で  は全応力テンソルである．これは二つの部分からなり，一つは等方性の圧力
p で，もう一つは粘性応力  である．
    p    
(3.1-12)ニュートン流体において粘性せん断応力は動的粘性係数μを介してεで表されるせん
断歪率に比例しなければならない．
    (3.1-13)
ここでv  v  2
   (  v) (3.1-14)3
x
x
等方性の圧力と粘性応力が分割し考える場合，エネルギー方程式は次のように書き換え  

られる．
Dep v   


 

Dt
 x
2

11

(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
 i v j  Dtxij 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 jj 1  jN

(3.1-19)

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

mj

j 1

j

i 

vi 

 N m j Wij   i vi   j 1  x  
xi
j
i



Wij

(3.1-20)式(3.1-20)を式(3.1-18)に代入すると式(3.1-21)を得る．
N m
Wij
D i
j
 i 
vij  Dtxij 1  j

ここで， vij  vi  v j である．
12

(3.1-21)

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

 v     
v 



 v   

x 
x 
 x

(3.1-22)

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

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

 
Wij
Dvi
1 N
j

mj

Dt
 i j 1 j xi

(3.1-24)

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

N

mj
j 1

 i Wij  i  N m j Wij 0 i  j xi
 i  j 1  j xi N i   
Wij
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)13

簡単な再配列の後，式(3.1-28)は次のように書き換えられる．
N    
 WijDv ij mj i2  2   Dt j  xi
j 1
 i

(3.1-29)式(3.1-12)を用いると，式(3.1-26)と式(3.1-29)はそれぞれ次のように書き換えられる．
NNpi  p j Wij
 i  i   j  
Wij
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
Dv i
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




v j

 2 N mj
 
 


v


W
j
i
ij

xi  3 j 1  jWij(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)

i

j 1

j

N

mj

j 1

j


N

i

i

vi

 j 1

j

i



Wij

 N mjmj  v   W  v      W   0
j 1

i

i

ij

i

j

 j 1

i

j

ij



(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)により加速度を計算することがでる．

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

p v 
p 
v  
p D

 2




2 
 
 x
 
x   Dt

式(3-37)の関係を考慮すると，関係を持つ粒子でp2(3.1-37)の値を求めながら，圧力仕事は式(3.1-18)や式(3.1-21)．式(3.1-23)の密度の連続性を用いることにより近似できる．例えば，
式(3.1-23)の密度の連続性を用いると，圧力仕事は次式のように表される．
Wij
pi N
p vi



m
v


j
ij
 xi  i2 j 1
xi

(3.1-38)

別の方法は次の恒等式を用いる．

p v   pv      p 

v
 


 x 
x    
x    

(3.1-39)

すると次の式を得る．



N
pj
Wij
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 
xi

関係を持つ粒子で求めたp

2(3.1-41)の値と式(3.1-21)を用いることで，式(3.1-37)の時間の増減率について近似すると次式を得る．
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)15また式(3.1-43)の変形を用いることで圧力仕事についての SPH 法の近似の別の形が求まる．
p j  Wijp vi
1 Nmvij  j
 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mvij  j
 xi 2 j 1
i  jxi(3.1-45)

圧力仕事を近似するには多くの形式があるので．内部エネルギーの発生のための SPH
法の公式は自然と多くの代替手段を持つことができる．
もっともよく使われるのは次の二式である．
 p
pj 
Wij
ei 1 N  m j  2i  2 vij    i  i  it
2 j 1   i
2 i
 j 
xi(3.1-46)pi  p j  Wijei 1 N mjvij    i  i  it
2 j 1i  j2 ixi(3.1-47)

また，密度，運動量，およびエネルギーを導くためのオイラー方程式に対応する SPH 法のそれぞれの方程式は以下のとおりになる．
連続の式
N
Wij
D i
  m j vij  
Dt
xi
j 1

(3.1-23)Npi  p j WijDvi  m jDt i  j xi
j 1

(3.1-48)N p
p j  WijDv i
  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 N mjvij  t
2 j 1i  jxi(3.1-51)16d) 人工粘性
SPH では人工粘性を考慮している。レーザー核融合で生じたプラズマと推進剤との衝突を SPH 法で記述しようとすると Fig.3-3 のように相対速度が音速を超える粒子群のぶつ
かり合いとして表される．

Fig.3-3

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

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

    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 はそれぞれ音速と粒子の速度ベクトルを表している．

17

この人工粘性を 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

18

(3.1-58)

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

近傍粒子の検索

平滑化計算

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

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

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

データ書き出し(OUTPUT)

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

19

3.2

3D hybrid code について

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

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

mi

dv i
 ZeE  v i  B 
dt(3.2-1)

dx i
 vi
dt

(3.2-2)

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

n e medv e
 ene E  v e  B   Pe
dt(3.2-3)

ここで v e は電子の速度， m e は電子の質量， n e は電子数密度， Pe は電子圧力である．
電子は慣性無しの流体としてみなすので，慣性による効果は無視し，(3.2-3)式の左辺は0 にできる( me  0 )．さらに電子圧力 Pe に関しては電子の状態方程式より以下のようにな
20

る．
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  Zen iを考慮し次式を得る．

E

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

21(3.2-11)但し，(3.2-11)式から得られる電場はプラズマが存在する領域内でのみ使用され，プラズマが存在しない真空領域での電場は後節で述べるラプラス方程式から得られる．
プラズマによる磁場を記述する方程式として，以下のファラデーの電磁誘導の法則を用いる．
B   E
t(3.2-12)3.2.3 計算方法
(3.2-1)式から(3.2-12)式を数値計算によって解くためには，基礎式を差分式に直し，数
値的に安定で，物理的に妥当な空間幅、時間幅で解かなければならない．
計算の大まかな手順を以下に示す．
[1]外部磁場および，イオンの粒子位置，速度の初期設定を行う．
[2]格子点上でのイオンの数密度，電流密度を計算する．
[3](3.11)式により，電場を求める．
[4][1]で設定した外部磁場と，[3]で求めたプラズマによる電場と磁場から，(3.2-1)式，(3.2-2)式を用いて新たなイオン粒子の位置、速度を求める。
[5](3.12)式により磁場を求める．
[2]～[5]を繰り返し，プラズマの時間的変化を追う．
３次元ハイブリッドコードに対する各方程式の差分式を以下の小節で示す．

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 
E  v i
t
mi 
2






(3.2-14)n  1 / 2 はある時刻 n から t / 2 だけ以前の値であることを示す．同様に， n  1 / 2 は
t / 2 だけ進んだ値であることを示す．上付き添え字が n ， n  1 / 2 である変数は既知量である．

22

プラズマ粒子の運動を求めるためには，まず(3.2-14)式において， v in 1 / 2 を左辺にまとめて，未知量である v in 1 / 2 を解く．次に，求めた v in 1 / 2 を(3.2-13)式の右辺に代入し，新たな粒子の位置 x in 1 を求めることが可能となる．
(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             
1   2       1   2 

23

(3.2-20)

3.2.5

イオン数密度、電流密度

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

ni (i, j , k ) V2
NV(3.2-21)

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

Fig.3-5 PIC 法

24(3.2-22)

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

B n1 / 2  B n (3.2-23)

また，(3.2-11)式は以下の式になる．

En 

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

(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の予測値 Enpred、 B npredを次式より求める．

1
Enpred En  2En1 / 2

1
B npred B n1 / 2 t
1  E npred2(3.2-25)(3.2-26)

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

3/ 2
3/ 2
3/ 2
3 / 2
れらによるイオンの運動を解き v inpred， ninpred， J inpredを求める．次に， B npred
を次式により求める．
3 / 2
1
B npred B npredt  E n1 / 2
2(3.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 タイムステップの値を求めるため
25

に， 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)式は，分母に数密度を含んでいるため， n i がゼロとなる領域ではこの式は発散してしまう．さらに，数値計算の特性上 n i があまり小さい値になると，発散もしくは計算誤差が生じる．このことを避けるために，数密度がある特定の値 n v 以上の場合においてのみこのような電場の計算を行う．

(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 が、あらかじめ決められた値以下になれば真空部分の電場の計算を終了する．

26

 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

(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)式は次式となる．
1
4

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

1 side1
 i , j ,k   icorner, j ,k12
48(3.2-37)

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

27

3.2.8 Hybrid code に課せられる制限
ハイブリッドコードにおいて，差分化の時間幅 t ，空間幅 x は次の条件を満たす必要がある．
CFL（Courant-Friedrichs-Lewy）条件
x(3.2-38)t 
VA
これは， t , x  0 で解が収束する条件である． V A はアルフベン速度であり次式で定義
する．
VA 

B0
 0 n i mi

また，付加的な数値的安定条件は以下である．
 ci t  0.2(3.2-39)(3.2-40)

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

(3.2-41)

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

28

第4章
4.1

計算結果

評価方法

シミュレーション結果を評価する指標として，3D hybrid code の計算結果による推進効率ηを用いた．
推進効率 η は以下の式で求めた．

推進効率  mv 全Z方向運動量 m v 全初期運動量
z

(4.1)

0

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

4.2

整形ターゲットの概念

①

②

③

Fig.4-1 膨張過程
整形ターゲットとは核融合ペレットの周囲に配置した推進剤を球形ではなく，ある形に整形したターゲットである．
上図は整形ターゲットが膨張する過程を示している．まずターゲットの中心に核融合ペレットを配置して周囲の推進剤を上図のように整形し，レーザーを照射し核融合反応を起こす（①）．核融合反応で生じたプラズマは周りの推進剤をプラズマ化していく．また核融合反応で生じた衝撃波が推進剤の中を伝搬し（②），推進剤が薄い部分から先に推進剤の外へ抜けて，その部分のプラズマが先に噴出し，最終的に（③）のような形にプラズマが膨張する．

29

4.3

球状ターゲットと整形ターゲットの比較

核融合ペレットの周りに配置した推進剤を整形していないものを球状ターゲット，整形したものを球状ターゲットとする．以下に球状ターゲット，整形ターゲットのそれぞれの概念
図，SPH code と 3D hybrid code での計算体系を示す．ただし，SPH code では核融合反応のシミュレーションはできないため，中心で核融合反応が起こり，プラズマが生成された状態から計算を始めている．

Fig.4-2 球状ターゲットと整形ターゲット
Table.4-1 計算パラメーター
Pellet mass[mg]0.141Moderator mass[mg]100Time step Δt [nsec.]

0.1

Number of particles

110000

Fig.4-3 膨脹後30

Fig.4-4 3D hybrid code での計算体系

Table 4-2 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

Fig.4-5 磁気ノズル中での挙動(上段：初期状態，中段：5.54μsec.後，下段：11.08μsec.後)32

球状ターゲットの場合の推進効率 η は 66 %，整形ターゲットの場合の推進効率は 78 %だった．推進効率が向上した理由としては，矢印 1 の部分の速度が大きい粒子はコイル近傍の強い磁場によってZ 方向に偏向され，矢印 2 の部分は速度が小さいので磁場によって容易にZ 方向に偏向されたためだと考えられる．
このことより，整形ターゲットを用いることで，プラズマの膨張をある程度制御できることが
わかった．

1

2

60000

平均速度 [m/s]

50000
40000
30000
20000
10000
0
4

28

52

76

100

124

148

XY平面とZ軸からなる角度 [度]Fig.4-6 速度分布

33

172

4.4

方向転換

ロケットの姿勢制御や方向転換を行う際には二通りの方法が考えられる．一つはメインエンジンとは別に，イオンスラスタ等の小型エンジンを取り付け，それを使用することで推進方向を制御する方法である(Fig.4-7)．もう一つの方法はメインエンジンにおいて，燃料の排出方向を傾けることで推進方向を制御する方法である(Fig.4-8)．前者の方法では小型エンジン本体や，その燃料，電源装置がメインエンジンとは別に必要があり，船体が重くなってしまう．しかし，プラズマの排出方向を非軸対称にすれば steering angle が発生し，
メインエンジンだけでの方向転換が可能になる．
その場合，前述の通り船体の軽量化につながるので，その分を物資の運搬に利用できるなど効率的にロケットを運用できる．
以前行われた計算では(6)，二つの超伝導コイルを使用し，片方を動かすことで磁気ノズルの形状を変え，プラズマの排出方向を非軸対称にしていた．しかし，それではコイル間の相互作用を支える機構，コイルを動かす大がかりな機構が必要になってしまい，船体の重量が増加してしまう．
そこで，整形ターゲットを用いることでの方向制御の可能性を探る．

Fig. 4-7 従来の方法での方向転換

Fig.4-8 排出ベクトルの制御を用いた方向転換
34

4.4.1 ターゲットの回転4.3 で使用したターゲットは回転させることでプラズマの膨張を非軸対称にすることができる．これによりプラズマの排出方向も非軸対称になり steering angle が発生する．

Fig.4-9 3D hybrid code での計算体系

Table 4-3 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

100000α [deg.]0~180

4.3 で使用した整形ターゲットのケースを 0 度として，α (Fig.4-9 参照)を 0 度から 180 度
まで 30 度ずつ変化させ計算を行い，それぞれの角度における推進効率 η，Steering
angle β を求めた．各種パラメーターは 4.3 のケースと同じである．
また Steering angle β は排出ベクトルと Z 軸のなす角度である．

35

Fig.4-10 α = 120°での計算結果(上段：初期状態，中段：5.54μsec.後，下段：11.08μsec.後)36

α = 30°の時，推進効率 η は 62 %で，steering angle β は 0.027 度だった．
α = 60°の時，推進効率 η は 58 %で，steering angle β は 0.35 度だった．
α = 90°の時，推進効率 η は 58 %で，steering angle β は 4.8 度だった．
α = 120°の時，推進効率 η は 66 %で，steering angle β は 4.9 度だった．
α = 150°の時，推進効率 η は 69 %で，steering angle β は 1.1 度だった．
α = 180°の時，推進効率 η は 60 %で，steering angle β は 0.34 度だった．
角度 α と推進効率 η，steering angle β の関係を Fig.4-11 に示す．

Fig.4-11 α と推進効率 η，Steering angle β の関係

上のグラフからもわかるように推進効率 η は角度 α が 0 度のケースが一番高かった．
Steering angle β は，α が 90°と 120°のケースが高い値を得たが，推進効率とのトレード
オフを考えると α が 120°のケースが最も方向転換に適していると考えられる．
Fig.4-10 のプラズマの挙動を見ると，多くの粒子が－X 方向に広がっており，－X 方向の速度成分を持つ粒子が多いことがわかる．このことにより全体的な排出ベクトルが－X方向に向かい，Steering angle が発生したと考えられる．

37

4.4.2 Cut モデル
次に Cut モデルを示す．Cut モデルは Fig.4-12 に示すとおり，球状ターゲットの一部分を切り落としたものである．切り落とした部分は他の部分よりも推進剤が薄くなるので 4.3で示したとおり，その部分のプラズマの膨張速度が大きくなり，軸対称性が失われる．

Fig.4-12 Cut モデルの概念図
Table 4-4 計算パラメーター
Pellet mass[mg]0.141Moderator mass[mg]100Time step Δt [nsec.]

0.1

Number of particles

110000

Table 4-5 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

10000038

Fig.4-13 計算結果
(上段：初期状態，中段：5.54μsec.後，下段：11.08μsec.後)推進効率 η は 67 %，Steering angle β は 3.96 度だった．

39

4.5

Egg モデル

これまでコイルの内側を通り－Z 方向に抜ける粒子を尐なくするため，磁場配位を変化させるなどの解析が行われた．しかし，コイルを二つ以上使用した場合，前述の通り，コイル間の相互作用を支える機構等で船体の質量が増加してしまう．
そこで，1 つのコイルでコイルの内側を通り－Z 方向に抜ける粒子を尐なくするために推進剤を整形する方法を考えた．
4.2 等で推進剤が厚い部分の粒子の膨張速度は小さいことがわかっているので，Fig.4-14 のようにコイル中心方向の推進剤を他の部分よりも厚くして，粒子の膨張速度を抑え，磁場の偏向による－Z 方向に抜ける粒子の減尐を狙った．

Fig.4-14 Egg ターゲットの概念図

Table 4-6 計算パラメーター
Pellet mass[mg]0.141Moderator mass[mg]100Time step Δt [nsec.]

0.1

Number of particles

11000040

Table 4-7 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

100000

Fig.4-15 計算結果
(上段左：初期状態，上段右：5.54μsec.後，下段左：11.08μsec.後，下段右：16.62μsec.後)

41

66 %
推進効率 [%](球状ターゲット)

80
70
60
50
40
30
20
10
0

経過時間 [µsec.]

Fig.4-16 推進効率の時間変化
これまでのケースでは，計算が 11.08μsec.経過する頃には推進効率はほぼ定常となり，それ以降上昇することはなかった．Fig.4-16 からもわかるように Egg モデルの場合も約8μsec.経過した頃から推進効率は一定になり始めた．しかし 10μsec.が経過した頃からまた推進効率は上昇し始めている．11.08μsec.後の計算を再び行ったところ，最終的に推進効率は 72 %まで上昇した．

Fig.4-17 球状ターゲットと Egg モデルのプラズマ挙動の比較(経過時間：11.08μsec.)
Fig.4-17 を見ると，球状ターゲットと Egg モデルでは 11.08μsec.が経過したプラズマの挙動に違いが見られる．Egg モデルのケースでは囲みの部分に多くの粒子が存在することを確認できるが球状ターゲットのケースではそのような粒子は存在しない．この後，これらの粒子がZ 方向に偏向されているのが Fig.4-15 の 16.62μsec.が経過した後の図からわかる．また，－Z 方向に抜けた粒子数を数えたところ，4.3 %だった(11.08 µsec.の時点では 2.6 %，球状ターゲットは 11.08 µsec.の時点で約 6 %)．これらのことが推進効率の向上につながったと考えられる．
42

次に，Fig.4-2 の整形ターゲットのように，側面を切り落して計算を行った．

Fig.4-18 計算体系

Table 4-8 計算パラメーター
Pellet mass[mg]0.141Moderator mass[mg]100Time step Δt [nsec.]

0.1

Number of particles

110000

Table 4-9 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

10000043

Fig.4-19 3D hybrid code での計算結果(上段左：初期状態，上段右：5.54μsec.後，下段左：11.08μsec.後，下段右：16.62μsec.後)

推進効率は 79 %だった．Fig.4-15 の経過時間 11.08μsec.の図と，今回の経過時間
11.08μsec.の図を比べると，原点に粒子が集中している様子が今回は見られない．側面を切り落したため，－Z 方向の粒子が前回のケースよりも尐なくなったためと考えられる．前回では－Z 方向の粒子は 9μsec.頃からZ 方向へ偏向し始めたが，今回は 6μsec.頃から
Z 方向へ偏向し始める様子が確認できた．また，切り落とした部分の質量を補うために，
その他の部分を尐し厚くした．それにより，粒子数が増えZ 方向の成分を持つ粒子が増加したことも一因だと考えられる．

44

4.6

衝撃点火

衝撃点火とは阪大レーザー研により提案された新たなレーザー核融合の点火方法である(13)．Fig.4-18 に衝撃点火の概念図を示す（同一の図中に初期構造と点火が起こる最終フェイズにおける燃料構造が描かれている．）．下側の主燃料部分は重水素・三重水素(DT)燃料殻にアブレーター(プラスチック)がコーティングされており，レーザーを照射し圧縮させ，最終的に上側のインパクター部分が中心で衝突し，その衝撃でインパクターの運動エネルギーが内部エネルギーに変化され，点火スポットと化す．

Fig.4-20 衝撃点火概念図
衝撃点火の利点は以下の三つが挙げられている．
(1) 主要物理がシンプルであること．
(2) 高利得設計が可能であること．
(3) 低コストであること．(PW 級のレーザーが不要)
そこで，cone guide の部分を推進剤で形成し，整形ターゲットを作り，プラズマの挙動を
調べた．

45

4.6.1 衝撃点火ターゲット
Fig.4-19 に衝撃点火の概念を取り入れた整形ターゲットの概念を示す．

Fig.4-21 衝撃点火ターゲット概念図
赤色の部分が核融合ペレット，青色の部分が推進剤により形成されている．今回使用している SPH code では爆縮の現象は計算できないため，中心で核融合反応が起こり，核融合プラズマができた状態から計算を開始する．Fig.4-20 に計算体系を示す．

Fig.4-22 SPH code での計算体系

Table 4-10 計算パラメーター
Pellet mass[mg]0.141Moderator mass[mg]100Time step Δt [nsec.]

0.1

Number of particles

11000046

Fig.4-23 SPH code での計算結果

Fig.4-24 3D hybrid code での計算体系

Table 4-11 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

10000047

Fig.4-25 計算結果(上段：初期状態，中段：5.54μsec.後，下段：11.08μsec.後)48

推進効率は 62 %だった．球状ターゲットよりも推進効率が低下した理由として，jet として噴出したプラズマ粒子の多くが，コイルの中心を通り－Z 方向へと抜けてしまったためと考えられる．
そこで，jet の噴出方向をコイル方向に向けて，コイル付近の強い磁場によりZ 方向に偏向されるように衝撃点火ターゲットを傾けた．

Fig.4-26 衝撃点火ターゲットへのレーザー照射
ここで，Fig.4-24 は衝撃点火ターゲットへのレーザーの照射を示している．レーザー核
融合ロケットは Fig.1-2 のようにコイル方向からのレーザーと，ミラーで反射させたレーザーにより核融合反応を起こす．衝撃点火はインパクターとなる上側の DT 燃料にレーザーを照射する必要があるので，Fig.4-24 のように，衝撃点火ターゲットを傾けるほうがレーザー照射は容易になる．

Fig.4-27 3D hybrid code での計算体系

49

Fig.4-28 3D hybrid code での計算結果(左：初期状態，右：11.08μsec.後)

推進効率は 61 %だった．この計算では－Z 方向への粒子があまり減尐しなかったこと，－X 方向への大きな速度成分を持つ粒子が増えたことにより推進効率が低下したと考えられる．そこで初期プラズマの位置を原点からX 方向へずらし，－Z 方向，－X 方向へ膨張する粒子の減尐を狙った．d は原点と初期プラズマの中心との距離である．

Fig.4-29 計算体系

50

Table 4-12 計算パラメーター
Coil radius [m]1.0Coil current [A]

3.57×106

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

100000d [m]0.1~0.4

Fig.4-30 X 方向に 10 cm ずらした場合の計算結果(左：初期状態，右：11.08μsec.後)
d=0.1 m の場合，推進効率は 63 %だった．
d=0.2 m の場合，推進効率は 61 %だった．
d=0.3 m の場合，推進効率は 59 %だった．
d=0.4 m の場合，推進効率は 56 %だった．
コイルの中心を通り，－Z 方向に抜ける粒子数を数えたところ，初期プラズマの中心が原点にあるケースでは全体の 4.5 %だったのに対して，X 方向に 10 cm ずらしたケースでは全体の 3.5 %だった．このことが，ターゲットを傾けただけのケースよりも推進効率が高くなったことにつながったと考えられる．
しかし，上の結果のとおり，原点と初期プラズマの中心との距離が長くなるほど推進効率は低下した．これはX 方向にプラズマをずらすと，コイル付近の強い磁場がある領域の上を通過するため，Z 方向に尐ししか偏向されなかったためと考えられる．

51

4.6.2 衝撃点火ターゲットの変形
Fig.4-26 を見ると，衝撃点火ターゲットを傾けただけのケースでは，－X 方向に膨張する粒子が，X 方向に膨張する粒子よりも多いことがわかる．そこで，Fig.4-29 のように－X
方向の推進剤を厚くし，粒子の膨張速度を抑え，Z 方向へ偏向する粒子の増加を狙っ
た．

X

Z

Fig.4-31 計算体系
Table 4-13 計算パラメーター
Pellet mass[mg]0.141Moderator mass[mg]100Time step Δt [nsec.]

0.1

Number of particles

110000Fig. 4-32 SPH code での計算結果
52

Fig.4-33 3D hybrid code での計算結果（上段左：初期状態，上段右：5.54μsec.後，下段左：11.08μsec.後，
下段右：Fig.4-25 での計算体系の 5.54μsec.後）

推進効率は 67%だった．Fig.4-25 の計算体系の 5.54μsec.後の図と今回の計算体系の5.54μsec.後の図を比較するとプラズマの膨張に大きな違いがみられる．推進剤が厚い部
分の粒子は膨張し始めるとすぐにコイル付近の磁場により，一度X 方向に偏向されたあと，再び磁場によりZ 方向に偏向されているように見える．このことにより推進効率が向上したと考えられる．

53

第5章

結論

これまでの計算結果から，整形ターゲットを用いることでプラズマの膨張をある程度制御でき，推進性能を向上させることができると言える．SPH code の初期状態において，周囲の推進剤は核融合反応で生じた中性子によりプラズマ化した粒子として速度 0 で配置し，中心の核融合反応プラズマは全方位に等方的に広がるような速度成分で配置されている．そのため，周囲の推進剤にエネルギーが分配される際に粒子が厚い部分は粒子数が多く，粒子 1 個あたりに分配されるエネルギーが尐なくなったと考えられる．
推進効率が向上する要因は，磁場により粒子が偏向された際に，いかに多くの粒子がZ 方向に偏向されたかによる．そこで，これまでの計算によりZ 方向の粒子が多くなる
要因として以下の 3 つが挙げられる．
1.…速度が大きい粒子はコイル近傍の磁場の強い方向に向けることで偏向される．
2.…速度が小さい粒子は弱い磁場でも容易に偏向される．
3.…－Z 方向に抜ける粒子を尐なくする．
1 番は Fig.4-2 の整形ターゲットの場合のように，推進剤が他の部分よりも薄い部分をつくることで速度が大きい粒子ができる．速度が大きいということは粒子が持つエネルギーが大きいので磁場にエネルギーを与えても磁場のエネルギーより小さくならず，そこで停止することなく進んでしまい，あまりZ 方向へと偏向されずにそのまま膨張してしまう．そこで，速度が大きい粒子をコイル方向に向ける．コイル近傍の磁場が持つエネルギーは大きいので，粒子のエネルギーが大きくてもZ 方向に偏向できる．
2 番は，1 番とは逆である．速度が小さいのでエネルギーも小さい．従って，磁場のエネ
ルギーが小さくてもZ 方向に偏向できる．
3 番は，－Z 方向に向ける粒子の速度を小さくすることで解決できる．コイル中心付近はコイル近傍よりも磁場が弱い．従って，今までの計算ではコイル中心から抜ける粒子数は計算体系には大きく左右されることはなかったが，Egg モデルでは－Z 方向の推進剤を厚くしたため粒子の速度が小さくなりコイル中心付近の弱い磁場でもZ 方向に偏向することができたと考えられる．
これらの条件を満たすターゲットは，Egg モデルが最も適していると考えられる．今後，Egg モデルを軸にターゲットを考案することで高い推進効率を望めると考えられる．
また，衝撃点火を用いたケースでは球状ターゲットのケースとの違いはあまり見られなかった．これは cone guide の部分の立体角が 0.02 程度であるため，球状ターゲットの結果とあまり変わらなかったと考えられる．従って，今のところ衝撃点火ターゲットの利点は衝撃点火という点火方式による利点だけであり，衝撃点火を整形ターゲットに応用し発生する
jet を活用する方法は難しいと考えられる．

54

最後に，今後の予定であるが，今年度は大阪大学レーザー研究所において本研究室のレーザー核融合ロケットグループにより，磁気ノズル内でのプラズマ挙動を観察する実験が行なわれた．そこで，永久磁石により作られた磁場中でのプラズマの挙動と，磁石が無い場合でのプラズマの挙動には明らかな違いが見られた．次の実験では，スラストスタ
ンドを設置し，磁気ノズルによる推力測定を行う予定である．

Fig.5-1 推力測定実験概念図
また，同時に Cone ターゲットの実験も行うことが考えられている．Cone ターゲットとは
Fig.5-2 のようなターゲットである(14)．核融合によって中性子が発生するが，その中性子がコイルや船体に衝突すると構造体が务化してしまう．そこで，核融合反応が起こる位置からコイル，船体が存在する角度上に推進剤を配置して中性子の減速に使用する．

Fig.5-2 Cone ターゲット概念図
そこで，その際に使用するターゲットのパラメーター設定を行った．使用する装置は大阪大学レーザー研究所の激光Ⅻ号である．レーザーのエネルギーは 6 kJ でエネルギー55

吸収を 10 %とするとペレットのエネルギーは 600 J である．以前行われた磁気ノズルの実
験(15)においては 3.5 J のエネルギーでプラズマ質量は 0.6×10-6g と設定された．よって簡単のためエネルギーを 200 倍としてプラズマ質量を 1.2×10-6g に設定した．また，これまでのシミュレーションにおいて核融合ペレットのエネルギーが 4MJ で推進剤質量を 0.1g に設定していた．よってエネルギーが 600 J の場合には推進剤質量を 1.5×10-5g に設定する．
そこで推進剤に水素を用いて Fig.5-2 のような形を作り，コイル方向の角度に推進剤を配置すると考えると，Cone 部分の大きさは約 0.7 mm 程度になる(Fig.5-3)．

Fig.5-3 ターゲット概念図

56

第6章

付録

以前，本研究室において，矩形コイルを用いたレーザー核融合ロケットの推進性能の向上を目的とした計算が行われた(16)．
Fig.A-1 の破線で囲まれた部分はコイルシステムの部分を拡大したものである．図のように小さな超伝導コイルを複数個使用し，円形になるように配置する．一つ一つの超伝導コイルに電流を流し，磁場を発生させ，合成磁場により磁気ノズルを形成する．

Fig.A-1 矩形コイルシステム概念図
単純に一つの円形の超伝導コイルを使用して磁気ノズルを形成する場合と違い，合成
磁場により形成されるので，矩形コイル一つあたりが発生する磁場は弱くて済む．つまり，矩形コイル一つに流す電流は小さくて済む．
今回，矩形コイルシステムは磁気ノズルの形状を変化させるものであるため，本文には掲載しなかったが，衝撃点火ターゲットにおいて，推進効率向上のために矩形コイルシス
テムを用いたので付録として紹介する．

Fig.A-2 矩形コイルシステムの計算体系

57

Fig.A-2 は衝撃点火ターゲットの一つのコイルによる計算体系(Fig.4-24 参照)の，コイル部分を矩形コイルシステムに変更した場合の計算体系である．
主な計算パラメーターは以下の通りである．
Table A-1
Rectangle coil size [m]Number of rectangle coil
Coil current [kA]
Plasma coordinate [m]Plasma mass [mg]Time step ⊿t [nsec.]
Number of particles

0.25×0.536
12
32(0,0,0)100
0.277
100000

Fig.A-3 矩形コイルシステムに対する 3D hybrid code での計算結果（左：初期状態，右：11.08μsec.後）
推進効率は 71 %だった．11.08μsec.が経過した図を見ると，コイル間の磁気ミラーに粒子が閉じ込められている様子がわかる．－Z 方向に抜けた粒子は全体の 0.6 %だった．このことが，推進効率が向上した理由だと考えられる．
次に，方向制御の可能性を考える．矩形コイルシステムは 12 個の矩形コイルを円形に配置したものであるため，一つ一つのコイルに流す電流値を変更することで，磁気ノズルの形状を自由に変えることができる．上側 6 個の矩形コイルには前回の計算と同じ 32 kAを流し，下側の 6 個の矩形コイルには半分の 16 kA を流す．これにより計算領域の－X方向の磁場が弱くなり，プラズマが－X 方向に多く排出され，Steering angle β が発生する．
以前行われた計算と違い(6)(4.4 参照)，矩形コイルシステムの上下で電流値を変化させるだけで磁気ノズルを非軸対称にすることができるので，ロケットを効率良く運用できると考えられる．
Fig.A-4 に計算体系を示す．
58

Fig.A-4 矩形コイルシステムの計算体系
Table A-1
Rectangle coil size [m]Number of rectangle coil
Upper coil current [kA]
Lower coil current [kA]
Plasma coordinate [m]Plasma mass [mg]Time step ⊿t [nsec.]
Number of particles

0.25×0.536
1232
16(0,0,0)100
0.277
100000

Fig.A-5 矩形コイルシステムに対する 3D hybrid code での計算結果（左：初期状態，右：11.08μsec.後）
推進効率は 74 %，Steering angle β は 4.5 度を得た．
これらのことより矩形コイルシステムを用いることでレーザー核融合ロケットの推進性能を向上でき，整形ターゲットと併用することでより一層の性能向上が可能と考えられる．

59

また，Fig.A-5 の 11.08μsec.が経過した図を見ると，矩形コイルの間を粒子が通っている
様子がわかる．そこで，粒子が矩形コイルへ衝突する可能性を確かめた．
超伝導電磁石用線材として NbTi，安定化材として Al を用いる．安全係数を考慮した
NbTi の臨界磁場強度(Jc)を 1.6×108 A/m2 とするとコイルの断面積はコイルに流す電流が32 kA なので，Scsc = Ic / Jc ≈ 2×10-4 m2 = 2 cm2
1 つのコイルを半径 40cm の円と近似すると
Mc = 2RcScsc×(NbTi  0.5Al)
NbTi = 7×103 kg/m3
Al = 2.7×103 kg/m3
よってコイル重量 Mc ≈ 4.2 kg．
また，コイルは自身が作る磁場によってローレンツ力を受ける．その力を吸収するのがコイルカンでありそれをコイルの周りに配置する．
コイル

コイルカン
コイルの断面半径 ac は Scsc = 2 cm2 より
ac ≈ 8×10-3 m
よってコイルに働く力はFhoop =

1
8R c 3μ0 Ic 2 ln−
= 3.38 kN
2ac
4

これよりコイルカンの断面積 Scc，および重量 Mcc は，材質として用いるグラスファイバーの密度gf = 1.89×103 kg/m3，使用応力 σu = 330 MPa = 330×106 N/m2 より，
Scc = Fhoop / 2σu ≈ 1.63×10-6 = 1.63 mm2
Mcc = 2RcSccgf = 7.7 g
これらの結果よりコイルカンも含めたコイルの半径は約 1cm になる．

60

Table A まとめ
コイル重量

4.2kg

コイル断面半径0.8 cm

コイルカン重量

7.7 g

コイルカン断面半径0.1 cm－Z 方向

Fig.A-6 矩形コイル内の粒子の位置
黒枠はコイルを表している．図から矩形コイルのZ 方向側を粒子が通っている様子がわかる．これらの粒子のラーマー半径は平均で 1 mm 程度であるためコイルには衝突する可能性は無いと考えられる．

61

参考文献
(1) R. A. Hyde, L. Wood, 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)川渕 亮, 九州大学 工学部 エネルギー科学科 卒業論文(2006)(7) C. D. Orth: UCRL-LR-110500 (2003)(8) G. R. Liu , M. B. Liu: “ Smoothed Particle Hydrodynamics a meshfree particle
method ” World Scientific Publishing Co. Pte. Ltd., (2003)(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) M. Murakami, H. Nagatomo: Nuclear Instruments and Methods in Physics Research,
544 (2005), 67(14) S. Şahin, H. M. Şahin: Annals of Nuclear Energy, 28(2001), 1413
(15) Konstantin V. Vchivkov, et.al.: Jpn. J. Appl. Phys. 42 (2003), 6590
(16)Y. Kajimura, R. Kawabuchi, N. Matsuda and H. Nakashima: Transaction of fusion
technology, 51 No.2T (2007), 229

62