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

修 士 論 文論文名
2 次元放射流体コードを用いた磁場中プラズマの挙動解析

氏

名

指導教員名

富永 宙志

中島 秀紀 教授, 山本 直嗣 准教授

第 1 章 ................................................................................ - 1 序論 .................................................................................... - 1 1.1 研究背景 ................................................................... - 1 1.2 目的 ........................................................................... - 2 第 2 章 ................................................................................ - 3 レーザー推進 ...................................................................... - 3 2.1 磁気スラストチャンバーシステムとは ..................... - 3 2.2 レーザー推進のメリット .......................................... - 4 第 3 章 ................................................................................ - 6 計算コード ......................................................................... - 6 3.1 放射流体コード ......................................................... - 6 3.1.1 基礎方程式 ........................................................... - 7 3.1.2 時間ステップ........................................................ - 7 3.1.3 計算領域, 計算格子 ............................................. - 8 3.1.4 物理量の更新........................................................ - 9 3.1.5 熱伝導 ................................................................ - 10 3.1.6 レーザー吸収...................................................... - 10 3.1.7 輻射輸送 ............................................................. - 11 3.1.8 リーマン問題...................................................... - 12 3.2 3 次元ハイブリッドコード ..................................... - 16 3.2.1 基礎方程式 ......................................................... - 16 3.2.2 計算方法 ............................................................. - 18 3.2.3 プラズマ粒子の運動の計算 ................................ - 18 3.2.4 イオン数密度, 電流密度の計算 ......................... - 19 3.2.5 電磁場の計算...................................................... - 20 3.2.6 Digital Filter ..................................................... - 22 -

3.2.7 計算制限-CFL 条件 ............................................ - 23 3.3 2 次元カップリングコード ..................................... - 23 第 4 章 .............................................................................. - 28 計算概要 ........................................................................... - 28 4.1 体系 ......................................................................... - 28 4.2 イオンビーム電流の計算 ........................................ - 30 4.3 力積の計算 .............................................................. - 31 第 5 章 .............................................................................. - 32 結果と考察 ....................................................................... - 32 5.1 イオンビーム電流比較 ............................................ - 32 5.1.1 実験, 1D, 2D の比較 ....................................... - 32 5.1.2 磁場による比較 .................................................. - 37 5.2 力積比較 ................................................................. - 41 第 6 章 .............................................................................. - 42 結論と今後の展望............................................................. - 42 参考文献 ........................................................................... - 43 謝辞 .................................................................................. - 45 -

第1章

序論
1.1 研究背景
1961 年, 人類が初めて宇宙へと進出してから約 50 年. 現在も様々な形で宇宙空間の利用, 宇宙開発は進められ, 気象衛星, GPS 衛星など我々の生活に身近なものから, 国際宇宙ステーションでの宇宙空間, 無重力下における様々な実験,
火星・木星探査まで幅広くある. 日本でも 7 年間の宇宙航行を経て, 小惑星「イトカワ」にてサンプル採取を行い地球に帰還した日本の無人宇宙探査機「はやぶさ」を題材に映画化されて話題を呼んだ. また 2013 年 9 月 14 日, 惑星分光観測
衛星（SPRINT-A）「ひさき（HISAKI）」を搭載したイプシロンロケット試験機が内之浦宇宙空間観測所から打ち上げられた. 「ひさき」は約 2 ヶ月間の初期運用
後に実際に観測や運用を開始し, 今年 2014 年にはハッブル宇宙望遠鏡や X 線天文衛星「すざく」,ハワイの地上望遠鏡などと一緒に, それぞれの観測機器で木星を観測する予定である. 「ひさき」では, 木星のプラズマ環境のエネルギーがどのように供給されているかを調べ, また地球型惑星の大気が宇宙空間に逃げ出すメカニズムを調査することで初期の太陽系で何が起こっていたかを知ることを目指している[1].
これらの宇宙開発の先には有人惑星飛行, さらには他惑星に眠る資源, 他惑
星に居住できるようにするテラフォーミングという目的がある. 現在使われているロケット, 探査機に使用されている推進装置は化学推進と電気推進である.
化学推進は大きな推力を出すことができるが, すぐに燃料を使い果たしてしまう. 電気推進は燃料利用効率が高く, 長時間使用することができるが推力が小さい. そのため上記の目的遂行のためには, 半世紀の間使われてきたこれらの推進装置では非常に長い時間がかかってしまう. 有人宇宙探査において, 宇宙
飛行に時間がかかるということは, 宇宙飛行士の身体的, 精神的負担につながる. すなわち, 長時間の宇宙航行による宇宙線被ばくやカルシウム損失, 狭小
空間での長期生活による心理的な問題などである. この問題を解決するためには, できるだけ短い時間でミッションを終わらせる必要があり, そのために高速推進可能な宇宙推進機の開発が求められる[2].
-1-1.2 目的
高速推進可能なロケットとして, 熱核ロケットや, 核パルス推進, 磁気セイルなど様々なタイプの推進方法が提案, 研究されている. その中で, レーザー
推進という, レーザーによって生成されたプラズマを制御することで推進する技術が注目されている. またこの制御には, 磁気スラストチャンバーシステムという磁場とプラズマを相互作用させて推力を得る方法が考えられている. メ
リットとして, ①高比推力＆大推力, ②ミッションに応じた設計, ③長寿命, の3 点が挙げられる. このようなレーザー推進については 2 章で詳しく述べる.
現在, 本研究室ではレーザー推進について, レーザー生成プラズマを用いた実験が行われ, 磁気スラストチャンバーシステム内のイオンの挙動について研究している. また, 磁場中のプラズマに関する数値解析を, プラズマ生成過程を計算する「放射流体コード」と, 磁場中のイオン挙動を計算する「3Ｄハイブ
リッドコード」を使用して行い, 妥当性検証のため結果を相互に比較している.
近年の実験で主に行われているのは, イオンの広がりと推力の計測である.
イオンの広がりを見るため, チャージコレクタによるイオンビーム電流計測が行われているが, 実験と計算では電流のオーダーがずれている. その理由として, 数値解析では今まで 1 次元の放射流体コードを使用しており, プラズマのアブレーション方向を仮定していたことが挙げられる. そのため, イオンの広がり方をうまく再現できていなかった. またスラストスタンドを用いてコイル
磁場とプラズマの相互作用による推力を実験で計測することができたが, ノイズが大きく正確な値を出すことが難しかった. より推力を大きくするため, 磁場強度を高めようと努めている. コイルの抵抗や大電圧をかけることによるリ
スクなど簡単ではないが, これより先の実験にてその達成を目指している.
そこで本研究では, 2 次元の放射流体コードを用いて次元を上げることで, プ
ラズマアブレーションをより精度良く計算した. これは数値解析におけるイオ
ン挙動の再現性向上を目的としたものである. また将来的な実験で目指す強磁場における磁気スラストチャンバー内プラズマの数値解析を行った.
本論文の構成は次のとおりである. 第 2 章では磁気スラストチャンバーを利用したレーザー推進について述べる. 第 3 章ではプラズマ挙動の数値解析に用いたコードについて述べる. 第 4 章では数値解析の体系を述べる. 第 5 章では数値解析の結果, 実験結果との比較について述べ, 第 6 章では本研究についての結論と今後の展望を述べる.
-2-

第2章

レーザー推進2.1 磁気スラストチャンバーシステムとは本研究グループで扱うレーザー生成プラズマと磁場の相互作用を利用するプ
ラズマ推進器の研究は 1970 年代に Kantrowitz により提唱され, その目的は惑星間航行を目指した高推力な推進システムであり, その後さまざまな具体的提案や基礎研究がなされてきた[3]. まずはどのようにして推力を得るのか, その原理について示す.
Fig.2.1 は推進のメカニズムを表したものである.
(a) コイルに電流を流し磁場を生じさせ, ターゲットにレーザーを照射して等方膨張プラズマを生成する.
(b) プラズマは磁場の凍結によって, コイルの電流によって形成されたコイル
磁場を押し広げる. そしてその磁場を打ち消そうと, プラズマの反磁性電流によって反磁性磁場が形成される.
(c) コイルとプラズマのそれぞれに反磁性磁場とコイルの電流, コイル磁場と反磁性電流のローレンツ力が生じ, プラズマは宇宙船の推進方向とは逆方向の力を受け, 排出される. 同様にして, 推進器も推力を得る.

Fig.2.1

磁気スラストチャンバー推進メカニズム

つまりプラズマは押しのけた磁場から張力を受けた形になり, バネのようにして押し出される. その反作用で推進器も力を受けて進むのである.
-3-

2.2 レーザー推進のメリット
磁気スラストチャンバーシステムを利用したレーザー推進は上記のように推力を得るが, この推進方法が惑星間航行を含めこれからの宇宙開発に必要であると言える理由を示す. 以下にメリットである 3 点をまとめる.
① 電気推進機と比べて高いプラズマ排出速度を有し, 高比推力と大推力を同時に達成できる. この高い比推力は, 推進剤質量の低減によるロケットの積載可能質量の増加につながり, コスト削減につながる[4]. また大推力はミッションの短縮に直結する.
② レーザー生成プラズマはレーザーパワーと燃料推進剤の質量は独立関係なので, プラズマの排出質量と排出速度を広範囲に制御することができ, ミッシ
ョンに応じた最適な設計が行える[3,4].
③ プラズマは磁場によって排出方向を制御されるので, ノズルなどの固体壁に接触損傷を与えず, 長寿命な推進器となる.
ミッション期間の短縮は, 有人飛行において大きな障害となる宇宙線による被爆期間を短くすることができる. また, ミッションに応じた最適な設計が行えるため, 惑星間航行だけではなく地球近傍におけるミッションにおいても活躍することができると期待されている[3,4]. つまり, この推進器は地球近傍から惑星間航行までのあらゆるミッションに対応できる万能な推進器である.
また, この推進システムは構造が単純で, さらに燃料ペレットの爆発力をダ
イレクトに推進力に変換することが可能である.
このレーザー推進の研究としては, 米国 Lawrence Livermore 国立研究所の Orthらのグループが Fig.2.2 に示す概念設計例「VISTA」を提案した[7]. しかし, その推進機構に関して十分な議論がなされているとは言い難く, 実現可能性を実証するには至っていない.

-4-

Fig.2.2 レーザー推進の概念設計例「VISTA」
上図（http://www.interstellarindex.com/wp-content/uploads/2012/04/vista.png）
下図（http://thestarshipman.files.wordpress.com/2013/11/kflongplate14.jpg）

本研究グループは大阪大学レーザーエネルギー学研究センターの共同研究施設を利用して, 磁気スラストチャンバーについての模擬実験を行い, この推進機構の実験的実証を行っている. 実験以前よりこのシステムについてのシミュ
レーションが行われてきたが, 実験結果の再現性を高めるためシミュレーショ
ンコードを同研究センターと共同で改良した結果を 3 章以降で示していく.

-5-

第3章

計算コード
この章では, 次章以降の数値解析に使用する, 計算コードについて説明する.
最初にプラズマのレーザーによる生成過程について計算する放射流体コード,
次に磁場中に置いたプラズマの挙動について計算する 3D ハイブリッドコードについて説明する. 最後に, その後に使用する 3D ハイブリッドコードに, データを受け渡すカップリングコードについて説明する.

Fig.3.1

シミュレーションの流れ

3.1 放射流体コード
本研究で取り扱う放射流体コード[8]は, レーザー核融合のシミュレーションを行う爆縮コードとして, 大阪大学 砂原敦 氏が開発・提供しているものである.
「アブレーション, 衝撃波, レーザー光線追跡, レーザー吸収, 電子加熱, 熱
伝導, 状態方程式, 原子過程」などの物理モデルを取り扱っている. 本研究では,レーザーによるプラズマ発生の過程を計算するために使用させて頂いた. 今回
イオン挙動の再現性向上のため, 2 次元の放射流体コードを使用した.
注意する点として, 2 次元放射流体コードでは実際の 3 次元での現象を 2 次元で計算シミュレーションするので, 解析しようとする対象には物理量に対称性があることが必要条件となる. 今回のシミュレーションの場合, レーザー照射・ターゲット形状・ターゲット密度・磁場に同じ軸の対称性が仮定できるように条件設定した. 出力された 2 次元データの対称性を仮定した, 3 次元への拡張
は「3.3 2 次元カップリングコード」で説明する.
-6-

レーザーアブレーションを記述する流体近似のモデルとして次に示す 1 流体 2
温度の輻射流体力学モデル, または単に 2 温度モデルが広く用いられている.

3.1.1 基礎方程式
以下の方程式は直交座標系および, CGS 単位系で記述されている.
イオン, 電子の 2 温度を仮定した放射流体の方程式[9]を以下に示す.
𝑑𝜌= −𝜌𝛻 ∙ 𝒖𝑑𝑡

(3.1.1)𝑑𝒖= −𝛻(𝑝𝑖 + 𝑝𝑒 + 𝑞)𝑑𝑡

(3.1.2)𝑑𝑇𝑖= −𝑝𝑖 𝛻 ∙ 𝒖 − 𝛻 ∙ (𝜅𝑖 𝛻𝑇𝑖 ) + 𝛼(𝑇𝑒 − 𝑇𝑖 )𝑑𝑡
+ 𝑄𝑖

(3.1.3)𝑑𝑇𝑒= −𝑝𝑒 𝛻 ∙ 𝒖 − 𝛻 ∙ (𝜅𝑒 𝛻𝑇𝑒 ) − 𝛼(𝑇𝑖 − 𝑇𝑒 )𝑑𝑡
+ 𝑄𝑒

(3.1.4)

連続の式
運動方程式

𝜌

イオンエネルギー式

𝜌𝑐𝑣𝑖

電子エネルギー式

𝜌𝑐𝑣𝑒

ここで ρ は密度, u は流体速度, Ti はイオン温度, Te は電子温度, pi はイオン
圧力, pe は電子圧力, q は人工粘性, cvi はイオン比熱, cve は電子比熱, κi・κe はそれぞれイオン・電子の Spitzer の熱伝導係数エラー! 参照元が見つかりません。, α は Spitzerの熱緩和係数, Qi は“人口粘性に伴う運動エネルギーからイオンの熱エネルギーへの変換分”, Qe は“レーザー吸収とプラズマ中の放射吸収により電子が受け取るエネルギー”である. （※添え字 i, e はそれぞれイオン, 電子の物理量である
ことを示す. ）3.1.2 時間ステップ
時 間 ス テ ッ プ 幅 は 初 期 状 態 を 1 × 10−16 ,𝑠𝑒𝑐- と し て 開 始 , 最 小 値 を 1 ×
10−15 ,𝑠𝑒𝑐-, 最大値を5 × 10−11 ,𝑠𝑒𝑐-として逐一計算していく. イオンと電子の温度, 密度, 圧力が増加した割合を求めて, その中で割合が最大の値を選んでステップ幅を決める. 例えば最大値がイオンの温度の場合,
𝑑𝑡(3.1.5)𝑑𝑡′ =(𝑑𝑇𝑖 ⁄(𝑇𝑖 + 𝑇𝑖 𝑚𝑖𝑛 )/0.25
𝑑𝑇𝑖 は n-1 ステップから n ステップにかけてのイオン温度の増分である.

-7-

3.1.3 計算領域, 計算格子
2 次元放射流体コードの計算領域の概略を Fig.3.2 に示す.

Fig.3.2

2 次元放射流体コードの計算領域概略

イメージとして各物理量は z 軸まわりの対称性を持つことになる. r は z 軸を中心とした半径方向である.
計算領域は z 方向, r 方向をそれぞれ 1 から nz, 1 から nr 個の計算格子に分けている. 計算格子の座標値は格子境界(ib, jb)で定義され, 下付き添え字 b を用いて zb(ib), rb(jb)で表される. また格子中心(ic, jc)での座標値を(3.1.6), (3.1.7)式のように定義している.
(3.1.6)𝑧𝑐 (𝑖𝑐 ) = ,𝑧𝑏 (𝑖𝑏 ) + 𝑧𝑏 (𝑖𝑏 − 1)-/2(3.1.7)𝑟𝑐 (𝑗𝑐 ) = ,𝑟𝑏 (𝑗𝑏 ) + 𝑟𝑏 (𝑗𝑏 − 1)-/2
物理量は格子中心で定義され, 速度は vzc(ic), vrc(jc), 密度は ρc(ic, jc), 温度はTic(ic, jc), Tec(ic, jc), 圧力は pic(ic, jc), pec(ic, jc)で表される.
境界条件の設定を容易にするために, 上下左右の端に２つずつダミーのメッシュを置いている. ダミーメッシュには内側のメッシュで定義された物理デー
タがコピーされる. 反射境界条件の場合には境界の速度を 0 として内側の格子での速度の符号を反転させたものを設定する.
メッシュの位置については補外され,-8-

𝑧𝑏 (𝑛𝑥𝑚2) = 2𝑧𝑏 (𝑛𝑥𝑚2 + 1) − 𝑧𝑏 (𝑛𝑥𝑚2 + 2)
𝑧𝑏 (𝑛𝑥𝑚2 − 1) = 2𝑧𝑏 (𝑛𝑥𝑚2 + 1) − 𝑧𝑏 (𝑛𝑥𝑚2 + 3)
𝑧𝑏 (𝑛𝑥𝑝2 − 1) = 2𝑧𝑏 (𝑛𝑥𝑝2 − 2) − 𝑧𝑏 (𝑛𝑥𝑚2 − 3)𝑧𝑏 (𝑛𝑥𝑝2) = 2𝑧𝑏 (𝑛𝑥𝑝2 − 2) − 𝑧𝑏 (𝑛𝑥𝑝2 − 4)

(3.1.8)(3.1.9)(3.1.10)(3.1.11)

ここで nxm2 はマイナス側の境界部分のダミーメッシュ番号（一番外側）, nxp2
は逆にプラス側の境界部分のダミーメッシュ番号を指す. 上から順に, 境界を
ib=0（マイナス側）, ib=nx（プラス側）とすると ib=-1, ib=-2, ib=nx+1, ib=nx+2の位置である. r 方向についても同じように書ける.
Fig.3.3 に z 方向の計算格子のイメージを示す

z 方向の計算格子のイメージ

Fig.3.3

3.1.4 物理量の更新
運動方程式に従って物理量を更新していく.
速度の更新
𝑛+

1

𝑛−

1

𝑢𝑖𝑏 2 − 𝑢𝑖𝑏 2𝑛1
1 = 𝑆𝑖𝑏𝑛+𝑛−𝑡 2−𝑡 2𝑛∆𝑚𝑖𝑏=

(𝑝 + 𝑞)𝑛𝑖𝑐+1 − (𝑝 + 𝑞)𝑛𝑖𝑐
𝑛
∆𝑚𝑖𝑏

1
𝑛
𝑛
(∆𝑚𝑖𝑐
+ ∆𝑚𝑖𝑐+1)2(3.1.12)(3.1.13)1𝑛∆𝑚𝑖𝑏=0= (∆𝑚1𝑛 )2(3.1.14)1𝑛𝑛∆𝑚𝑖𝑏=𝑛𝑥= (∆𝑚𝑛𝑥)2(3.1.15)𝑛𝑆𝑖𝑏：メッシュの面積-9-密度の更新
𝑛
𝜌𝑖𝑐
=

∆𝑚𝑖𝑐
∆𝑉𝑖𝑐𝑛+1(3.1.16)

𝑉𝑖𝑐𝑛 ：メッシュの体積

3.1.5 熱伝導
熱伝導は古典的な Spitzer-Harm モデルで考えている. またある条件（長波長,低温）の時, レーザー吸収により加熱された電子は Maxwell 分布から外れることが知られている. 特に平均自由行程が長い場合, 電子による非局所的な熱伝導が顕著になり, 古典的な近似で現象を記述できなくなる. その場合はFokker-Planck 方程式による解析を行う予定であるがまだコードに組み込んでいない.

3.1.6 レーザー吸収
続いてレーザー吸収について記す. レーザーは計算領域端からプラスチック
ターゲットに垂直に入射される. プラスチックは金属のように導電体ではなく誘電体であるために自由電子は多光子吸収やトンネル電離などで生じる. ひと
たび自由電子が生成されると, レーザーの電場で揺り動かされた（加熱された）自由電子が中性粒子またはイオン粒子と三体衝突することにより振動エネルギ
ーを熱化してレーザーの光子エネルギーを吸収して電子が加速する. この過程
を制動放射（電子がイオン近傍で減速した時にエネルギーを放出する現象）の逆過程であるため, 名の通り逆制動放射（Inverse-bremsstrahlung）という. レー
ザーエネルギー吸収は逆制動放射を仮定し, レーザー光線に沿って, 局所的な電子密度 ne, 平均電離度 Z*, 電子温度 Te から求められる吸収係数を用いて各位置におけるレーザーエネルギーの吸収を計算した.
3

1 1 𝑚3⁄2 4𝜋𝑒 2
𝜅= 2() 𝑛𝑒2 𝑧 ∗ 𝑙𝑛 𝛬 (𝑘𝑇𝑒 )−2⁄3
𝜔 3 2𝜋 3⁄2 𝑚
ここで,
𝜅：単位長さあたりの吸収係数 [1/cm]
ln Λ：温度に対するクーロン対数
𝑚𝑣2

𝑓̃(𝑣) = e− 2𝑘𝑇である.

- 10 -

1

̃
2 𝑓0 (𝑣 = 0)𝜔𝑝𝑒𝑐 (1 − 2 )𝜔(3.1.17)

3.1.7 輻射輸送
次に輻射輸送について記す. 輻射スペクトルは量子力学的な準位間の遷移過程によって決まるので, 電子の束縛状態間の遷移は, 線スペクトルとなってエ
ネルギー空間における変動が激しい. ここで特定の光子エネルギーへのすべての輻射遷移過程を逐次評価する手法を line-by-line 法と呼ぶ. この手法は輻射ス
ペクトルの詳細を知りたい時には有効であるが, 輻射流体計算には膨大な時間を要するため不向きである. そこでエネルギー空間をある幅を持った領域で分割し, その領域ごとに輻射輸送方程式を解くという手法が取られる. これを多群近似, またはマルチバンドモデルと呼ぶ[10]. また, 輻射輸送は一般的に空間
3 次元, 運動量空間 3 次元の 6 次元の位相空間でボルツマン方程式を解くことになるが, 計算負荷が大きく現実的ではないため近似を用いる[11]. 放射流体コードにおいても熱束制限多群拡散近似エラー! 参照元が見つかりません。により光子のエネルギーを多数のエネルギー群に分け, それぞれのエネルギーについて拡散近似 エラー!
参照元が見つかりません。
により反復して解いて輻射エネルギー密度を求める. ２次元計算では計算負荷が大きいため光子のエネルギー群は 40 グループに分割し, EUV 領
域では 0.4 [eV]刻みである. 輻射輸送の 0 次と 1 次のモーメント式を(3.1.18),(3.1.19)式に示す.
0次𝜕𝐸𝜈+ 𝛻 ∙ 𝑭𝜈 = 4𝜋𝜂𝜈 − 𝑐𝜒𝜈 𝐸𝜈
𝜕𝑡(3.1.18)

𝜕𝑭𝜈
(3.1.19)+ 𝛻 ∙ 𝑷𝜈 = −𝑐(𝜒𝜈 + 𝜍𝑠 )𝐸𝜈𝜕𝑡
ここでEν は周波数 ν に対するエネルギー密度, 𝑭𝜈 は輻射エネルギー流速, 𝑷𝜈
は輻射テンソル, Dν は拡散係数, χν は吸収率, 𝜂𝜈 は放射率, c は光速である. 輻
1次

1

射場を等方的とすると, 𝑷𝜈 = 3 ∇𝐸𝜈 である. また(3.1.19)式左辺第１項を無視する. (3.1.20)式に最も簡単に近似した拡散型の輻射輸送方程式を示す.
𝜕 𝐸𝜈
(3.1.20)
𝜌 ( ) − 𝛻 ∙ (𝐷𝜈 𝛻𝐸𝜈 ) = 4𝜋𝜂𝜈 − 𝑐𝜒𝜈 𝐸𝜈
𝜕𝑡 𝜌(3.1.20)式の右辺は流体からの放射によるエネルギー損失を表しており, 積分
した値は電子のエネルギー式の生成項𝑄𝑒 に一部として代入される. また(3.1.18),(3.1.19)式において輻射エネルギー流速は(3.1.21)𝑭𝜈 = −𝐷𝜈 𝛻𝐸𝜈
としている. ここで拡散係数は,
𝐷𝜈 =

𝑐
3𝜒𝜈

(3.1.22)

であるが, 近似において𝜕𝐹𝜈 ⁄𝜕𝑡 = 0としているため, ∇Eν が大きくなると輻射エ
ネルギー流速は追随してどこまでも大きくなり得る. そこで, 𝑭𝜈 = −𝑐𝐸𝜈 を超え- 11 -

ると非物理的で解が正確ではなくなるため, 拡散係数を(3.1.22)のように書き換
える.
𝐷𝜈 =

𝑐
|𝛻𝐸 |
3𝜒𝜈 + 𝐸 𝜈𝜈

(3.1.23)

これにより輻射場が等方的極限|∇𝐸𝜈 |⁄𝐸𝜈 ≪ 1 と非等方的極限で物理的に正しい解を得ることができる. このモデルを熱束制限拡散近似と言う.
(3.1.1)式~(3.1.23)式において𝑝𝑖 , 𝑝𝑒 , 𝑐𝑣𝑖 , 𝑐𝑣𝑒 , 𝜂𝜈 , 𝜒𝜈 は密度 ρ と温度 T の関数として表される. 本研究では, これらの値は米国の Los Alamos 国立研究所で作成された”SESAME”と呼ばれるデータテーブルを引用し, 状態方程式を解くことで求めている[12].

3.1.8 リーマン問題
衝撃波管問題とは, 一定の断面積をもった十分長い管の中に比熱比が等しく
熱力学的状態の異なる流体を 1 枚の仕切りによって分けて入れておき, その仕切りを取り去った後の流体の状態を求める問題である[13].
隔膜（仕切り）が瞬間的に破られると, 衝撃波（Ｓ）が発生する. 高圧部から低圧部へ移動. 入射衝撃波と呼ばれる急勾配の衝撃波面が形成される. この衝撃波は, 低圧の流体の温度と圧力を急激に上昇させ, 衝撃波と同じ方向に向かう流体の流れを作る. また同時に膨張波（Ｅ）が高圧側へ進行. 低圧と高圧の流体の断面は接触面（Ｃ）と呼ばれ, 衝撃波を追って急速に移動する. このような現象の流体の状態を求める問題をリーマン問題と呼ぶ.

Fig.3.4

リーマン問題の概略図- 12 -

この解には, 衝撃波, 膨張波, 接触面（接触不連続）という典型的な物理現象が含まれ, 厳密解も求められているため, 流体計算コードのテストにもよく用いられる.
リーマン問題を処理するため流体の式を解いていくが, その方法に関しては様々な解法が考えられている. その中で今回採用したのは, HLLC 法である. これはリーマン解法を近似的に解くため, 近似リーマン解法と呼ばれている中の一つである. HLLC 法について記述する前に尐々HLL 法について触れる.
HLL 法は Fig.3.4 において最大と最小の特性速度（両端に広がる波の速度）だ
けを考え, その内側の領域は中間状態として構造を無視して一様と仮定する.
そうすることで２つの特性速度が分かれば式を解くことができ, 導出を非常に簡単にした解法である. Harten, Lax, van Leer らによって提案された解法であり, 彼らの頭文字からとって HLL 法という名前がついている. ただ, この解法には問題点があり, それは中間の領域すべてを一様としているため, 接触不連
続面を正確に捕えることができない. そこで, その接触面の特性速度を考え適度に近似したのが HLLC 法である. この C は「Contact Surface」の C である. 以
下で簡単に解法について記述する[14,15,16]. 簡単のために一次式を考える.
(a)概略図- 13 -

(b)計算格子

HLLC 法

Fig.3.5

まず流体の方程式を保存系で次のように書く.
𝜕𝜕
(3.1.24)
𝑈+
𝐹=0
𝜕𝑡
𝜕𝑥
ここで𝑈は保存量, 𝐹は流束である. (3.1.23)式を差分化して次のような形に書き換える. Fig.3.5 (b)に示しているように, ステップ𝑛の値を使って次のステップの保存量を求める.
𝑈𝑖𝑛+1 = 𝑈𝑖𝑛 −

∆𝑡 𝑛
𝑛
(𝐹
− 𝐹𝑖−1⁄2 )∆𝑥 𝑖+1⁄2(3.1.25)

この式の中で既知であるのは 𝑈𝑖𝑛 であり, 𝐹 𝑛 を求めればステップ 𝑛 + 1の値が求まる.
(3.1.24)式を積分して値を代入して整理していくと, 下式のようになる.
(3.1.26)𝑆𝐿 𝑈𝐿∗ − 𝐹𝐿∗ = 𝑆𝐿 𝑈𝐿 − 𝐹𝐿∗∗
(3.1.27)
𝑆𝑅 𝑈𝑅 − 𝐹𝑅 = 𝑆𝑅 𝑈𝑅 − 𝐹𝑅
ここで𝑆𝐿 は膨張波の特性速度, 𝑆𝑅 は衝撃波の特性速度, 𝑈𝐿∗ , 𝑈𝑅∗ は中間状態の保存量,
𝐹𝐿∗ , 𝐹𝑅∗ は中間状態の流束である. 例えば運動量を考えて, 𝑈 = 𝜌𝑢, 𝐹 = 𝜌𝑢2 + 𝑃として(3.1.26), (3.1.27)式に代入すると,
𝜌𝐿∗ =𝑆𝐿 − 𝑢𝐿𝑆𝐿 − 𝑆𝑀(3.1.28)𝜌𝑅∗ =𝑆𝑅 − 𝑢𝑅
𝑆𝑅 − 𝑆𝑀(3.1.29)

𝑆𝑀 は接触面の特性速度. このように中間状態の密度を各特性速度と流体の速度- 14 -

という既知の値から導くことができた. 𝑢𝐿∗ = 𝑢𝑅∗ = 𝑆𝑀 としているため, これで時間ステップ𝑛の値がすべて分かり, 次ステップの保存量を求めることができる. 上
記の式を見ても分かる通り, HLLC 法は非常にシンプルで計算コストも小さい上,厳密解にも近い値を出せることがよく知られている. その例を Fig.3.6 に示す.
(a)

(b)

Fig.3.6

HLLC 法の使用例 (a):[14],(b):[17]

上記に示すのは HLLC 法と厳密解の比較である. (a)は細いラインが厳密解で緑と赤の太いラインが HLL 法と HLLC 法による近似解である. 接触不連続面以外に関しては膨張波と衝撃波とよく一致しているのが分かる. また(b)は HLL,
HLLC 法のほかに Linde や Roe といったスキームとも比較している. 右の図において HLLC 法が HLL 法よりも厳密解に近いことが分かる. また, 計算時間はHLL と HLLC 法が他のスキームよりも短いことが報告されている.

- 15 -

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

3.2.1 基礎方程式計算に使用する方程式を導出するための基礎方程式は次のように与えられる.
式は直交座標系および, MKSA 単位系で記述されている.
イオンの運動方程式

𝑚𝑖

電子の運動方程式𝑚𝑖𝑑𝒗𝑖= 𝑍𝑒(𝑬 + 𝒗𝑖 × 𝑩)
𝑑𝑡𝑑𝒗𝒊
= 𝑍𝑒(𝑬 + 𝒗𝑖 × 𝑩)𝑑𝑡
𝛻 × 𝑩𝑝 = 𝜇0 (𝑱𝑒 + 𝑱𝑖 )アンペアの法則

𝜕𝑩= −𝛻 × 𝑬𝜕𝑡電子の状態方程式
𝑃𝑒 = 𝑛𝑒 𝑇𝑒
以後計算に使用する方程式の導出を説明する.
まず粒子として取り扱うイオンの運動方程式を示す.
ファラデーの電磁誘導の法則(3.2.1)(3.2.2)(3.2.3)(3.2.4)(3.2.5)𝑑𝒗𝑖= 𝑍𝑒(𝑬 + 𝒗𝑖 × 𝑩)
𝑑𝑡

(3.2.6)𝑑𝒙𝑖 𝑑𝒚𝑖 𝑑𝒛𝑖
𝒗𝑖 = (,,)
𝑑𝑡 𝑑𝑡 𝑑𝑡

(3.2.7)

𝑚𝑖

𝑑𝒙𝑖
𝑑𝒚𝑖
𝑑𝒛𝑖
(3.2.8), 𝒗𝑦𝑖 =
, 𝒗𝑧𝑖 =
𝑑𝑡
𝑑𝑡
𝑑𝑡
ここで vi はイオンの速度, vxi, vyi, vzi はそれぞれ x, y, z 方向のイオン速度, Zはイオンの価数, e は電気素量, xi, yi, zi はイオンの位置, mi はイオンの質量,
また E , B はそれぞれ電場, 磁場を表している.
𝒗𝑥𝑖 =- 16 -

流体として扱う電子の運動方程式を電場の解法に用いる. 電子の運動方程式は以下の式で表される.
𝑑𝒗𝑒(3.2.9)= −𝑒𝑛𝑒 (𝑬 + 𝒗𝑒 × 𝑩) − 𝛻𝑃𝑒𝑑𝑡
ここで ve は電子の速度, me は電子の質量, ne は電子数密度, Pe は電子圧力である.
電子は慣性無しの流体としてみなすので, 慣性による効果は無視し, (3.2.9)式の左辺は 0 にできる(me→0).
さらに電子圧力 Pe に関しては電子の状態方程式より以下のようになる.
(3.2.10)
Pe = ne Te
𝑛𝑒 𝑚𝑒

ここで Te は電子温度である. (3.2.10)式を(3.2.9)式に代入して Pe を消去する.
(3.2.11)0 = −𝑒𝑛𝑒 (𝑬 + 𝒗𝑒 × 𝑩) − 𝛻(𝑛𝑒 𝑇𝑒 )荷電準中性条件および, ダーウィン近似を用いたアンペアの法則により以下の式が得られる.
𝛻 × 𝑩𝑝 = 𝜇0 (𝑱𝑒 + 𝑱𝑖 )(3.2.12)
ここで μ0 は真空中の透磁率であり, Ji , Je はそれぞれイオン電流密度, 電子電流密度である. また, Bp はプラズマにより生成される磁場で以下のようになる.
𝑩𝑝 = 𝑩 − 𝑩𝑒𝑥𝑡
(3.2.13)Ji はイオンの数密度と電流密度から計算することが出来る. Je は電子電流密度の定義により以下のようになる.
(3.2.14)
𝑱𝑒 = −𝑒𝑛𝑒 𝒗𝑒(3.2.12),(3.2.14)式を用いて Je を消去し ve について求めると以下のようになる.
𝒗𝑒 = −

1 1( 𝛻 × 𝑩𝑝 − 𝑱𝑖 )𝑒𝑛𝑒 𝜇0(3.2.15)

これを(3.2.11)式に代入して電場 E について整理する. 荷電準中性条件 ne＝niより次式を得る.
1 1
𝑬={ (𝛻 × 𝑩𝑝 ) × 𝑩 − 𝑱𝑖 × 𝑩 − 𝑍𝑇𝑒 𝛻𝑛𝑖 }(3.2.16)
𝑍𝑒𝑛𝑖 𝜇0
但し, (3.2.16)式から得られる電場はプラズマが存在する領域内でのみ使用さ
れ, プラズマが存在しない真空領域での電場は後節で述べるラプラス方程式から得られる.
プラズマによる磁場を記述する方程式として, 以下のファラデーの電磁誘導の法則を用いる.
𝜕𝑩= −𝛻 × 𝑬
𝜕𝑡- 17 -(3.2.17)

3.2.2 計算方法
(3.2.6)式から(3.2.16)式を数値計算によって解くためには, 基礎式を差分式に直し, 数値的に安定で, 物理的に妥当な空間幅, 時間幅で解かなければならな
い.
計算の手順を以下に示す.
[1] 外部磁場および, イオンの粒子位置, 速度の初期設定を行う.
[2] 格子点上でのイオンの数密度, 電流密度を計算する.
[3] (3.2.15)式により, 電場を求める.
[4] [1]で設定した外部磁場と, [3]で求めたプラズマによる電場と磁場から,(3.1.6)式, (3.2.7)式を用いて新たなイオン粒子の位置, 速度を求める.
[5] (3.2.16)式により磁場を求める.
[2]～[5]を繰り返し, プラズマの時間的変化を追う.
3 次元ハイブリッドコードに対する各方程式の差分式を以下の小節で示す.

3.2.3 プラズマ粒子の運動の計算
3 次元ハイブリッドコードでは計算負荷を減らすために, いくつかのイオン粒子を１つの大きな電荷, 大きな質量をもった粒子と見なす. このような粒子を超粒子と呼ぶ. 超粒子を個別に運動させることでプラズマ粒子の運動を追う.
今回の計算では超粒子数を 100 万個としてシミュレーションを行った. (3.2.7)式の x 方向成分に, leap-frog 法[20]による時間中心差分を施すと次式となる.
𝑥𝑖𝑛+1 − 𝑥𝑖𝑛
𝑛+1/2(3.2.1)= 𝑣𝑖∆𝑡
ここで各変数の上付き添え字は時間を表し, 下付き添え字は i 番目のイオンを表す. また, Δt は時間幅を表す.
(3.2.6)式に時間中心差分を施すと以下の式となる.
𝑛+1/2

𝒗𝑖

𝑛−1/2

− 𝒗𝑖
𝑍𝑒 𝑛 1 𝑛+1/2
𝑛−1/2(3.2.2)={𝑬 + (𝒗𝑖+ 𝒗𝑖) × 𝑩𝑛 }∆𝑡𝑚𝑖2
n −1/ 2 はある時刻 n から Δt / 2 だけ以前の値であることを示す. 同様に, n +1/
2 は Δt / 2 だけ進んだ値であることを示す. 上付き添え字が n , n −1/ 2 である変数は既知量である.
n+1/2

プラズマ粒子の運動を求めるためには, まず(3.2.1)式において, 𝐯i
を左辺
n+1/2
n+1/2
にまとめて, 未知量である𝐯iを解く. 次に, 求めた𝐯iを(3. 20)式の右辺
n+1 n+1 n+1に代入し, 新たな粒子の位置(xi , yi , zi ) を求めることが可能となる.
(3.2.1)式を変形すると次式となる.
(3.2.3)
𝑨𝒗𝑛+1/2 = 𝑺- 18 -

ここで行列 A , ベクトル S , α , β , γ , δ を次のように定義している.
1
𝑨=( 𝛾
−𝛽

𝛾
1
𝛼

𝛽
−𝛼 )
1

(3.2.4)𝑚𝑣𝑖𝑥+ 𝑐𝐸𝑥𝑛 + 𝛾𝑣𝑦𝑚 − 𝛽𝑣𝑧𝑚𝑚
𝑺 = (𝑣𝑖𝑦 + 𝑐𝐸𝑦𝑛 + 𝛾𝑣𝑧𝑚 − 𝛽𝑣𝑥𝑚 )𝑚𝑣𝑖𝑧+ 𝑐𝐸𝑧𝑛 + 𝛾𝑣𝑥𝑚 − 𝛽𝑣𝑦𝑚

(3.2.5)

𝛼=

𝑞∆𝑡 𝑛
𝐵
2𝑚𝑖 𝑥

(3.2.6)

𝛽=

𝑞∆𝑡 𝑛
𝐵
2𝑚𝑖 𝑦

(3.2.7)

𝛾=

𝑞∆𝑡 𝑛
𝐵
2𝑚𝑖 𝑧

(3.2.8)

𝛿 = 1 + 𝛼 2 + 𝛽2 + 𝛾 2(3.2.9)なお, c = q Δt /mi , m = n −1/ 2 とおいている. また E , B は粒子位置におけるそれぞれの値を示している.
さらに, (3.2.3)式を次のように変形して, 速度を求める式が得られる.
(3.2.10)𝒗𝑛+1/2 = 𝑨−1 𝑺
ここで, A−1 は以下のように得られる.
1 + 𝛼21
𝑨−1 = (𝛼𝛽 − 𝛾𝛿𝛼𝛾 + 𝛽𝛼𝛽 + 𝛾1 + 𝛽2𝛽𝛾 − 𝛼𝛼𝛾 − 𝛽𝛽𝛾 + 𝛼 )1 + 𝛾2(3.2.11)3.2.4 イオン数密度, 電流密度の計算
イオン数密度を求めるためには, 空間内に広がっているプラズマを格子点上に割り付けなければならない. そのために PIC（Particle In Cell）法を用いる.
今, １つの超粒子の中に N 個のイオンが含まれている場合を考える. Fig.3.7に示すように, 格子点内に存在する超粒子の位置によって, 8 個の格子点によって囲まれた体積 V の立方体を, V1 , V2 , ～ V8 の体積を持つ直方体に分割す
る. 格子点(i,j,k)に割り振られる超粒子に含まれるイオンの数は次式で与えられる.
𝑉2(3.2.1)𝑁𝑉
すべての粒子についての和を δni についてとれば, ある格子点上の ni を定義することが出来る.
イオン電流密度 Ji は, 数密度 ni を用いて次式により求めることが出来る.
𝛿𝑛𝑖 (𝑖, 𝑗, 𝑘) =- 19 -

𝑱𝑖 = 𝑍𝑒𝑛𝑖 𝒗𝑖

Fig.3.7(3.2.2)

PIC(Particle In Cell)法

3.2.5 電磁場の計算(1) プラズマ領域での計算
電磁場を求めるためには, (3.2.15)式, (3.2.16)式を用いる. (3.2.16)式を片側差分して, 𝐵 𝑛+1/2 だけを左辺に残すと次式となる.
𝑩𝑛+1/2 = 𝑩𝑛 −∆𝑡
𝛻 × 𝑬𝑛2(3.2.1)

また, (3.2.15)式は以下の式になる.
1 1
1 𝑛
𝑇𝑒
𝑬𝑛 = {(𝛻 × 𝑩𝑛𝑝 ) × 𝑩𝑛 −𝑱𝑖 × 𝑩𝑛 − 𝛻𝑛𝑖𝑛 }𝑛𝑖 𝜇0 𝑍𝑒𝑍𝑒𝑒
ここで, 𝐁 n , 𝐄n は既知であるとすると,(3.2.2)(3.2.31)式から𝐁 n+1/2 を求めるこ
n+1/2とが出来る. 次に, 𝐁 n+1/2 , 𝐉 n+1/2 , 𝐯 n+1/2 , niを(3.2.32)式に代入し
n∗1
n+1
n+1/2
て 𝐄
を求め, 線形外挿より電磁場の予測値 𝐄pred , 𝐁pred
を次式より
求める.
𝑛
𝑛+1/2
𝑬𝑛+1
𝑝𝑟𝑒𝑑 = −𝑬 + 2𝑬

(3.2.3)∆𝑡(3.2.4)
𝛻 × 𝑬𝑛+1
𝑝𝑟𝑒𝑑
2
n+1
n+1
ここで, 下付き添え字 pred は予測値であることを表している. 𝐄pred, 𝐁pred
𝑛+1/2
𝑩𝑛+1
−
𝑝𝑟𝑒𝑑 = 𝑩

n+3/2

n+3/2

n+3/2

を求めた後, それらによるイオンの運動を解き, 𝐯i pred , ni pred , 𝐉i pred を- 20 -

n+3/2求める. 次に, 𝐁i pred を次式により求める.
𝑛+3/2

𝑩𝑝𝑟𝑒𝑑 = 𝑩𝑛+1
𝑝𝑟𝑒𝑑 −∆𝑡𝛻 × 𝑬𝑛+1/2
2
n+3/2

そして, (3.2.32)式で 𝐄n を求めた時と同様にして, 𝐄predを求める. 最後に𝐄n+1 , 𝐁 n+1を求める.

次式に示す線形内挿から次のタイムステップの値,
𝑬𝑛+1 =(3.2.5)

1 𝑛+1/2 1 𝑛+3/2
𝑬+ 𝑬𝑝𝑟𝑒𝑑22(3.2.6)∆𝑡
𝛻 × 𝑬𝑛+12(3.2.7)

𝑩𝑛+1 = 𝑩𝑛+1/2 −
n+3/2

n+3/2

n+3/2

これらの値𝐄n+1 , 𝐁 n+1 より, 𝐯i
, ni
, 𝐉i
を求めることが出来る.
ここで, 電場を求める際に用いた一連の手法は, 予測子－修正子法と呼ばれる収束計算の一種である. この方法では, n タイムステップから n +1 タイムステ
ップの値を求めるために, n + 3/ 2 タイムステップの予測値を求め内挿計算を行う. そして, 求めた n +1 タイムステップの値を再度予測値として用いる. 具体
n+1
的には, (3.2.37)式により求められた𝐁 n+1 を次の予測値𝐁predとし, (3.2.33)式から(3.2.37)式までを次式で定義する条件[21]を満足するまで, 繰り返し計算を行う.
𝑛+1(𝑞+1)

𝑚𝑎𝑥𝑖,𝑗,𝑘 |𝑩𝑖,𝑗,𝑘𝑛+1(𝑞)− 𝑩𝑖,𝑗,𝑘𝑛+1(𝑞)

𝑚𝑎𝑥𝑖,𝑗,𝑘 |𝑩𝑖,𝑗,𝑘

|

|

≤𝜖

(3.2.8)

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

(2) 真空領域での電場の計算
真空領域では, ラプラス方程式∇2 𝐄 = 0 を SOR（Successive Over Relaxation）
法を用いて解く. SOR における加速係数 α は正方形領域の場合は次のように表
される.
2
𝛼=
(3.2.9)𝜋
1 + 𝑠𝑖𝑛( ⁄𝑁𝐺𝑀𝑋)- 21 -

ここで NGMX は x, y または z 方向の格子数である. α の値は通常 1.8 から
1.9 付近である. また, α= 1 の場合には, SOR 法はガウス－ザイデル法に等しい.
SOR 法では, 次式によって変数の収束を加速する.
𝑬𝑝 (𝑖, 𝑗, 𝑘)
𝑬𝑝 (𝑖 − 1, 𝑗, 𝑘) + 𝑬𝑝 (𝑖 + 1, 𝑗, 𝑘) + 𝑬𝑝 (𝑖, 𝑗 − 1, 𝑘)= 𝛼{}⁄6+𝑬𝑝 (𝑖, 𝑗 + 1, 𝑘) + 𝑬𝑝 (𝑖, 𝑗, 𝑘 − 1) + 𝑬𝑝 (𝑖, 𝑗, 𝑘 + 1)
(3.2.10)+ (1 − 𝛼)𝑬𝑝 (𝑖, 𝑗, 𝑘)(𝑝 = 𝑥, 𝑦, 𝑧)
次式で定義される残差𝜉𝑥,𝑦,𝑧 が, あらかじめ決められた値以下になれば真空部
分の電場の計算を終了する.
𝑬𝑝 (𝑖 − 1, 𝑗, 𝑘) + 𝑬𝑝 (𝑖 + 1, 𝑗, 𝑘) + 𝑬𝑝 (𝑖, 𝑗 − 1, 𝑘)+𝑬𝑝 (𝑖, 𝑗 + 1, 𝑘) + 𝑬𝑝 (𝑖, 𝑗, 𝑘 − 1)
𝜉𝑝 = {}(𝑖,(𝑖,+𝑬𝑝 𝑗, 𝑘 + 1) − 6𝑬𝑝 𝑗, 𝑘)(3.2.11)(𝑝 = 𝑥, 𝑦, 𝑧)

3.2.6 Digital Filter
PIC 法で格子点上の量を決める場合, 例えば粒子の位置を格子上に割り振る
場合には, メッシュ内に存在する粒子の数が尐ない領域ではその付近の格子に比べて数密度が小さくなり, 数密度の変化が不連続になる可能性がある. このような現象を避けるために Digital Filter と呼ばれる一種の補間を施すことがある. Digital Filter は, 格子上に定義される値, 例えば電場, 数密度等を計算する際に発生する高周波の雑音を除去する目的にも用いられる[22]. Digital Filter の式は以下のように表される.
𝑠𝑖𝑑𝑒
𝑐𝑜𝑟𝑛𝑒𝑟𝑀𝛷𝑖,𝑗,𝑘 + 𝑆𝛷𝑖,𝑗,𝑘+ 𝐾𝛷𝑖,𝑗,𝑘(3.2.1)
𝑀 + 6𝑆 + 12𝐾
new
Φi,j,k を格子点上に定義される値とし, Φi,j,k
は Digital Filter を通した新しい𝑛𝑒𝑤𝛷𝑖,𝑗,𝑘←side
cornerΦi,j,k を表す. またΦi,j,kとΦi,j,kは side term と corner term であり,それぞれ以下の式で表される.
𝑠𝑖𝑑𝑒
𝛷𝑖,𝑗,𝑘= 𝛷𝑖−1,𝑗,𝑘 + 𝛷𝑖+1,𝑗,𝑘 + 𝛷𝑖,𝑗−1,𝑘 + 𝛷𝑖,𝑗+1,𝑘 + 𝛷𝑖,𝑗,𝑘−1+ 𝛷𝑖,𝑗,𝑘+1
𝑐𝑜𝑟𝑛𝑒𝑟𝛷𝑖,𝑗,𝑘 = 𝛷𝑖,𝑗−1,𝑘−1 + 𝛷𝑖,𝑗+1,𝑘−1 + 𝛷𝑖,𝑗+1,𝑘+1 + 𝛷𝑖,𝑗−1,𝑘+1+ 𝛷𝑖−1,𝑗,𝑘−1 + 𝛷𝑖−1,𝑗,𝑘+1 + 𝛷𝑖+1,𝑗,𝑘+1+ 𝛷𝑖+1,𝑗,𝑘−1 + 𝛷𝑖−1,𝑗−1,𝑘 + 𝛷𝑖−1,𝑗+1,𝑘+ 𝛷𝑖+1,𝑗+1,𝑘 + 𝛷𝑖+1,𝑗−1,𝑘- 22 -(3.2.2)(3.2.3)

ここで, M, S, K の値をそれぞれ 12, 4, 1 とすると(3.2.1)式は次式となる.
1
1 𝑠𝑖𝑑𝑒1 𝑐𝑜𝑟𝑛𝑒𝑟(3.2.4)
𝛷𝑖,𝑗,𝑘 + 𝛷𝑖,𝑗,𝑘+ 𝛷𝑖,𝑗,𝑘4
1248
数密度, イオン電流密度, 電場を計算した後に, それぞれの値にこの Digital
Filter を作用させる.
𝑛𝑒𝑤
𝛷𝑖,𝑗,𝑘
←3.2.7 計算制限-CFL 条件
ハイブリッドコードにおいて, 差分化の時間幅 Δt , 空間幅 Δx は次の条件を満たす必要がある.
CFL（Courant-Friedrichs-Lewy）条件∆𝑥∆𝑡 <(3.2.1)
𝑽𝐴
これは, Δt, Δx →0 で解が収束する条件である. 𝐕A はアルフベン速度であり次式で定義する.
𝑩0
𝑽𝐴 ≡(3.2.2)𝜇0 𝑛𝑖 𝑚𝑖また, 付加的な数値的安定条件は以下である.
(3.2.3)𝜔𝑐𝑖 ∆𝑡 < 0.2
ここで, ωci はイオンサイクロトロン周波数である.
メッシュ幅の最小値を求める条件として以下の式がある.
∆𝑥 > 𝑐 ⁄𝜔𝑝𝑖
(3.2.4)このc⁄ωpi より∆x が小さい場合にもハイブリッドコードは発散することなく計算を行えるが, その解が正しいものである保証は無い. ここで, ωpi はイオ
ンプラズマ周波数である.

3.3 2 次元カップリングコード
このコードは放射流体コードのデータを 3 次元ハイブリッドコードで扱えるようにするためにデータの変換をし, 2 つのコード間を“連結”する役割を持つ.
具体的には単位系, 座標系, 次元数, またイオンの取り扱い方などの違いに対するデータ変換を行う.
本研究において, カップリングで 3 次元ハイブリッドコードに受け渡しが必
要なデータは「イオン位置座標とイオン速度, イオン電荷, イオン質量」である.
ここではイオンに関連する物理量のみで電子に関する物理量が含まれない.
(3.2.15)式からわかるように 3 次元ハイブリッドの計算で必要な電子の物理量は電子温度 Te のみである. しかし過去の計算結果から電子温度 Te は 3 次元ハイブリッドコードの計算結果にあまり影響を与えないことがわかっているので,- 23 -

計算時間短縮のため Te=0 とした. そのため電子に関する物理量の受け渡しは必
要ない.
これらの必要なデータは 2 次元放射流体コードで求められた各計算格子の境界で定義される位置座標と, 各計算格子の中心で定義される速度, 質量密度,
イオン数密度, イオン電荷を変換して得る.

行ったデータ変換の流れは以下のようになる.
① 2 次元放射流体コードデータを CGS 単位系から MKS 単位系に変換する.
② 流体であったイオンを粒子に置き換えて、イオン粒子を“超粒子”として扱う.
③ 座標系を変換する.
以下にその詳細を示す.
①2 次元放射流体コードのデータの CGS 単位系から MKS 単位系への変換
z [cm]→

z×10-2 [m]r [cm]→

r×10-2 [m]vz [cm/s]→

vz×10-2 [m/s]vr [cm/s]→

vr×10-2 [m/s]質量密度p [g/cm3]→

ρ×103 [kg/m3]

イオン数密度ni [cm-3]→

ni×106 [m-3]

位置座標

速度

表1 ：Coupling code --単位系の変換②イオン粒子の“超粒子”としての取り扱い
前述したように 3 次元ハイブリッドコードでは計算負荷を減らすために, い
くつかのイオン粒子を超粒子とおく. 実際の粒子数を n, 超粒子数 ns をとすると粒子数比 ration=ns/n を用いて計算を行った. 流体コードから粒子にするために,イオンの数密度を用いて粒子として配置していく.
また超粒子という仮定をおいたため, 超粒子数の電荷・質量を定める必要があ
る. 2 次元放射流体コードの i 番目の格子の質量密度 ρ(i), イオン数密度 ni(i), イ
オン電荷 valence(i), 体積 v(i) [(3.3.3)式より]から(3.3.1), (3.3.2)式のように超粒子の平均電荷 valences・平均質量 masss を求めた.
- 24 -

𝑛𝑣𝑎𝑙𝑒𝑛𝑐𝑒𝑠 = [∑ 𝑛𝑖 (𝑖) × 𝑣(𝑖) × 𝑣𝑎𝑙𝑒𝑛𝑐𝑒(𝑖)] /𝑛𝑠(3.3.1)𝑖=1𝑛𝑚𝑎𝑠𝑠𝑠 = [∑ 𝜌(𝑖) × 𝑣(𝑖)] /𝑛𝑠(3.3.2)𝑖=1③2 次元座標系 z-r から 3 次元座標系 x-y-z への変換
前述したように物理量は軸対称性(z 軸対称性)を念頭に設定してあるので, 物
理量の軸対称性を仮定すると, 2 次元放射流体コードの計算データは Fig.3.8 のように円柱形に拡張できる. ある一つの計算格子(i,j)の拡張を考えると Fig.3.9 のようになる.

Fig.3.8次元放射流体コードの計算データ拡張のイメージ

Fig.3.9 計算格子(i, j)とその拡張のイメージ- 25 -

この要素の体積 v(i, j), イオン粒子数 n(i, j)は(3.3.3)式, (3.3.4)式のようになる.
(3.3.3)𝑣(𝑖, 𝑗) = 𝜋,𝑟𝑏 2 (𝑗 + 1) − rb2 (𝑗)- × ,zb(𝑖 + 1) − zb(𝑖)(3.3.4)
𝑛𝑠 (𝑖, 𝑗) = 𝑛𝑖𝑐 (𝑖, 𝑗) × 𝑣(𝑖, 𝑗) × 𝑟𝑎𝑡𝑖𝑜
ここでπ は周期率で, zb(i), rb(j) は z 方向 i 番目・r 方向 j 番目の格子の境界位
置座標(Fig.3.9), nic (i, j)はその格子中心で定義されるイオン数密度[m-3], ratio は
前述した粒子数比である.
今ここで整理すると, Fig.3.9 のような格子を回転させたリングごとの超粒子数 ns, 格子境界の z, r 方向それぞれの速度 vz(i, j) , vr(i, j), 上述した格子の位
置座標が分かっている.
3 次元ハイブリッドコードではイオンを粒子として取り扱うので, 各格子の粒子数だけ(3.3.8)式のようにイオンの位置座標(x, y, z), イオン速度(vx, vy, vz)
を設定する. ここで分布が均等になるように乱数(call random number)を使って拡張回転角 θm を(3.3.7)のように定義した. 下付き添え字 m は超粒子数 ns の中の m
番目の粒子であることを示している.
𝑧𝑚 = 𝑧𝑏(𝑖 − 1) + *𝑧𝑏(𝑖) − 𝑧𝑏(𝑖 − 1)+ × 𝑟𝑎𝑛𝑑𝑜𝑚 𝑛𝑢𝑚𝑏𝑒𝑟𝑚
𝑟𝑚 = 𝑟𝑏(𝑗 − 1) + *𝑟𝑏(𝑗) − 𝑧𝑏(𝑗 − 1)+ × 𝑟𝑎𝑛𝑑𝑜𝑚 𝑛𝑢𝑚𝑏𝑒𝑟𝑚
𝜃𝑚 = 2𝜋 × 𝑟𝑎𝑛𝑑𝑜𝑚 𝑛𝑢𝑚𝑏𝑒𝑟𝑚(3.3.5)(3.3.6)(3.3.7)x𝑚 = r𝑚 × cosθ𝑚{ ym = rm × sinθmzm = zm(3.3.8)v𝑥𝑚 = vr (i, j) × cosθi{ vym = vr (i, j) × sinθivzm = vz (𝑖, 𝑗)(3.3.9)

ここでこれらの式のイメージ図を下記する.
Fig.3.10 は式(3.3.5), (3.3.6)を示しており, z と r 方向の超粒子の位置を乱数で散りばめている. さらに Fig.3.11 は式(3.3.7), (3.3.8)を示しており, r 方向の超粒子の位置を乱数によって決めた角度 θ を使って x と y に分散させている. 前述したように, 式(3.48)によって各リングに入る超粒子の数が決まっているため, その数の超粒子の位置がこれらの式によって決められ, Fig.3.12 のように超粒子配置されることになる. そして次の格子（リング）でまた同じように超粒子の位置を決定する. 1次元の場合は z 方向の情報をランダムな角度である φ を加えて、リングではなく球殻にイオンを配置するイメージとなる。
次に速度に関しては z 方向はそのままであるが, x, y 方向は Fig.3.13 のように乱数で決めた角度 θ を使って分散させている.
- 26 -

Fig.3.10 z-r 方向の粒子位置振り分け

Fig.3.11 θ 方向（x-y 平面）の粒子位置振り分け

Fig.3.12

カップリング後粒子位置の概念図

Fig.3.13

x-y 方向速度ベクトルの概念図- 27 -

第4章

計算概要
4.1 体系
磁気スラストチャンバーシステムのイオン計測のため, ファラデーカップ型のチャージコレクタを用いて大阪大学の EUV データベースレーザー施設にて実験が行われた. 数値解析においてその結果と比較するため, 計算体系を実験と同様にした. その体系を説明するため, まずはその実験概略図を Fig.4.1 に示す.

Fig.4.1

イオンビーム電流計測実験概略図

Fig.4.1 において, アルファベットの長方形がチャージコレクタであり, 中心の円がターゲット位置. その隣の長方形がコイルまたは磁石を表している. 図の左側からレーザーを照射して, アブレーションしたプラズマを下流側へ押し返して推力とする（レーザー入射方向を上流, 反対方向を下流と以下表現する）.
最初に計測点について説明する. 実験ではチャンバー内にチャージコレクタ
を[A, A’, B, E, F, G, H]の 7 点に置いた. ターゲットからの距離は A’点のみ
200mm で, 他点は 150mm である. また, レーザー軸（A と A’の間）を0°とし（以下基準線）, 上側をマイナス, 下側をプラスと考えると各点それぞれ,[−14°, 14°, 45°, 60°, 75°, 104°, −162°-である. 数値解析においては, ターゲットからの距離 150mm, 基準線からの角度[30°, 90°, 120°, 135°, 150°]の点をさらに加えて計算した.

- 28 -

続いて計算体系の概略図を Fig.4.1 に示す. 実験において磁場をネオジム磁石とコイルによって作ったため, シミュレーションにおいてもそれに倣い, ほぼ同じ条件の初期磁場をセットする. ネオジム磁石は中がくり抜かれた中空のものを使用し, 外径が16𝑚𝑚, 内径が8𝑚𝑚, 長さは20𝑚𝑚である. コイルは外径が80𝑚𝑚, 内径が52𝑚𝑚, 長さは10𝑚𝑚で, 磁石, コイル共に X-Y 平面に平行である. ターゲットは半径250𝜇𝑚の球状ポリスチレンを使用. 磁石やコイルとターゲット（原点）との距離は実験値に合わせ, 磁石が8𝑚𝑚, コイルが14𝑚𝑚離れている. また, ターゲット（原点）における磁場強度はそれぞれ0.12𝑇と0.11𝑇である. また同じパラメータの磁石を使い, ターゲットから35𝑚𝑚離すことでターゲット（原点）位置で磁場がほぼ無いとした条件（以下「磁場なし」）を設けている. さらに今年 2014 年 2 月末に行われる実験ではターゲット位置で 1T の磁場強度となるようコイルを新たに作成する. これに合わせ, より大きな磁場中のイオン挙動を見るためターゲット位置で 1T, 10T の条件でも数値解析を行う.
このコイルは外径が50𝑚𝑚, 内径が26𝑚𝑚, 長さは10𝑚𝑚である. 今回行った計算パラメータと条件（コイル, 磁石とターゲットとの距離, ターゲット位置における磁場強度）をそれぞれ Table. 1 と Table.2 にまとめた.

Fig.4.2 計算体系- 29 -

Table.1

計算パラメータTime step 数

10000

Time step 幅 [ps]50

計算時間 [ns]

5000

超粒子数

1000000

Table.2

シミュレーションの磁場条件

磁場発生源

距離 [mm]

磁場強度 [T]

コイル

14

0. 11

コイル

13

1. 0

コイル

13

10

磁石

8

0. 12

磁石

35

磁場なし

次にレーザーパラメータを Table.3 に示す.
Table.3 レーザーパラメータ

4.2

エネルギー

4 ,𝐽-

波長

1.064 ,𝑛𝑚-

パルス幅

15 ,𝑛𝑠-

集光径

500 ,𝜇𝑚-

イオンビーム電流の計算

実験と同様の位置においてシミュレーション内でイオンビーム電流を計算するために, 各計測点に面（チャージコレクタの大きさ程の円）を仮定した. その面をある時間内に通り過ぎた超粒子の数を計算し, そこから実際のイオン粒子の数を出して電流を算出した（Fig.4.3 参照）. 面の向きはターゲットを向いており, 単位は A である.

- 30 -

イオンビーム電流の式は次の通りである.

𝐼𝑜𝑛 𝐵𝑒𝑎𝑚 𝐶𝑢𝑟𝑟𝑒𝑛𝑡 =𝑍 ∙ 𝑒 ∙ 𝑛𝑢𝑚𝑏𝑒𝑟
𝑡𝑖𝑚𝑒

(4.1)

Fig.4.3 イオンビーム電流計算の概略図

4.3 力積の計算
推力は以下のように Impulse bit で表す. これは 3D ハイブリッドでの初期運動量を 0 とした運動量変化である. 単位は [Ns].
impulse bit = ∆ [∑ mvz ] = ∑ mvz − ∑ mvz0
∑ mvz ：推進方向の運動量総和∑ mvz0 ：推進方向の初期運動量総和- 31 -

(4.1)

第5章

結果と考察
5.1 イオンビーム電流比較実験においてイオン粒子が磁場と作用してどのように広がっているのか, 直接イオン粒子を可視化することはできないため, 本研究室ではチャージコレク
タを用いて飛んできたイオン粒子の量を計測して, 間接的にどの範囲に広がっているのかを調べている. チャージコレクタは膨張するプラズマイオンを導体で受けることにより, イオン電流を時間の関数として記録できる. また小型で容易に量産でき, 多チャンネルによる空間分布の計測が行えるというメリットを有し, レーザー生成プラズマ計測のための一般的計測器として普及した[23].
そこで数値解析においても実験結果と比較していくためにイオンビーム電流を各測定点で計算した. その結果をこの節で示していく.

5.1.1 実験, 1D, 2D の比較
まず, 実験結果と数値解析結果の比較について述べる. 上述したように実験
ではどの角度の測定点（チャージコレクタ）にイオンが良く飛んできているかを調べることで, イオン粒子の広がり方を探っている. そこで, 各測定点でのイ
オンビーム電流の積算値を求め, それを比較する. 以下に示すのは, 磁石を使用し, ターゲットとの距離を 8 [mm]と設定した時の結果である. ターゲット位
置での磁場強度は 0.1 [T]である.

- 32 -

Fig.5.1

Fig.5.2

数値解析における各計測点

各測定点におけるイオンビーム電流比較

Fig.5.2 は縦軸をイオンビーム電流（単位：[A]）で対数目盛を使用している. 横
軸は基準線（0 [°]）からの角度を表している（単位：[°]）. 青の丸が 1 次元を使用した結果, 四角い茶色のマーカーが２次元を使用した結果, 三角の灰色マー
カーが実験結果である. Fig.5.2 を見て分かる通り, 実験結果においては角度の絶対値が小さい領域, これを以下低角度域と表記するが（具体的には Fig.5.1 参照）,
この領域に飛んできているイオン粒子が他の角度に比べて尐ない. これは, 低
角度域がターゲット表面でレーザーが照射されていない裏側にあたるため, 照- 33 -

射側よりもアブレーションするプラズマの量や速度が非常に小さくなっているためである（以下, 背面照射）. 数値解析の結果において、180 [°]付近で値が下がっているのは、磁石にイオンが当たったため、チャージコレクタまで到達していないからである。Fig.5.3 を見るとよく分かる。
これに対して 1 次元の放射流体コードを使用した計算結果を見てみると, 同
じ低角度域に飛んできているイオン粒子は他と比べてあまり大差がない. つまり背面照射が再現できていないという事である. その理由は, 第 3 章でも尐し触れているが, 1 次元放射流体コードではカップリング時に 1 次元軸上のイオン
粒子を 360 度の様々な方向に散りばめ, レーザーの照射方向を考えずに計算しているからである. そのため, 飛んできているイオン粒子の量に角度依存性がない.
ここで 2 次元放射流体コードを用いた計算結果を見ると, 他の角度に比べて低角度域にはイオンが飛んできていないことが分かる. 使用する放射流体コードを２次元にしたことでレーザーの照射方向を考える事ができ, 背面照射を再現することができている. よって 1 次元放射流体コードから２次元放射流体コードに変更したことで, イオン粒子挙動の角度依存性を数値解析でも表現し再現性向上を実現できた. Fig.5.2 に数値解析によるイオン粒子の動きを示した. これを見ると 1 次元と 2 次元両者の違いがよく分かる.
(a) 1 次元(b) 2 次元- 34 -

Fig.5.3

磁石（ターゲットとの距離 8 [mm]）条件下でのイオン粒子の動き

Fig.5.3 の左側と右側はそれぞれ１次元と２次元の放射流体コードを使った結果である. 黄色い長方形は磁石, 青い点はイオンビーム電流の計測点である. レー
ザーは左側X = 0 ,m-の部分から入ってきている. 真ん中右側の突き出した部分は,レーザー照射の真裏にあたるアブレーションプラズマが広がったものである.
ここで, 特に 1 番右下の図で粒子がきれいに縦に割れている理由を考察する.
この理由として考えられるのは大きく 2 つある.
1 つ目は単純に放射流体コードにおけるメッシュ数が尐なかったことである. 1
次元コードでは, 半径方向に 500 のメッシュを切っていて, イオン粒子が約 10[mm]まで広がっている. 一方で２次元コードでは計算領域が Z 方向に 5 [mm]でそれを 100 のメッシュで切ってある. 同じコードでメッシュ間隔に幅はあるにせよ,
単純に 1 次元コードの方がメッシュに関して密である. ただ, 計算コストを考えると, 2 次元コードで 1 次元並みのメッシュ幅にするとかなりの時間を要することが想定される. メッシュの切り方に関しては, 今後議論, 試行の余地があると言える.
2 つ目の理由はカップリング時におけるイオン粒子の散りばめ方の違いである.
イメージに関しては第 3 章で既に記述しているが, 1 次元はある距離の格子幅を持った球殻に粒子を乱数で散らせる. 2 次元は格子を格子幅だけ持ったリング状に散らせる. つまり 1 次元はどの方向から見てもある程度滑らかに粒子が分散しているが, 2 次元は X-Y 平面で見ればある程度滑らかだが, X または Y-Z 平面で見ると,100 個のメッシュで区切られている. これらの理由から, 右下の図のような現象が起きていると考察できる. Fig.5.4 と Fig.5.5 にカップリングした後のイオン粒子
分布を示す.

- 35 -

Fig.5.4

カップリング後の Z-X 平面のイオン粒子位置分布

Fig.5.5

カップリング後の X-Y 平面のイオン粒子位置分布

Fig.5.4 をよく見ると, X 軸方向（縦）に多くのラインが入っているのが分かる. Z
軸方向に関してはラインが入っていない. Fig.5.5 を見ると, まんべんなく粒子が散っているのが分かる. カップリングコードでデータを渡す際にメッシュを切りなおして間のデータ（イオン粒子数密度）を補間すればもう尐し滑らかになるはずなので, これから改善の余地がある.

- 36 -

5.1.2 磁場による比較
次に数値解析における磁場依存性を見る. 実際に磁気スラストチャンバーシ
ステムを使うためにはシミュレーションコードの開発は非常に重要であり, その妥当性検証に実験との比較は必要不可欠である. その比較のため, 磁場依存性を見ることもまた重要になってくる. 本研究室ではターゲット位置の磁場強度が 0.1 [T]または磁場がない（ほぼ 0 [T]）という条件に加えてターゲットとコ
イルや磁石の距離を変えることで多尐磁場強度を変化させて実験を行った. そしてこれから先の実験において磁場強度 1 [T]以上における実験を目指している.
そこで数値解析においても, 磁場の強さを変えて結果を出してイオンビーム電流の磁場依存性を調べる. この結果を後の実験結果と比べていくこととなる.
以下に示すのは, 0.001 [T], 0.1 [T], 1 [T], 10 [T]におけるイオン粒子の位置分布である. 0.001 [T]以外はコイルの条件である（第 4 章参照）.

0.001 [T]0.1 [T]

1 [T]

10 [T]

Fig.5.6２次元放射流体コードを使用した結果- 37 -

Fig.5.6 を見て明らかなように, ターゲット位置の磁場強度が高くなるごとにイオン粒子の拡散が抑えられている. これは以前の研究における数値解析の結果とも合致し, ２次元コードを使用した場合にも同じであることが判明した.
0.1 [T]のコイルは昨年用いたコイル, 1 [T] , 10 [T]のコイルは今年２月に使用予定のコイルである. 0.001 [T]において上流側（黄色い磁石の左側）に粒子がないのは, イオンが磁石に当たっているからである. シミュレーション内では磁石やコイルに当たると跳ね返らずにそのまま消滅するとしている.
そこでこれらイオン粒子の広がりを実験結果と比較していくため, イオンビーム電流として出した結果を以下に示す.

30 [°]

162 [°]- 38 -

Fig.5.7

イオンビーム電流の時間発展における磁場依存性

縦軸が電流, 横軸が時間である. イオン飛来の時間発展における特徴をより早く見るためにチャージコレクタの位置を 150 [mm]ではなく, 50 [mm]としている. Fig.5.7 上図（30 [°]）を見ると, 1 [T]においては一度多くのイオンが飛来した後, もう一度同じように飛来している. これは最初の山は, 速度が大きく,
また初期の速度ベクトルが下流側に向いていたイオンが比較的早い段階でチャージコレクタまで到達していることを示している. そして 2 つ目の山は初期の速度ベクトルが上流側であったイオンが磁場に押し返されて下流側へ流れた,
または速度の遅いイオンの飛来を示している. よってこのグラフにより, 磁場の押し返しが起きているつまり磁気スラストチャンバーシステムが機能していることが分かる. 低角度域においてこのように 2 つの山ができる傾向は, 昨年
度の実験においても見られたもので, 合致している.
一方で 162 [°]のグラフでは, 0.001 [T]では非常にイオンが多く飛来しているが,他の条件では尐ない. やはりここからも磁場がほぼない状態であれば上流側に飛来しているはずのイオンが磁場の影響で押し返されていると言える.
この２つのグラフにより, 磁場がほとんどかかっていない状態と比べて高角
度域から低角度域にイオンが絞られていることが考察できる. 30 [°]や 162 [°]だ
けではなく, 他の角度についても同様の結果を出すことで実験結果と比較していくことが可能となり, そこからイオン粒子の広がりも分かることとなる.

- 39 -

Fig.5.8

イオンビーム電流積算の磁場依存性（2 次元放射流体コード使用）

Fig.5.8 は 2 次元放射流体コードを使ったイオンビーム電流を計算終了までの値を積算し, 角度比較を磁場ごとに示したものである. 0 [°]に近づくほどイオン粒子が広がっていないと言える. 1T や 10T の場合は明らかに粒子の広がり方が0.001 [T]と比べて狭い. 次に 0.1 [T]の条件では 0.001 [T]と比べて 75 [°]～120 [°]のイオン粒子が増えているため, 磁場で押し返されていると言える. 1 番右端の部分の値が高いのは, Fig.5.6 を見て分かる通りコイルの穴を通って上流側へ流れている粒子が多いためである.
以上のようなことがイオンビーム電流の結果から分かった. ただし, 定性的には実験結果からも分かる背面照射が再現できているが, 定量的にはワンオー
ダーずれている. それについて考察する. 今回の実験では発生した遅いプラズ
マはノイズに埋まりほとんど計測できていない. しかし数値解析では遅いプラ
ズマも電流として計算している. 最大イオン速度は約 50～200 [km/s]ほどで, この値は数値解析と実験とオーダーが合っており矛盾していない. そのため, チャージコレクタの位置まで飛んできているイオン粒子の量が数値解析よりも実験が小さいと考えられる. よって今回の結果でも実験の方が数値解析よりワン
オーダー電流値が低いと考察できる.
さらに, Fig.5.2 からも分かるように実験の方が数値解析（2 次元）よりも低角度域に多尐絞られている点についても考察していく. 数値解析（3Ｄハイブリッ
ドコード）では第 3 章でも記述したようにイオンを粒子として, 電子を流体として解くため, クーロン力は考えていない. だが実験では, 磁場によって押し返- 40 -

されたのではなく, 磁場にトラップされた電子とのクーロン力でイオンの初期
速度ベクトルが変化し, イオンの広がりが抑えられている. 数値解析よりも実験の方が低角度域にイオンが飛来しているものと考えられる. これは特に弱い磁場中で見られる傾向である.
これらの事から, 3Ｄハイブリッドコードだけではなく, 弱い磁場中で数値解析を行う場合などは, 今後 Full-PIC にしていくことを視野に入れる必要がある.
また, イオンビーム電流の計算に関しては, 現在の方法だけではなく, 各格子で計算される電流密度を利用した方法を使えば, 「点」における電流を出すことができ, 「面」に入ってくるイオン粒子を計算するよりも厳密な値が出せると考えられる. 今後このような改善を行っていくことで, よりシミュレーションの再現性が高まるものと考える.

5.2 力積比較
２次元コードを使った場合でも強い磁場をかけるほどイオン粒子を押し返せることは前節にて示されているが, ここではレーザー推進としてどれほどの力積を出せるのか, Impulse bit の磁場依存性をグラフとして具体的に Fig.5.9 に示す.
横軸はターゲット位置における磁場強度（単位：[T]）で, 縦軸は計算終了時点（5000 [ns]）での力積（単位：[Ns]）である. このグラフを見ると 1 [T]と 10 [T]では値があまり変わらない. よって今後の実験において, まずはターゲット位
置の磁場強度 1 [T]で行うことを目指すべきだと言える.

- 41 -

Fig.5.9

力積（Impulse bit）の磁場依存性

第6章

結論と今後の展望
最後に本論文のまとめ, 結論を記す.
本研究では磁気スラストチャンバーシステム内のイオン挙動について「放射流体コード」と「3Ｄハイブリッドコード」の 2 つを用いてシミュレーションを行った. 2 次元の放射流体コードを用いて次元を上げることで, プラズマアブレー
ションをより精度良く計算した. これは数値解析におけるイオン挙動の再現性- 42 -

向上を目的とし, 今後実験との比較を進めていくことを目指したものである.
以下に結論と展望を列挙する.


放射流体コードを 1 次元から 2 次元にしたことで背面照射を再現でき, 実験結
果に定性的に近づくことができた

磁場強度を上げる程イオンの広がりを抑えて力積につなげられ, 昨年度の実験でも確認できた磁場の押し返しを, 2 次元コードを用いた場合でも確認することができた. ただし 1 [T]と 10 [T]ではこれらの磁場をターゲット位置にかけても力積はほとんど変わらないという傾向が見えたため, 実験においてまずは 1T の磁場を目指すべきであることが判明した

粒子の挙動の中で不自然な列をなすのを解消するために, 2 次元放射流体コードから 3Ｄハイブリッドコードへの受け渡しの際に, 新たなメッシュを設けて前後のイオン数密度を反映させた傾斜を作るなどの工夫が必要である

参考文献[1] http://www.jaxa.jp/projects/sat/sprint_a/index_j.html[2] 中島秀紀, ：プラズマ・核融合学会誌, Vol. 83, 264 pp., (2007).
[3] 内田成明：レーザー推進の地上打ち上げおよび宇宙空間での応用, プラズ
マ・核融合学会誌, Vol. 83, 271pp, (2007).
[4] 安藤晃：レーザー推進, プラズマ・核融合学会誌, Vol. 83, 276 pp. , (2007).
[5] 内田成明：レーザー推進, プラズマ・核融合学会誌, Vol. 81, 186pp, (2005).
[6] NASA Jet Propulsion Laboratory,
Voyager The Interstellar Mission,- 43 -

http://voyager.jpl.nasa.gov/index.html,(2012,

December).

[7] Orth, C. D.: VISTA - A Vehicle for Interplanetary Space Transport Application
Powered by Inertial Confinement Fusion, UCRL-LR-110500, (2003).
[8] H. Nagatomo, T. Jozaki, A. Sunahara, K. Mima:”Integrated Implosion
Code-Radiation Hydrodynamic Code”: J. Plasma Fusion Res Vol82, No3,
(2006), p141-144[9] 砂原淳：次世代半導体リソグラフィー用極端紫外発光（Extreme Ultra-Violet）の放射流体シミュレーション, 第 10 回シミュレーション・サイエンス・シン
ポジウム＆平成 17 年度核融合科学研究所共同研究「大型シミュレーション研
究」合同研究会, (2006)[10] 大西直文：輻射流体シミュレーションの基礎, プラズマ・核融合学会誌 88,467 (2012)[11] 砂原淳, ：EUV 高原や ICF におけるレーザープラズマ相互作用の初期過程―レーザーアブレーション―, プラズマ核融合学会誌 89, 419 (2013)[12] S.P.Lyon and J.D.Johnson, Los Alamos National Laboratory Tech. Rep. No.
LA-UR-92-3407, (1992)[13] 山本和之：数値天体物理学サマースクール（流体・磁気流体コース）講義テキ
スト 「近似的リーマン解法－流体計算－」, 千葉大学, (2001)[14] 村田貴紀：銀河風研究のための流体シミュレーションコード開発, 2013 年度
第 43 回天文・天体物理若手夏の学校 資料, (2013)[15] 水田晃：相対論的流体方程式の数値的解法Ⅱ, 講義録[16] Takahiro Miyoshi, Kanya Kusano:”A multi-state HLL approximate Riemann solver
for ideal magnetohydorodynamics”, Journal of Computational Physics 208, p
317-323, (2005)[17] Shengtai Li, Journal of Computational Physics 203, 355, (2005)[18] 川渕亮：“磁場による慣性核融合炉ビームポート保護の可能性”：九州大学総合理工学府修士論文(2008)[19] 松田伸夫：“整形ターゲットを用いたレーザー核融合ロケットの推力ベクトルの制御技術”：九州大学総合理工学府修士論文(2009)[20] E.Horowitz, et al, J.Comp.Phys., 84, 279(1989)[21] A.Mankofsky, R.N.Sudan, J.Denavit:J.Comp.Phys:70(1987), 89[22] (20) C.K.Bridsall, A.B.Langdon:”Plasma Physics via Computer Simulaiton”:IOP
Publishing Ltd, (1997)[23] 日永智之, 磁気スラストチャンバーの原理検証実験,
大学院, (2013).

- 44 -

修士論文,

九州大学

謝辞
本研究において、あらゆることを御指導頂いた中島秀紀教授に深く感謝致します。
研究室生活において色々と忙しい中、身近で御指導して頂いた山本直嗣准教授に深く感謝致します。快適な研究環境を提供して下さった馬渡隆子秘書に深く感謝致します。
本研究にて使用させてもらい、また多くの助言など頂いた、放射流体計算コードに携わっていらっしゃる砂原淳先生、城崎知至先生に深く感謝致します。

- 45 -

同じ LFR 組として研究ついての多くのアドバイスを頂いた安永将人氏、秀平靖磨
氏、日永智之氏、三輪祐大氏に深く感謝いたします。大阪大学での実験を共に行い、長い大阪生活を共にした齋藤直哉氏、川島諒祐氏に深く感謝致します。研究についてのアドバイスをしていただいたり、楽しい研究室生活にしたりして頂いた諸先輩
方、下川床潤氏、杉田健策氏、平岡義章氏、中村祐輔氏、平野賢治氏、廣池匠哉氏に深く感謝致します。
研究にソフトボールやサッカー、飲み会など日々苦楽を共にした伊藤匠氏、坂井寿
光氏、木邊厚視氏、楠田将士氏、田中啓介氏、冨永真人氏、牛尾康一氏、前田優士
氏、山口敦氏に深く感謝致します。
遠いインドの地からこの研究室に来てくださった Indranuj Dey 氏に深く感謝いたします。また台湾から来て下さった李後毅氏に感謝致します。
同じ学年としてアドバイスし合い、共に頑張ってきた豊田裕司氏、笹川裕太郎氏に深く感謝致します。
忙しい仕事の合間を縫って私のコードや研究についてアドバイス頂いた前野旭弘
氏に深く感謝致します。
最後に、これまで支えてくださった家族、友人に深く感謝致します。

- 46 -