J. Plasma Fusion Res. Vol.92, No.11 (2016)７
８５‐７８６

プロジェクトレビュー

核融合科学研究所における数値実験炉研究プロジェクト
Numerical Simulation Reactor Research Projectat the National Institute for Fusion Science

１．はじめに
堀内利得
自然科学研究機構核融合科学研究所（原稿受付：２
０
１
６年７月１
３日）

１）数値実験炉研究プロジェクト体制について

平成元年（１９８９年）５月２９日核融合科学研究所の創設に際し，その研究推進の母体となる内部組織として研究所設

数値実験炉研究プロジェクトでは，これまでの炉心プラ立の主目的である大型ヘリカル装置（LHD）を建設し，そズマから周辺プラズマまでの研究を，実験や理論と連携しの装置を用いて実験的に核融合プラズマの学理の構築に貢ながらさらに発展させ，核融合磁場閉じ込め装置全体のプ

献する「大型ヘリカル研究部」と並立して，特定の閉じこ

ラズマ挙動を予測することのできるヘリカル型核融合装置

め概念にとらわれず広い視野に立って核融合プラズマの理に対する数値実験炉の構築をめざした大規模シミュレー

論体 系 化 を め ざ す「理 論・シ ミ ュ レ ー ショ ン 研 究 セ ンション研究をその目標として研究活動を展開してきた．

ター」が設立された．

このプロジェクトの推進に当たっては，大きな課題が３

この設立理念を基調としながらも，その後のシミュレー

つある．まず，（１）炉心プラズマからプラズマ対向壁でのション研究を取り巻く環境の変化や要請にこたえる最適のプラズマ物質相互作用に至るまでの核融合プラズマを構成

研究環境と研究体制を構築しつつ，シミュレーション研究する各種物理要素の導入とその物理の解明である．次に，を推進してきた．平成１６年度の法人化では，大学利用機関（２）数値実験炉を構成する各要素物理の数値モデルの高精

法人自然科学研究機構の一員として，中期計画・中期目標

度化と要素物理を組み上げて数値実験炉の構築へ向けたを策定し，共同利用・共同研究の推進に努めることを求め

コード体系化に向けた研究が必要となる．
（３）この開発研られるようになった．さらに，平成２２年度から始まった第究において，他の研究グループで開発されている同種のシ

２期中期計画では，核融合原型炉に向けた学理の体系化を

ミュレーションコードとのベンチマークテストを通じての加速するために，研究者のライン組織として単一研究部の理論検証，および，具体的な実験データとの比較による数

ヘリカル研究部に一本化し，それぞれのタスクを持った大値モデルの正当性・適応性の検証が重要となる．

型ヘリカル装置（LHD）計画，数値実験研究，核融合工学

これらの課題の解決に向けた研究活動を柔軟にかつ機能研究からなる３研究プロジェクト体制を導入した．

的に遂行できる体制として，図１に示した研究体制を構築

この第２期中期計画において，「核融合プラズマ閉じ込した［２，３］．まず，研究推進の基礎母体として機動性と柔軟めの物理機構解明，その体系化及び数値実験炉の構築」を性を持ち合わせた８つのタスクグループを組織し，このグ

数値実験研究プロジェクトの目標として掲げ研究活動を開

ループを中心となって組織・運営を行う幹事を配置した

始した．平成２６（２０１４）年７月１日には，核融合原型炉の基（表１）．これらのグループ活動では，複数のグループへの本設計に数値シミュレーションの立場から貢献することを参加，グループ間の研究交流，実験グループとの研究交流目的とした数値実験炉の研究をより一層活加速するため，が推奨されている．次に，プロジェクト活動の全体の検

数値実験研究プロジェクトの名称を「数値実験炉研究プロ

証・調整，各タスクグループの進捗状態のモニターリング

ジェクト」と変更し，第２期中期計画期間の研究活動を続と必要に応じた目標設定や研究体制の修正等を行う場としけた．第２期中期計画期間の最終年である平成２
７年度に，て推進会議を設け，定期的に必要な議論・決定を行ってい第２期中期計画全体に亘るプロジェクト研究活動に対し外る．推進会議は，議長としての研究総主幹，シミュレー

部評価を受けた［１］．本稿では，第２期中期計画期間全体

ション研究系の二人の研究主幹，およびタスクグループのでどのような活動がなされ，どのような研究成果が得られ幹事で構成され，必要に応じて他のメンバーの参加も要請たかを報告する．

している．また，プロジェクトの推進には，LHD 計画プロ

1. Introduction HORIUCHI Ritoku
National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: horiuchi.ritoku@nifs.ac.jp

７８５!2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

表１

タスクグループ体制．

タスクグループ名
プラズマ流体平衡・安定性
高エネルギー粒子
統合輸送シミュレーション
新古典・乱流輸送シミュレーション
周辺プラズマ輸送
プラズマ壁相互作用
多階層複合物理
シミュレーション科学基盤
図１

幹事
市口 勝治
藤堂 泰
横山 雅之
菅野 龍太郎
鈴木 康浩
中村 浩章
三浦 英昭
大谷 寛明

メンバー
１
４名
５名
２
１名
１
４名
９名
１
８名
１
８名
１
６名

数値実験炉研究プロジェクトの研究体制．

ジェクトや核融合工学研究プロジェクトとの連携が重要とする核融合プラズマや基礎プラズマの数値解析，中小規模なる．この２つのプロジェクトとの連携の窓口として，担のシミュレーション，装置設計やデータ処理に関する共同当者を配置し，研究交流や連携研究の推進を図っている．

研究を支援することを目的として運用される LHD 数値解

また，プロジェクト全体での意見交換・情報交換，グルー析サーバ（１４４ノード）と，これを支えるフロントシステ

プ研究活動報告，推進会議報告，出張報告や国際会議の報

ム，データ解析を目的としたデータ解析サーバ，可視化処

告，LHD 計画プロジェクトや核融合工学研究プロジェクト

理サーバ等で構成されている．日常的なプラズマシミュとの研究交流の場として，全体会議を設定している．さらレータの運用，関連研究環境の整備や共同研究者への様々には，国際連携・国内連携活動においても担当者の配置やなシステム利用支援は，プログラム開発支援室の活用や計国際委員等の派遣を行うなど，COE として活動を積極的に

算機作業班等の活動を通じて行っている．その結果，プラ推進している．これら共同研究に関する活動状況の概要にズマシミュレータを用いた共同研究に参加する共同研究者

ついては，平成２５（２０１３）年１月号のプロジェクトレビュー．
の人数は常に１５０名余となっている（図２）の中で既に報告してある［３］．

以下の２章と３章では，核融合プラズマシミュレーショ

２）プラズマシミュレータの整備

ンとその支援研究である基礎物理シミュレーションに分け

数値シミュレーションがその主たる研究手段である数値て，大規模並列型計算サーバを用い国際・国内共同研究を

実験炉研究プロジェクトにとって，プロジェクト推進に最通じて得られた研究成果の概要を紹介する．４章では全体適なスーパーコンピュータの導入・更新と関連研究環境ののまとめと今後の展望を記す．

整備が必須となる．これまで，共同研究用のマシーンである大規模並列型計算サーバ（プラズマシミュレータを含参考文献

む）の定期的な更新・性能向上により，最先端のスーパー［１］核融合科学研究所運営会議外部評価委員会：「平成２７年
度外部評価報告書」
（２
０
１
６年３月）
．
［２］R. Horiuchi et al., Plasma Fusion Res. 6, 2101055 (2011).
［３］R. Horiuchi, J. Plasma Fusion Res. 90, 2 (2014).

コンピュータを用いた研究環境を数値実験炉研究プロジェ
クトを中心とする共同研究に提供してきた．表２に示しているように，平成２４（２０１２）年１０月の中間レベルアップでは，日立製作所製 SR１６０００モデル L２１２８ノードを同社製
SR１６０００モデル M１（３２２ノード構成）へと約４倍の性能向上を，平成２７（２０１５）年６月の更新 で は，さ ら に 約８倍の2.62 Pflops の演算性能と８１TB の主記憶を持った富士通製
PRIMEHPC FX１００（２５９２ノード構成）を中心とするシステ
ムを実現した．
その結果，世界のスーパーコンピュータ性能の指標であるTop５００において世界２７位，もう一つの指標であるHPCGにおいて世界１２位を記録した（ともに平成２７（２０１５）年７月時点）．大規模並列型計算サーバは，数値実験研究プロ
ジェクトの大規模並列計算のための計算機 プラズマシ
ミュレータ （２４４８ノード），LHD プラズマを中心に，関連

図２
表２
期間
H２
１．
３〜H２
４．
８
H２
４．
１
０〜H２
７．
３
H２
７．
６〜H３
１．
５

機種
SR16000 モデル L2
１
２
８ノード
SR16000 モデル M1
３
２
２ノード
PRIMEHPC FX100
２
５
９
２ノード

プラズマシミュレータを用いた共同研究の利用者数．

大規模並列型計算サーバの主要性能．
演算性能

主記憶容量

外部記憶装置容量

77 TFlops

16 TB

0.5 PB

315 TFlops

40 TB

2.0 PB

2.62 PFlops

81 TB

10.0 PB

７８６

J. Plasma Fusion Res. Vol.92, No.11 (2016)７
８７‐７９３

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

２．核融合プラズマシミュレーション研究
2. Simulation Researches in Fusion Plasmas

２．
１ MHD 平衡と安定性
2.1 MHD Equilibrium and Stability
市 口 勝 治１，２），NICOLAS Timothée＊１），佐 藤 雅 彦１），石 澤 明 宏３），石 崎 龍 一１）
ICHIGUCHI Katsuji１，２）, NICOLAS Timothée＊１）, SATO Masahiko１）, ISHIZAWA Akihiro３）and ISHIZAKI Ryuichi１）
１）

自然科学研究機構核融合科学研究所，２）総合研究大学院大学，３）京都大学エネルギー科学研究科（原稿受付：２
０
１
６年７月１
３日）

流体モデルのシミュレーションによって，プラズマのマクロ的挙動の解明が大きく進んだ．まず，大型ヘリ
カル装置（LHD）の平衡及び圧力駆動型不安定性に関して，解析を進めた．独自に開発したマルチスケール手法
を用いて磁気軸スイング実験での圧力崩壊現象を再現し，そのメカニズムを解明した．また，共鳴摂動磁場を印加した場合での圧力分布の変化とそれに伴うモード構造の変化を明らかにし，それによって，実験で得られている崩壊現象の局在を説明した．さらに，交換型モードの安定性に対する反磁性効果の影響を導入することにより，線型固有値解析及び非線型発展の両方の観点から，モードが電子反磁性効果に回転することを示し，実験結果を説明した．次に，微視的効果を含めることにより，磁気島に対する分極電流の駆動効果を新たに発見し，イオン
温度勾配不安定性による乱流粘性の値を初めて評価した．また，ペレットアブレーションによって形成されるプ
ラズモイドの移動方向が，その生成位置に大きく依存することを見出した．
Keywords:
magnetohydrodynamics (MHD), equilibrium, stability, interchange mode, diamagnetic effect, magnetic island,pellet ablation

２．
１．
１ LHD 磁 気 軸 ス イ ン グ 実 験 の シ ミ ュ レ ーション

こと，及び初期条件として圧力分布には実験に対応した分布を利用し回転変換分布には交換型モードが非線型飽和したときにえられる分布を採用したことが挙げられる．

LHD 実験では，真空磁気軸位置（!$#%）に対する安定限界を調べるために，磁気軸スイング放電が行われた［１］．ここのシミュレーションの結果，図１に示すように，磁気の実験では，放電中に垂直磁場を実時間制御して，対応す

軸スイングのない場合にはプラズマは安定であるが，磁気

る真空磁気軸位置を !$#%!3.6 m から 3.5 m まで磁気軸を内軸スイングを行った場合には崩壊現象が生じるという，実側へスイングさせている．このとき，!$#%!3.54 m 付近に

験結果を再現する結果が得られた．磁気軸スイングを行っ

おいて，急激な中心領域でのコア崩壊現象が観測された．

た場合では，コア領域の圧力崩壊により，中心ベータ値が!
#!!の磁場揺動が顕著に観測されていまた同時に，"!"

大きく低下するが，シミュレーションでも !$#%"3.55 m

る．そこで，このコア崩壊現象のメカニズムを解明するたで中心ベータの急激な低下が得られており，少し不安定なめに，簡約化 MHD 方程式に基づく NORM コード［２］を用

結果となっているが，実験とよく一致していると言える．

いた非線型 MHD シミュレーションを行った．このシミュ

また，図２に示すように，この崩壊現象前後での圧力分布

レーションにおいては，磁気軸スイングに対応する背景磁の変化においては，定性的には一致しているが，実験に比場の変化を取り入れることが本質的となる．ところが，こべ中心圧力低下がやや小さい結果が得られた．さらに，実の背景磁場はプラズマのダイナミクスに比べてはるかに長験では，加熱効果によってベータ値は再度上昇し崩壊が繰い時間スケールで変化している．そこで，この変化を取りり返されることが得られており，この現象も今回定性的に入れるために，これまで開発してきた平衡磁場計算と非線再現できた．ただし，崩壊の繰り返しの周期は一致させる

型ダイナミクス計算を交互に繰り返すマルチスケール手法

ことはできなかった．

を，さらに拡張した［３］．今回拡張した点は，平衡計算の際!
#!!成分が支配的
この崩壊を引き起こした摂動が "!"にダイナミクス計算による回転変換の変化分を取り入れたなインファーナルモードであることも得られている．この＊現所属：ローザンヌ連邦工科大学

corresponding author’s e-mail: ichiguch@nifs.ac.jp

７８７!2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

図２

図１

磁気軸スイング時の崩壊前と崩壊後の圧力分布．

磁気軸スイングなし（a）とあり（b）
での中心ベータ値（実線）及び対応する真空磁気軸位置（破線）の時間発展．青線
は実験値，赤線はシミュレーション結果を示す．

モード数は観測された磁場揺動モード数に対応している．
このインファーナルモードは圧力駆動型モードであり，磁気丘の影響を強く受ける．このような内寄せ磁気軸スイン
グの場合には磁気丘が強くなるため，この磁気丘の変化がインファーナルモードを不安定化したと考えられる
［３］．
図３

２．
１．
２ LHD プラズマの圧力駆動型モードに対す
る誤差磁場の影響
LHD では，誤差磁場の存在が実験的に確認されており，

RMP 印加時の平衡圧力分布の鳥瞰図と等高線．文献［５］‘Three-dimensional MHD analysis of heliotron plasma
with RMP’, K. Ichiguchi et al., NUCLEAR FUSION, vol.55,(2015) 073023 (10 pp), Figure 2 (a), Published 23 June 2015
by IOP Publishing より転載．

"!!
"#!!の成分が支配的となっている．そこで本研究で!
!"の磁場配位において，共鳴摂動
は，!%#&!3.6 m，#$!!

配位は，誤差磁場がない平衡では，交換型モードに対して磁場を含む LHD プラズマのMHD特性を三次元コードを用

不安定になる事が得られている．この場合，ポロイダル断

いて解析した．誤差磁場が存在する場合には，それによっ面において，摂動圧力のモードパターンはポロイダル方向て磁気面形状に擾乱が生じる．したがって，VMEC のようにほぼ一様に分布する．これに対し，誤差磁場が存在するな入れ子状の磁気面を仮定した平衡計算コードは適用でき平衡においては，図４に示すように，摂動が磁気島の X 点ない．また，この磁気島構造と矛盾のない圧力分布を持つ近傍に局在するようなモード構造を持つモードが最も線型平衡解を得なければならない．一方，HINT コードでは，入不安定になることを見出した．磁気島形状において，O 点

れ子状の磁気面の存在を仮定していないため磁気島を含む

近傍では圧力勾配が小さいのに対し，X 点近傍では圧力勾平衡解を計算でき，同時にこの磁場構造に対応する圧力分

配が大きい．したがって，この構造を持つことによって，布が自動的に決まる．そこで，本研究では，HINT コードを最も効率よく不安定エネルギーを利用することができる．

用いて平衡計算を行った．水平方向に一様な誤差磁場を仮さらに非線型ダイナミクスにおいては，線型モード構造

定 し て 平 衡 計 算 を 行 っ た 結 果，図３に 示 す よ う に，から予想されるように，摂動磁場によって形成される磁気"!!
"#!!の構造を持つ磁気島が存在する平衡が得られ島の X 点近傍から崩壊が生じて全体に広がっていく結果がた．このとき，圧力分布は磁気島内部で局所的な平坦構造

得られている．したがって，モード構造のトロイダル及びを持っている．

ポロイダル方向の空間的位相は，磁気島の幾何学的形状に

この平衡解における三次元安定性を MIPS コードを用い

よって固定される．一方，誤差磁場がない場合でも交換型て解析した．MIPS コードでは，HINT コードと同じ円柱座

モードが崩壊現象を引き起こすが，その場合にはそのよう標を用いているため，HINT コードによる平衡解の安定性な特別な位相は存在しない．LHD 実験においては，誤差磁をそのまま計算できるという利点がある．今回用いた磁場場がキャンセルされる場合とそのまま残っている場合につ７８８

Project Review

2.1 MHD Equilibrium and Stability

図５
図４

RMP 印加時の摂動圧力の線型モードパターン（赤及び青の領域）の三次元図と，横長ポロイダル断面でのモードパ
ターン及び磁力線ポアンカレプロット．モードパターンにおいて赤及び青は摂動圧力の符号が異なることを示す．文献［５］‘Three-dimensional MHD analysis of heliotron
plasma with RMP’, K. Ichiguchi et al., NUCLEAR FUSION,
vol.55,(2015) 073023 (10 pp), Figure 7, Published 23 June2015 by IOP Publishing より転載．

K. Ichiguchi et al.

FFHR-d1（FIS-HA-MPD）の線型成長率の粘性依存性，及び，$= 10−6 及び 3×10−4 での摂動圧力のモードパターン．

であることが示唆されている．しかし，トロイダル配位での異常粘性の評価手法はまだ確立されておらず，これよりも大きい可能性がある．その場合，成長率を低減する可能性がある．また，プラズマ中には，径電場による E×B フ
ローや反磁性フローが存在し得るが，本解析では，これらのシアフローの効果が含まれていない．以下の解析結果にいてプラズマ圧力の崩壊現象が比較されている［４］．誤差も示されているように，反磁性フローが存在すると，交換磁場がキャンセルされている場合には，崩壊時のモード構型モードに対し て 安 定 化 の 寄 与 を 持 つ 可 能 性 が あ る．

造がトロイダル方向に散乱して分布しているのに対し，誤E×B フローにも同様の安定化寄与が期待される．した差磁場が存在する場合には，トロイダル方向において局在がって，今後は，粘性及びシアフローの効果を適切に取りすることが観測されている．したがって，このモードの局入れる手法を確立し，FFHR-d1 に対して精度の高い予測を在は，本研究で見出されたメカニズムによるものと考えら行うことが必要であると考えている．

れる［５］．ただし，誤差磁場が無い場合の不安定モードの"'##であるのに対し本研究
モード数が，実験では &##２．
１．
４ LHD における理想交換型モードの電子反磁性方向回転メカニズムの解明"'#$であり，必ずしも一致していない．このでは，&#$
差異は，支配的なモードを規定する散逸パラメータの選択

本研究では，LHD における交換型不安定性の成長率と実が困難であることによるものであり，今後，数値計算手法

周波数に対する反磁性効果や熱拡散の効果を解析した．そのさらなる改善を行うことにより，より実験に即した解析の結果，これまでにない新しい性質を見出した．LHD におを進める予定である．

ける交換型モードの安定性は完全には解明がなされていな
い．特に，モードの反磁性回転による影響はその中の一つ

２．
１．
３ FFHR における三次元 MHD 安定性解析

である．LHD 実験においては，モードの回転がなくなるとヘリオトロン型 DEMO 炉として設計が進められている

同時にモードが急激に成長する現象が観測されている．そFFHR-d1 の三次元線型安定性の解析を行った．今回は，文の現象はモードロッキングと呼ばれて，圧力崩壊を引き起献［６］に示されている，FIS-HA-MPD の平衡の安定性を調こす．また，これまで知られている理論によると，反磁性

べた．三次元平衡計算には HINT コードを用い，安定性解

回転方向はイオン反磁性方向になるはずであるが，実験で析には MIPS コードを用いた．この安定性解析では，電気は，電子反磁性方向に回転していることが報告されてい抵抗，粘性，熱伝導の効果を考慮した．電気抵抗には，主る．したがって，反磁性回転による安定性への寄与とモー半径 "+!.,*+/-.#15.6 m，磁場強度 !+!.,*+/-.#4.7 T，及び文献

ドの回転方向の解明は重要な課題である．本研究では反磁［６］の Figure 2 に示されている電子密度 ',#2.5×1020 m−3，性効果を含む安定性解析を行い，交換型不安定性の成長率

電子温度 #%#11 keV に対応する値を用いた．また，磁場に及び実周波数が熱拡散係数によって変化しうることを明ら

垂直方向の熱拡散係数は，#& #1 m2/sec を仮定し，磁場にかにした［７］．

平衡方向の熱拡散係数 #'に対しては，数値安定性が確保で

ここでは，線型簡約化 MHD コードを用いて，理想交換

きる大きな値として，#'#$!
'("#"&#& を用いた．このときの MHD 不安定性の成長率の粘性 $に対する依存性を図５

型モードが不安定な平衡での線型固有値解析を行った．特に示す．ここで，成長率及び粘性は，それぞれ，アルフベ値に与える影響について注目した．まず，熱拡散が与える

ン速度 $) を用いて，$)""+!.,*+/-.及び "+!.,*+/-.$)で規格化さ影響について調べた．一般に，熱拡散は安定化の寄与を持れている．また，$#10−6−4

及び 3×10に，磁力線に垂直方向の熱拡散，粘性，反磁性効果が固有での，摂動圧力のち，その係数が大きくなるにつれて第１固有値の成長率はモードパターンも示している．この図より，$$%"#"!&

小さくなることが知られている．今回，第２固有値も同時で，交換型モードが不安定である結果が得られている．ス

に調べることによって，もっと複雑な振る舞いがあること

ラブ配位での流体乱流解析では，異常粘性として $%#"!'がわかった．図６に示すように，熱拡散係数を大きくする
７８９

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016と，第１固有値の成長率は減少するが，同時に第２固有値値は同じで符号が異なる値を持つ．ここで，正及び負の値の成長率は増加する．そして，熱拡散係数 #" がある限界

は，それぞれ，イオン反磁性方向及び電子反磁性方向に対値，#!，以上になると，両者の成長率は融合することが得応する．#" "#!では，成長率が等しく反対方向の回転周波られた．一方，モードの回転周波数に対応する固有値の実数を持つ２つのモードが共存するため，モードの回転方向

周波数は，#" !#!では，第１固有値と第２固有値は理想は決まらない．

モードのようにどちらもゼロである．しかし，#" !#!で成この状況において，反磁性効果を導入した．この場合，長率の融合が生じると同時に実周波数の分岐が生じ，絶対

図７（a）に示すように，#" "#!において２つの固有値の成長率の融合が外れて両者は分離する．また，図７（b）に示すように，実周波数も #" !#!においても分離している．この場合，第１固有値の実周波数は常に正の値を持ち，第２固有値の実周波数は常に負の値を持つ．即ち，最大成長率を持つモードはイオン反磁性方向に回転することになる．
さらに，十分大きな粘性を加えると，この状況は変化す
る．図７（d）に示すように，実周波数の傾向は粘性が無い
場合と同じである．しかし，図７（c）に示すように，熱拡散係数がある値より大きくなると，２つのモードの成長率
は，大小関係が逆転する．したがって，この大きな熱拡散係数の領域では，電子反磁性方向に回転するモードが支配的となる．このように，熱拡散と粘性とが同時に加わることによって，実験において電子反磁性方向に回転するモードが観測されていると考えられる．full-MHD 方程式を解く
MIPS コードを用いた計算結果でも，同様の結果が得られている［８］．

図６

図７

交換型不安定性の第１固有値（緑）と第２固有値（赤）の成長率（$，実線）及び実周波数（&r，破線）の熱拡散係数#⊥に対する依存性．文献［７］‘Dissipation and diamagnetic
effects on the interchange mode growth rate and rotationin reduced magnetohydrodynamics applied to stellarators’, T.Nicolas and K.Ichiguchi, NUCLEAR FUSION vol.56(2016) 026008 (15pp) Figure 3, Published 30 December
2015 by IOP Publishing より転載．

２．
１．
５ 抵抗性インターチェンジモードの伝播のシミュレーション研究
二流体効果は MHD 不安定性の安定性に影響を及ぼす
他，モードの回転を引き起こしうる．線形解析によれば，成長率や実周波数に対する反磁性効果の影響．粘性%がない場合の成長率（a）と実周波数（b），十分大きな粘性が有る場合の成長率（c）と実周波数（d）．粘性がない場合は，イオン反磁性方向に回転するモードが支配的であるが，十分大きな粘性が有る場合，成長率は逆転し電子反磁性方向に回転するモードが支配的になりうる．文献［７］‘Dissipation and diamagnetic effects on the interchange mode
growth rate and rotation in reduced magnetohydrodynamics applied to stellarators’, T.Nicolas and K.Ichiguchi, NUCLEAR FUSION vol.56 (2016) 026008 (15pp) Figure11 and Figure 14, Published 30 December 2015 by IOP Publishing より転載．

７９０

Project Review

2.1 MHD Equilibrium and Stability

K. Ichiguchi et al.

抵抗性インターチェンジモードは電子の反磁性ドリフト方向に回転し，その速度は電子の流体速度よりも小さいことがわかっている．一方，LHD 実験によれば抵抗性インター
チェンジモードは電子の流体速度と同じ速度で回転していることが観測されている．したがって，実験で観測されるモードの回転速度は線形解析だけで説明することはでき
ず，非線形的に決定されるものであることを示唆していると考えられる．実際，二流体効果を含まない３場モデルによる抵抗性インターチェンジモードのシミュレーションでは，非線形的に磁気島が発生する結果が得られており，このような磁気島の形成はモードの回転に強く影響を及ぼしていると考えられる．そこで本研究では，LHD 実験で観測されているモードの回転メカニズムを明らかにするため
図８ （a）#= 0.08 と（b）#= 0.02 における飽和状態での磁気面．

に，二流体効果を含んだ４場モデルを用いて抵抗性イン
ターチェンジモードの非線形シミュレーションを行った．
""#$平 面 で の ス ラ ブ プ ラ ズ マ を 考 え，""!お よ び#"""に固定境界を課 し，#方 向 は 周 期 境 界 条 件 を課し!
$で磁場強度が０となた．初期平衡状態においては，""!
るような tanh 型の磁場強度分布を仮定した．ただし，電流勾配は０とすることにより，テアリングモードに対しては安定な系を考えている．圧力分布に対しても tanh 型の分布を仮定し，""!にて圧力がもっとも高く，"が増加するにつれて単調減少するような分布を仮定した．また，磁場の曲率に由来する実効的な重力が "方向に存在すると仮定し!
$でのベータ値を 0.5％ とし，抵抗性インターチェ
た．""!
ンジモードが線形不安定である場合を考えた．このとき，線形シミュレーション解析から抵抗性インターチェンジ
モードの実周波数が，電子の流体速度の約６０％であることがわかった．また，高次モードは二流体効果と有限の散逸の効果により安定化されることがわかった．

図９!
$）での磁気フ線形段階においては，磁気中性面（""!
ラックスの摂動成分は０であり，磁気島は生じない．しかしながら，非線形段階においてモードカップリングにより磁気島の幅に対する磁気島の伝播速度の依存性．電子の流体 速 度 及 び E!B フ ロ ー 速 度 も プ ロ ッ ト し て い る．正（負）の速度は電子（イオン）の反磁性ドリフトの方向に流れていることに対応する．

低次モードの磁気島が生ずる．図８では，#"!
!
!%およびは，非線形的に磁気島が形成されモード構造が変化するこ#"!
!
!#の場合に対する非線形飽和状態での磁気面の様子とで説明できることを明らかにした．

を表したものである．ここで #は装置のスケール長に対するイオンのスキン長の比である．#"!
!
!%では２つの磁気２．
１．
６ 磁気島の成長とヒーリングの数値シミュ
レーション

島が，#"!
!
!#では１つの磁気島が形成されていることがLHD における !""磁気島の成長及びヒーリングの機構

わかる．
非線形状態においては，磁気島内部の圧力分布が平坦化とともに磁気島の閉じ込め悪化への影響を，数値シミュされるため，磁気中性面での電子の反磁性ドリフト速度はレーションおよび理論により理解および予測することを目減少する．一方，イオンの反磁性ドリフト方向に E×B フ

指して研究を進めている．特に，磁気島周辺の径電場，温ローが形成されるようになり，電子の流体速度はこれらの度および密度分布に着目し，従来の理論モデルでは難し

速度の和となる．飽和状態での磁気島の幅は，二流体効果

かった磁気島成長およびヒーリングの理解を行った．

のパラメータ #に依存する．#の値を変えた非線形シミュ

第一に，従来の磁気島成長の数値シミュレーションではレーションを行い，得られた飽和状態での磁気島の幅と磁イオン温度を０と仮定していたため反磁性効果による磁気気島の伝播速度の関係を図９に示す．磁気島の幅が小さい

島周りの流れ生成が正しく評価できなかった点を改良しときは，磁気島の伝播速度は電子の流体速度（電子の反磁た．新しい数値シミュレーションでは有限のイオン温度を性ドリフト速度とE!Bフローの和）よりも小さいままであ考慮して数値シミュレーションを行い，有限イオン温度にる．しかしながら，磁気島の幅が大きくなるにつれて，伝起因する新たな磁気島成長の機構を明らかにした．この新播速度は電子の流体速度に近くなる．このように，LHD

たな成長機構は分極電流による磁気島駆動効果である．分

実験でみられる抵抗性インターチェンジモードの回転速度

極電流は磁気島の成長またはヒーリングを決定する中心的７９１

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

役割を果たす．分極電流を評価するためには磁気島と流れて磁気島を成長させることを新たに明らかにした（図１１）との相互作用を理解しなければならない．このときに有限［１０］．そして，イオン温度勾配不安定性乱流により巨視的

イオン温度による反磁性流が決定的な役割を果たす．JIFT

な MHD 現象に対する乱流粘性の値を新たに評価すること

研究者派遣に基づいてテキサス大学と共同研究を行い，磁に成功した（図１２）［１０］．乱流粘性はジャイロボームで 0.2

気島成長に重要な分極電流による磁気島成長効果を明らか（MHD 規格で 10−5）程度である．この乱流粘性値は今後のにした（図１０）［９］．さらに，分極電流効果がヒーリングにMHD シミュレーションにおける粘性の一つの参照値とし寄与する場合であっても，磁気島の背景にイオン温度勾配て非常に有用である．また，高温プラズマにおける磁気島

不安定性駆動乱流がある場合，乱流は分極電流効果を通じ生成に重要となる電子慣性効果による磁気リコネクションの運動論シミュレーションを行い，新たなエネルギー変換
過程を明らかにした［１１］．

２．
１．
７ LHD に お け る ペ レ ッ ト 入 射 シ ミ ュ レ ーション
ペレット入射により生成される低温高密度のプラズモイ
ドの挙動はトカマクと LHD では異なる．これは磁場の結合長が異なるためであることがこれまでの研究で明らかになった．トカマクでは常にプラズモイドが低磁場方向にドリフトするのに対して，LHD においては，プラズモイドの位置によってドリフト方向が変化する［１２］．このメカニズ
ムを解明するために，ペレット溶発 MHD シミュレーショ
ンを行った［１３］．
図１０ 分極電流による規格化磁気島駆動自由エネルギー !pol の磁
気島幅 W 依存 性．こ こ で $s は イ オ ン の ラ ー モ ア 半 径．
!pol ＞ 0 が不安定，!pol ＜ 0 が安定であり，その絶対値が大
きいほど不安定効果または安定化効果が大きい．図は磁気
島の成長およびヒーリングにおける分極電流効果の重要性を示す．冷たいイオン近似を用いた場合 Ti = 0，磁気島の幅によらず分極電流は磁気島をヒーリングする．一方，冷たいイオン近似を用いない場合（Finite Ti），磁気島幅がラーマー半径の５倍以上であると磁気島が不安定になる．文献
［９］
より転載．

LHD にはトロイディシティによる磁場強度とヘリシティによる磁場強度がある．図１３（a）と（b）は横長ポロイダ
ル断面と縦長ポロイダル断面を示しており，それぞれにプ
ラズモイド１，２が位置している．プラズモイド１はトロ
イディシティ高磁場側，かつヘリシティ低磁場側に位置しており，プラズモイド２はトロイディシティ高磁場側，かつヘリシティ高磁場側に位置している．これらのプラズモ
イドの磁気面 平 均 密 度 の 時 間 発 展 を 図１４（a）と（b）に示

図１
１ イオン温度勾配乱流に影響された磁気島の成長を調べた．
乱流による規格化磁気島駆動自由エネルギー !の粘性 #依存性．!＞0 は乱流が作る分極電流が磁気島を駆動すること，!＜0 は乱流が磁気島をヒーリングすることを示す．そして，その絶対値が大きいほど不安定効果または安定化効
果が大きい．規格化温度勾配 "i が０の場合，イオン温度勾
配不安定性は安定で，この場合，!＜0 で磁気島は安定な状況である．"i ＝ 2.5 の場合，イオン温度勾配不安定性が不安
定で，!＞0 となり磁気島はイオン温度勾配乱流によって駆動される．また，磁気面平均磁場局率κが有限な場合は，イオン温度勾配不安定性がさらに不安定となり，乱流は分極電流効果により磁気島を強く駆動する．文献［１
０］
より転
載．

図１
２ 磁気島が乱流粘性を通じて流れから受ける抵抗−df /du を計算し，イオン温度勾配不安定性乱流による乱流粘性の値を評価した．規格化温度勾配 "i が０の場合，イオン温度勾
配不安定性は安定で流れは層流となり，抵抗は粘性 #に比例する．"i = 2.5 の場合，イオン温度勾配不安定性が不安定で，乱流による抵抗が磁気島に作用し，その抵抗値は粘性#にほとんど依存しない．図の黒線と青線を #が大きいほうに外挿し赤線に交わる点の粘性 #は 0.2程度である．このことは磁気面平均磁場局率κが有限であっても０であっても乱流粘性はジャイロボームで 0.2（MHD 規格で 10−5）程度であることを示す．文献［１
０］
より転載．

７９２

Project Review

2.1 MHD Equilibrium and Stability

K. Ichiguchi et al.

図１
３ （a）
横長ポロイダル断面におけるt = 20 μsのプラズモイド１の密度分布．（b）
縦長ポロイダル断面における t = 20 μs のプラズモイド２の密度分布．密度は背景密度によって規格化されている．

す．横軸は規格化小半径である．プラズモイド１は小半径方向（大半径負方向），つまり低磁場方向にドリフトしているが，プラズモイド２は小半径負方向にドリフトした
後，正方向（大半径負方向），つまり高磁場方向にドリフト
図１
４ （a）
Boozer 座標におけるプラズモイドの磁気面平均密度の時間発展．横軸は規格化小半径，縦軸は背景密度による規格化密度，（a）（b）
はぞれぞれプラズモイド１，
２を示す．

している．このことは，プラズモイドの挙動は初期位置によって異なることを意味する．プラズモイドには，主に1/R による力と磁場の復元力による力が働く．結合長が長いトカマクにおいては前者が相対的に大きくなり，常に大半径方向へドリフトする．一方で，プラズモイド１，２に

参考文献

おいては結合長が短く，後者がプラズモイドの運動を支配［１］S. Sakakibara et al., Proc. 23rd Fusion Energy Conf. Oct.1116, 2010, Daejeon, EXS/P5-13.
［２］K. Ichiguchi et al., Nucl. Fusion 43, 1101 (2003).
［３］K. Ichiguchi et al., Nucl. Fusion 55, 043019 (2015).
［４］S. Sakakibara et al., Nucl. Fusion 53, 043010 (2013).
［５］K. Ichiguchi et al., Nucl. Fusion 55, 073023 (2015).
［６］J. Miyazawa et al., Nucl. Fusion 54, 043010 (2014).
［７］T. Nicolas, and K. Ichiguchi, Nucl. Fusion 56, 026008 (2016).
［８］T. Nicolas et al., Plasma Fusion Res. 10, 3403018 (2015).
［９］A. Ishizawa et al., Phys. Plasmas 19, 072312 (2012).
［１０］A. Ishizawa and F.L. Waelbroeck, Phys. Plasmas 20,
122301 (2013).
［１１］A. Ishizawa and T.-H. Watanabe, Phys. Plasmas 20,
102116 (2013).
［１２］R. Sakamoto et al., Nucl. Fusion 44, 624 (2004).
［１３］R. Ishizaki and N. Nakajima, Plasma Fusion Res. 9,3403130 (1〜5) (2014).

する．後者，つまりこの磁場の復元力は，プラズモイドの位置における低磁場方向の向きに加えて，磁力線に沿った
磁場強度分布によって決まる．プラズモイド１においては，低磁場方向は大半径負方向であり，磁力線に沿った磁場強度は極小になっている．溶発により磁力線方向に膨張するプラズモイドは，磁力線に沿って磁場強度が大きくなるため，その運動が妨げられる．これが結果として大半径負方向へのドリフトを生じさせる．一方で，プラズモイド２においては，低磁場方向は大半径正方向であり，磁力線に沿った磁場強度は極大になっている．つまりドリフト方向を決定する磁場の二つの効果がプラズモイド１とは全く逆であるため，逆の逆は同じとなりプラズモイド１と同じ大半径負方向にドリフトする．これが高磁場方向にドリフトするメカニズムである．今後は，実験でのペレット入射に対応したモデルに拡張し，RMP の影響などを評価していく．

７９３

J. Plasma Fusion Res. Vol.92, No.11 (2016)７
９４‐８０５

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

２．核融合プラズマシミュレーション研究
2. Simulation Researches in Fusion Plasmas

２．
２ コアプラズマの輸送
2.2 Transport Simulation of Core Plasmas
菅 野 龍 太 郎１，２）， 佐 竹 真 介１，２）， 沼 波 政 倫１，２）， 仲 田 資 季１，２）， 登 田 慎 一 郎１），
山 岸 統１）， 洲 鎌 英 雄１，２）， 渡 邉 智 彦３）， 石 澤 明 宏４）
KANNO Ryutaro１，２）, SATAKE Shinsuke１，２）, NUNAMI Masanori１，２）,
NAKATA Motoki１，２）, TODA Shin-ichiro１）, YAMAGISHI Osamu１）,
SUGAMA Hideo１，２）, WATANABE Tomo-Hiko３）and ISHIZAWA Akihiro４）
１）

自然科学研究機構核融合科学研究所，２）総合研究大学院大学，３）名古屋大学，４）京都大学（原稿受付：２
０
１
６年７月１
３日）

数値実験炉研究プロジェクト新古典・乱流輸送シミュレーショングループにおけるコアプラズマの輸送研究に関する最新の成果について報告する．新古典・乱流輸送シミュレーショングループでは，ヘリカルプラズマをはじめとするトーラスプラズマの密度分布，温度分布，閉じ込め時間の予測に有用な知見を得るため，コア領域における輸送現象に対して，ドリフト運動論およびジャイロ運動論に基づく輸送モデルや大規模計算機シミュ
レーションコードを開発し，また実験結果との比較によるコードの検証を行っている．新古典輸送の局所近似理論モデルの構築，磁場摂動の新古典粘性や熱拡散への影響の評価，装置間・コード間比較検証について報告する．
また，ジャイロ運動論シミュレーションに基づく，大型ヘリカル装置（LHD）や JT-60U における乱流輸送の再現と予測，輸送モデルの構築やシミュレーションコードの拡張，新しいゾーナルフロークロージャーモデルに基づくジャイロ流体シミュレーションなどの成果を紹介する．
Keywords:
drift-kinetic simulation, gyro-kinetic simulation, gyro-fluid simulation, integrated simulation, neoclassical transport,
turbulent transport

２．
２．
１ はじめに比較による予測性能の検証が複数行われ，開発した輸送シ

磁場閉じ込めトーラスプラズマ中の輸送現象の物理機構

ミュレーションコードの実証研究も進展している．次節か解明や輸送特性の定量評価は将来の核融合炉における炉心らは，新古典輸送に関しては，（１）局所近似理論モデルの性能の評価や予測・最適化にとって不可欠であり，実験解構築，（２）新古典粘性のプラズマ回転に与える影響の実験

析に加えて理論・シミュレーション研究が重要な役割を結果との比較，（３）装置間・コード間国際比較研究，およ担っている．本章では，数値実験炉研究プロジェクト新古

び（４）新古典輸送現象への共鳴摂動磁場の影響について報典・乱流輸送シミュレーショングループにおけるコアプラ

告する．また，乱流輸送に関しては，（５）LHD 有限ベータズマの輸送研究に関する平成２５年度から２７年度の期間におプラズマの乱流輸送に対する理解と予測の進捗状況，（６）

いて主に得られた成果について報告する．新古典・乱流輸

JT-60U トカマク実平衡における ITG-TEM 乱流の解析結送シミュレーショングループでは，ヘリカルプラズマをは

果，（７）ジャイロ運動論解析に基づく乱流輸送モデルの統じめとする磁場閉じ込めトーラスプラズマの密度分布，温合コードへの適用，（８）ジャイロ運動論的シミュレーショ

度分布，閉じ込め時間の予測に有用な知見を得るため，コンのためのフラックス・チューブ連結モデルの開発，（９）

ア領域における新古典および乱流輸送現象に対して，ドリ燃焼プラズマにおける乱流輸送解析に向けた多粒子種モデ

フト運動論およびジャイロ運動論に基づく理論モデルや計ルの開発，（１０）国際ベンチマークの推進，（１１）乱流計測シ

算機シミュレーションコードの開発や拡張を進め，将来のミュレータによる数値乱流計測，および（１２）新しいゾーナ

燃焼プラズマを見据えた大規模な理論・シミュレーション

ルフロークロージャーモデルを用いたジャイロ流体 ITG研究を展開している．また，LHD や JT-60U をはじめとす

乱流シミュレーションについて成果を紹介する．

るヘリカル型およびトカマク型実験における計測結果とのcorresponding author’s e-mail: kanno@nifs.ac.jp

７９４!2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchProject Review

2.2 Transport Simulation of Core Plasmas

R. Kanno et al.

"""!における粒子フラックスの極端なピークが５次元グ

２．
２．
２ 新古典輸送計算における磁気ドリフト項の重要性

ローバル計算においては抑制され，またピークの幅が拡

ヘリカルプラズマにおける新古典輸送は，３次元磁場配がって全体に ""!!側にシフトしていることがわかる．ま位中における荷電粒子のガイディングセンタードリフト運

た，磁場ドリフトの磁気面接線成分を残した ZOW モデル

動と，特にヘリカル磁場リップルに捕捉された粒子のクーだけが，グローバル計算に見られるこの傾向をよく再現す

ロン散乱によって決まるが，ヘリカル装置が作る磁場配位ることが判った．ヘリカル磁場リップルに捕捉された捕捉の特徴や，衝突周波数領域，両極性径電場による "!!

粒子のポロイダル方向の歳差ドリフトが，磁場ドリフト項

ドリフト速度の大小など，様々な因子に強く依存する．そを完全に無視してしまうと """!において完全にゼロとなのため，ガイディングセンター運動に対して様々な近似をり（ポロイダル共鳴），全ての粒子が大きな径方向輸送に行ってこれを解く従来の数値解法には近似に起因する誤差寄与することになってしまう．一方，磁場ドリフトの磁気が含まれるが，それがどの程度定量的な精度に影響するか

面接線方向成分の有限性を残すと，このポロイダル共鳴はについては，これまでほとんど議論されてこなかった．そ"!! ドリフトと磁場ドリフトのポロイダル成分の釣り合

うした近似法として広く用いられているものには，ガイいによって決まるため，イオンの場合 ""!!側にシフト

ディングセンター軌道の小半径方向ドリフトを高次の微小し，かつ個々の粒子の速度によって釣り合いの条件が変わ

量として無視する局所近似，荷電粒子の運動エネルギーのるために，共鳴による径方向フラックスが低く抑えられる

変化を無視した mono-energy 近似などがある．特に従来ののである．こうして，ヘリカル系の新古典輸送計算におい局所近似では，磁場の勾配ドリフト及び曲率ドリフト（磁て磁場ドリフト項の磁気面接線方向成分が従来の局所近似

場ドリフト項）を小半径方向だけでなく磁気面接線成分もとグローバル計算の差の主な要因であることが明らかと含め全部落とす近似が一般的に使われてきた．このようなった．

に，本来５次元位相空間におけるドリフト運動論方程式にエネルギーの２つの次元をパラメータ化することで３次元

２．
２．
３ 新古典粘性がプラズマ回転分布に与える影響に簡約化した局所近似・mono-energy 近似の新古典輸送計

新古典輸送理論は径方向輸送やブートストラップ電流を算法（例えば DKES コード，GSRAKE コードなど）と，そ

扱うだけではなく，トーラス磁場の磁気面上での磁場強度のような近似を用いず５次元ドリフト運動論方程式を直接の不均一性に起因する新古典粘性の評価にも用いられる．

解くグローバル新古典輸送計算コード（FORTEC-3D）の新古典粘性はプラズマのトロイダル・ポロイダル回転分布

計算結果には，特に低衝突領域で "!! ドリフト速度が大の決定に重要な役割を果たすもので，最近のトカマク研究

きくない場合に大きな差が生じることがこれまでの研究でにおいても注目されている物理の一つである．共鳴摂動磁

わかってきた．しかしながら，両者によって求められる新場（RMP）印加実験やトロイダルコイルが作るリップル磁

よって扱われる新古典輸送現象を，小半径方向座標と運動

古典粒子・熱フラックスには様々な近似に基づく差が混在化しており，具体的にどのような因子が最も大きく差異を生じさせているのか特定しづらかった．
そこで，両モデルの間をつなぐモデルとして，磁場ドリフト項の磁気面垂直成分だけを高次の微小量として無視
し，磁気面接線成分を残し，かつ運動エネルギーの次元を残すことで結果的に DKES で用いられているような非圧縮"!! フローの近似を使わなくてよい，新しい４次元簡約化ドリフト運動論方程式と，それを FORTEC-3D をベース
とした #
［１］．
#モンテカルロ法で解く数値手法を開発した
この，磁気面垂直成分のドリフト項のみを落としたモデル
を Zero-Orbit-Width（ZOW）近似と呼び，接線成分も含めすべての磁場ドリフト項を落としてしまう従来の ZeroMagnetic-Drift（ZMD）近 似 と 区 別 す る．こ う し て，
FORTEC-3D コードは５次元のグローバルモデルから，４
図１次元の ZOW または ZMD 近似，さらに mono-energy 近似と非 圧 縮 "!! フ ロ ー の 近 似 ま で 取 り 入 れ た３次 元 のDKES 極限と，様々な段階の近似を用いた新古典輸送計算
を１つのコードで行うことが可能となり，それらの比較も統一的に行うことができるようになった．
様々なドリフト運動論モデルに基づいて LHD 配位におけるイオンの新古典粒子フラックスの径電場依存性を比較したものが図１である．ZMD や DKES モデルに見られる，
７９５

様々な近似モデルのドリフト運動論方程式による新古典径方向粒子束の径電場依存性の比較．ZOW，ZMD，DKESlim．はそれぞ れ Zero-orbit-width，Zero-magnetic-drift，
DKES 極限（ZMD＋非圧縮 E ×B フロー）の局所近似計算法
を 表 す．GLOBAL は５次 元 FORTEC-3D コ ー ド の 結 果，
GSRAKE は局所近似＋バウンス平均の計算法である．磁場ドリフト項を考慮しているGLOBALとZOWモデルだけが，Er = 0 における極端な "i のピークが緩和されており，磁場ドリフト項の有無が新古典輸送計算に大きく影響することが示されている．

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

場成分など，実際のトカマク磁場に存在する微弱な非軸対コードによる新古典粘性のシミュレーション結果と，バイ

称磁場摂動による新古典トロイダル粘性が，プラズマのト

アス電圧‐バイアス電流関係の観測から推定された新古典

ロイダル回転分布に与える影響を予測するためには高精度ポロイダル粘性がピーク値をとるバイアス電圧の値の間にな新古典トロイダル粘性の計算法が必要とされる．

有意な一致を示すことに成功した．

そこで，FORTEC-3D コードを用いて新古典ポロイダ

２．
２．
４ 新古典輸送計算の装置間・コード間国際比較検証研究ル・トロイダル粘性を評価できるようにコードの拡張が行われた．新古典粘性は圧力テンソルの磁力線平行・垂直方
向成分の非等方性によって生じるので，ドリフト運動論に２．
２．
２節で述べたように，新古典輸送計算法における近従っ て 分 布 関 数 を 計 算 す る FORTEC-3D コ ー ド に お い

似，特に従来の局所近似ドリフト運動論方程式モデルにおて，これを直接的に数値的に評価する手法を開発した．そける磁場ドリフト項の磁気面接線成分を無視することが，して，そのシミュレーションの応用研究として，日本原子新古典径方向粒子・熱輸送評価の定量的信頼性に対する影力研究開発機構（原研機構，現 QST）との共同研究におい響が大きいことがわかってきた．こうした局所近似計算て，JT-60U 装置のトロイダル磁場コイルが作るリップルと，FORTEC-3D コードのように近似を用いないグローバ

磁場成分による新古典トロイダル粘性を FORTEC-3D コー

ルコードとの計算結果の差が，実際のヘリカル装置におけドで計算し，その結果を原研機構で開発された統合輸送

る両極性径電場の予測や，両極性条件における粒子・熱フ

コード TOPICS に渡すことで，従来の輸送シミュレーショ

ラックスの計算結果にどの程度見られ，それが乱流輸送レンで無視されていた新古典トロイダル粘性を含むプラズマ

ベルと比べて無視できるのか否か，という問いに答えるた

回転分布の予測計算を行うことを可能とした．その結果，

めにまず，新古典輸送コード間のベンチマークによる各

図２に示すように新古典トロイダル粘性の影響を含めた計

コードの実装が正しいかの検証（Verification）と，新古典算結果が，それ以前に行われていた NBI トルク入力と高速輸送計算が現実のプラズマの輸送現象をどの程度再現できイオンの軌道損失による #!! トルクのみを考慮していたているかの検証（Validation）を様々なヘリカル装置で行TOPICS コードによるシミュレーション結果に比べ，観測うことでシミュレーション自体の信頼性を証明する必要が結果により近い回転分布を再現できることを示した［２］．

あると考えた．

今後，さらにトカマクや LHD における RMP 印加実験でのそこで核融合科学研究所の新古典輸送の研究グループが新古典粘性計算と，プラズマ回転分布への影響に関する研主導して，欧米のヘリカル装置を持つ研究機関の研究者に究を進める予定である．

呼びかけて新古典輸送計算の装置間・コード間国際比較検

その他の新古典粘性計算の応用例として，LHD における証研究をこの数年間精力的に進めてきた．この共同研

バイアス実験において，バイアス電場によって "!! 回転

究には，スペイン CIEMAT，ドイツ Max-Planck 研，米国を外部から駆動した場合の，新古典粘性のポロイダルフ

Wisconsin 大 か ら の 協 力 が 得 ら れ た．そ し て，LHD やロー速度に対する非線形な応答に関する数値シミュレー

TJ-II，W7-AS，HSX 装置において，新古典輸送シミュレー

ションと実験観測の比較研究が行われた［３］．FORTEC-3D

ションの検証が行えるように，電子・イオン温度分布，密度分布，並びに HIBP や CXRS などによる "!! 回転計測
を 揃 え た 実 験 デ ー タ を 取 得 し，FORTEC-3D やGSRAKE，DKES コードで新古典輸送計算を行った．またそのシミュレーション結果を，加熱分布と温度分布の時間変化から評価したパワーバランス解析による熱輸送と比較した．新古典輸送シミュレーションから得られた両極性径
電場とプラズマ回転観測の比較も行った［４］．
新古典輸送計算の Validation の難しさの一つは，実際のプラズマ中の粒子・熱輸送が新古典輸送のみで決まっているわけではない点が挙げられる．ただし，微視的乱流輸送
はヘリカルプラズマの両極性径電場の決定にはあまり関与していないとみなされており，まず径電場分布を実験とシ
ミュレーションで 比 較 することが Validation の一つとなる．
LHD における，$""$!のイオンルート及び $"!$!の 電 子ルートが想定される実験の解析例を図３に示す．３つのシ

図２

ミュレーションで求めた両極性径電場分布は観測結果とあJT-60U の NBI Co，Counter，Balance 入射実験のそれぞれにおいて，FORTEC-3D で計算した新古典トロイダル粘性（NTV）の影響を含める・含めないでトロイダル回転分布の計算を TOPICS コードで行った結果（実線及び破線）と，ト
ロイダルフローの計測結果（点）の比較．NTV を考慮した計算の方がトロイダル回転の観測結果をよく再現している．

る程度よい傾向の一致を見た．ただし，FORTEC-3D コー
ドと GSRAKE，DKES コードでは "#分布に有意な差が出ている．これは，２．
２．
２節で見た磁場ドリフト項の有限性による新古典粒子フラックスの径電場依存性によるものである．また，図３（a）のショットでの電子・イオン熱フラック
７９６

Project Review

2.2 Transport Simulation of Core Plasmas

R. Kanno et al.

図３

LHD 実験における，（a）
Ti ! Te，（b）
Ti ＜ Te のプラズマにおける両極性径電場のシミュレーション結果と観測結果の比較．なお，（b）
において局所近似コード（DKES，GSRAKE）は３根の両極性解を示しているが，グローバル FORTEC-3D は径電場分布を時間発展で解くため，電子ルート（Er ＞ 0）のみが解として得られている．HIBP はコア側，CXRS はエッジ側において観測精度が高い．

図４

LHD 実験の Ti ! Te プラズマにおける３つの新古典シミュレーションによる解析によるトータル熱フラックスの評価値（Exp.）の比較．

新古典イオン・電子熱フラックスとパワーバランス

スの計算結果とパワーバランス解析の比較結果を図４に示の少ない FORTEC-3D コードに置き換 え る こ と に よ っ

す．このショットではイオン新古典熱フラックス !"がトーて，様々な誤差が入りやすいコアプラズマの輸送解析におタルの約半分ほどを占めているが，電子熱フラックスにおいて少しでも不確定要因を減らすことが，特に新古典輸送

ける新古典輸送の割合は低かった．乱流熱輸送レベルはが全体に占める割合が大きいケースになるほど重要とな

Q（Exp）‐Q（Neoclassical）によって評価されるが，新古典る．

輸送の占める割合が大きい場合，新古典輸送計算自体の近２．
２．
５ 新古典輸送現象への共鳴摂動磁場の影響

似による誤差が乱流輸送の評価にも影響を与えることになる．この例では３つの新古典コードの計算結果には大きな共鳴摂動磁場（RMP）を静的に加えた領域（例えば，プ

差がなかったが，図５に示した TJ-II のイオンルートプラ

ラズマエッジのエルゴディックな磁力線領域）におけるプズマの実験解析結果を見ると，FORTEC-3D と DKES コー

ラズマの輸送現象の物理的理解は，LHD プラズマのエッジ

ドでは約２倍の差が !"，!!に生じており，また新古典輸送における粒子・熱閉じ込めや，ダイバータへの熱負荷を制がトータルの熱輸送に占める割合もかなり高い．このよう御するためITERなどで計画されているようなRMPを用いな状況は LHD でも特にイオンルートプラズマにおいて見た Edge-Localized Mode（ELM）の抑制を行った際の閉じられる傾向であり，高い新古典輸送計算の精度が要求され込めの性能評価に重要である．近年の実験研究，例えば，

る．実際，図５（b）において DKES コードの !!はパワーバ

DIII-D トカマクにおける RMP を用いた ELM 抑制実験の結ランスによるトータル熱輸送の評価値を超えてしまってお果，低衝突なプラズマエッジ（衝突頻度 !!! 0.1）におけるり，この結果を信じるなら，パワーバランス解析による熱

粒 子・熱 輸 送 は，従 来 の 標 準 的 な 理 論 モ デ ル で あ る

輸送が過少評価されたか，乱流熱輸送が内向きになっていRechester-Rosenbluth の磁力線拡散理論（以下，R-R 理論

ると解釈することになってしまう．新古典輸送計算を近似と呼ぶ）に基づいた予想とは，かけ離れていることが明ら
７９７

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

図５

TJ-II プラズマにおける新古典イオン・電子熱フラックスとパワーバランス解析によるトータル熱フラックスの評価値（青線）の比較．FORTEC-3D の計算結果が赤線．DKES の結果（緑）には温度密度分布計測の誤差伝搬によるエラーバーが評価されている．

かになった．ここで，ELM 抑制実験については，文献［５］を参照．この実験事実により，このようなプラズマエッジ
における輸送の基本的性質に対する物理的理解の見直しが早急に必要となっている．R-R 理論と同様に電場を無視する仮定の下で実行した "ドリフト運動論的シミュレーショ
!
ン（KEATS コード）によって，（１）案内中心分布の準定常状態における径方向のイオン熱拡散係数は，クーロン散乱
効果によって R-R 理論値と比べて極めて小さくなることを確認し，（２）その物理パラメータ依存性については，R-R
理論モデルとほぼ同じであることがわかった．また，（３）
RMP の効果は，衝突頻度が大きくなるに従い弱められることが明らかになった．図６は，イオン熱拡散係数に対す
る RMP の影響について，まとめたものである．詳細は，文献［６］を参照．実験と定性的に矛盾しない結果が得られた
ため，定量的な解析ができるよう，電子およびイオンの粒
子フラックスに対する両極性条件から決まる径電場の効果も含めて輸送解析ができるようにコード開発を行った．初期的な計算結果として，正の径電場によって電子の粒子フ

図６

ラックスは抑制されるということを見出した．以上の初期成果により，両極性条件を満たす径電場を考慮して衝突輸送に対する RMP の影響を解析する準備が整った．

２．
２．
６ LHD 有限ベータプラズマの乱流輸送の理
解と予測

円形トカマク磁場に RMP を加え，磁力線構造がエルゴダイ
ズした領域において，RMP が径方向のイオン熱拡散係数＊が大きくなるに従い弱めら!r に与える影響は，衝突頻度 #れる．図中の●は新古典熱拡散係数，■は RMP を加えた場
合の熱拡散係数である．この計算例では，低衝突領域にお＊! 0.04 であっても），RMP を加えた場合の熱拡いても（#散係数は，新古典熱拡散係数に対しておよそ１０倍大きくな
る程度であった［６］．

ヘリカルプラズマのジャイロ運動論シミュレーション研て定量評価可能となった．また，磁場揺動の導入により

究は磁場配位が複雑な３次元構造をもつために計算科学のMHD 不安定性に対する最適化と乱流輸送に対する最適化観点からも挑戦的な課題である．そして核融合炉をめざすの関係も理解できる展望が開けた．したがって，ヘリカル上で有限ベータ乱流輸送の理解および予測は不可欠であるプラズマにおける乱流輸送の理解と予測に対する非常に大にも関わらず，これを解析できるジャイロ運動論シミュきな進展が期待される．ヘリカル・ステラレータ系の有限レーションコードは開発されていなかった．我々は有限ベータ乱流輸送に対する非線形ジャイロ運動論シミュレー

ベータプラズマに適用可能な電磁的ジャイロ運動論シミュ

ションは他に例がないほど困難である．またトカマクに対レーションコードを開発した［７］．その結果，これまでプラして行われたシミュレーションでは，ゾーナル流が弱くな

ズマベータが０の極限に限られていたヘリカルプラズマのり，その結果，飽和が得られない又は輸送レベルが極めてジャイロ運動論シミュレーションの適用範囲が有限ベータ

高くなることが指摘されている．このように，ほとんど未領域へと広がると共に電子・イオンの粒子・熱輸送がすべ知の，有限ベータ乱流輸送に対する非線形ジャイロ運動論
７９８

Project Review

2.2 Transport Simulation of Core Plasmas

R. Kanno et al.

シミュレーションを進め，乱流の飽和機構を調べた．そのいる．一方，図８（d）に示すように，横軸に示した温度勾配

結果，ヘリカルプラズマ特有の乱流抑制機構を明らかにしが２０％程度下がるだけで，縦軸にある熱輸送（計算結果をた（図７）．この抑制機構は非軸対称性磁場に捕捉された実験値で規格化したもの）は 0.4 倍にも低下することが，同粒子に起因し，高ベータで不安定になる運動論的バルーニ

時に明らかとなった．すなわち，熱輸送はイオン温度勾配

ングモードの空間構造がトカマクと比較して大きく異なるの値に非常に敏感であり，実験における熱輸送を再現する

点にある［７］．ヘリカルプラズマの場合，ゾーナル流が弱温度勾配の評価は２０％程度の誤差であることが定量的に示くとも，不安定なバルーニングモード同士がこのモード構された．実験で観測されるプラズマは，上述のイオン温度

造を持つがゆえにお互いに引きちぎりあって飽和する［８］勾配による不安定性を引き起こす臨界勾配状態に近い密度（図７）．この飽和機構は，"!ジャイロ運動論方程式が持つ及び温度分布を持つと考えられる．さらにプラズマが乱流

２次の保存量［９］を利用して，エントロピー密度が不安定に遷移すると，不安定性理論の予測を超えて，輸送が生じモードから高波数へ非線形効果によって移送されることをる臨界状態は変化することがある．特に，プラズマ乱流が示すことによって明らかになった．この新たな乱流飽和機

作り出すゾーナル流と乱流渦の非線形相互作用に大きく影

構を含めトカマクおよびヘリカルプラズマにおける電磁的響されることが，これまでの理論・数値シミュレーション

乱流のジャイロ運動論シミュレーション研究をレビュー論

によって指摘されてきた．しかし，今回のシミュレーショ

文としてまとめた［１０］．

ンが示すところは，この従来のパラダイムからの帰結とはLHD 実験で得られた磁場形状やプラズマ密度と温度を

必ずしも一致せず，結果として，LHD プラズマにおいて乱計算モデルに忠実に反映して乱流輸送のシミュレーション

流輸送を駆動する臨界温度勾配は，線形理論が与える評価を行う場合，有限ベータ効果よりも運動論的電子効果が大に近いことがわかった．これは LHD プラズマではトカマきな影響を与える．運動論的電子効果の導入はイオン温度

ク型プラズマと比較して，乱流による流れの自己組織化効

勾配不安定性の成長率を増大させるだけでなく従来の電子果が弱いことを示唆しており，高性能ヘリカル型プラズマの断熱応答を仮定した解析ではできなかった乱流による粒の乱流輸送モデル構築に向けた重要な結果である［１１］．

子輸送および電子熱輸送の評価を可能にする．今回の計算
験炉において有望視されている高性能なヘリカル型プラズ

２．
２．
７ JT-60U トカマク実平衡における ITG-TEM
乱流輸送解析

マの一つで，この放電におけるプラズマはイオン温度勾配

ITER や DEMO といった次期核融合実験炉における乱流不安定性に駆動される乱流が輸送を引き起こすことが知ら

輸送特性の定量評価及び予測には，ジャイロ運動論に基づれている．運動論的電子効果を含んだ計算結果を図８に示いた第一原理乱流シミュレーション解析・予測が有効であす．図中で Time＝４０以降は乱れを含んだ準定常状態が得

り，燃焼プラズマ中の乱流輸送解析に向けてコード開発やられており，この状態が実験で観測される状態に対応す

V&V が進められている．本研究では，ジャイロ運動論的乱る．図８（c）は粒子がトーラス中心方向に輸送されている

流シミュレーションコード GKV に実プラズマ形状と平衡

対象は LHD における高イオン温度放電と呼ばれ将来の実ことを示す．図８（b）は電子のエネルギー輸送の時間発展分布データを導入するインターフェースを新たに構築した（赤線）を示す．黒色実線は実験観測値を示し，両者はよく［１２］．これにより，JT-60U 実験プラズマに対する初めて一致していることがわかる．図８（a）はイオン熱輸送の時の第一原理乱流シミュレーション解析が実現され，微視的間発展を示す．赤線がシミュレーション結果，黒線が実験

不安定性特性や乱流スペクトル構造の同定に加え，イオン

観測値を示す．シミュレーションによるイオンの熱流束

及び電子に関する熱輸送レベルの比較が行われ，GKVによは，このケースでは実験観測値のおよそ３倍の値を与えてる乱流輸送シミュレーションの定量予測性能の検証が成さ

図７

大型ヘリカル装置における非軸対称トーラスプラズマ内の乱流の静電ポテンシャル揺動分布を示す．左図は低ベータで乱流輸送を支．右図は高ベータで乱流輸送を支配する運動論的バルーニングモード乱流を配するイオン温度勾配不安定性乱流を示す（!＝ 0.3％）
示す（!＝ 1.7％）．高ベータでは乱流の空間構造が低ベータと異なり，この構造がヘリカルプラズマ特有の乱流の飽和機構を生み出す．

７９９

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

図８

LHD実験解析結果：（a）
イオンエネルギー流束，（b）
電子エネルギー流束，（c）
粒子流束，（d）
エネルギー流束と温度勾配長の関係図．

れた．ITER や DEMO における規格化プラズマサイズは!!
"#!"!!（!：プラズマ半径，"#：イオンジャイロ半径）を示すため，そのような将来の大型装置に対する予測性能"#""!!）をを評価する観点から，JT-60U L-mode 放電（!!
対象に乱流シミュレーション解析を行った．イオン種に加え，通過／捕捉電子もジャイロ運動論に従う無撞着な取り扱いを実現することにより，半径位置に応じた異なる微視的不安定性特性を同定し，プラズマ内側領域ではイオン温度勾配（ITG）
不安定性が支配的となり，外側領域に向かって捕捉電子不安定性（TEM）へと遷移していく様相を明らかにした［１３］．さらに，それらの微視的不安定性から駆動されるイオンスケール乱流中の熱・粒子輸送レベルを評価し，実験計測結果と比較を行った［１４］．GKV による乱流輸
送シミュレー シ ョ ン は，従 来 の 準 線 形 流 体 輸 送 モデル（TGLF）
では定量性を欠く領域に対しても，イオン及び電子熱輸送レベルに関する実験結果を比較的良く再現することに成功している（図９）．これらの知見はシミュレーショ
ンモデルのさらなる高度化や，建設中の次期実験装置に対する予測研究の基盤構築に資する研究成果となっている．

２．
２．
８ ジャイロ運動論解析に基づく乱流輸送モデルの輸送シミュレーションへの適用方法
図９ （a）
イオン及び（b）
電子熱輸送レベルに関する JT-60U 実
験ー乱流シミュレーション比較．赤線はGKV，青線はTGLF
モデルによる計算結果を示している［１４］．

ジャイロ運動論方程式を解く GKV コードを用いて，ヘリカルプラズマにおける乱流輸送について，LHD でのイオン温度勾配不安定性（ITG）モードや帯状流について研究する．イオン熱拡散係数はジャイロボーム熱拡散係数とあジャイロ運動論解析から計算される混合長概算と帯状流崩る関数の積とす る．あ る 関 数 は ITG 成 長 率 の 混 合長概

壊時間の関数が簡約化モデルとして導出されている．統合

算，磁気シアパラメータ，帯状流の崩壊時間（the decay

コード TASK3D では，熱乱流拡散係数についてはジャイtime）に依存する．イオンの熱拡散係数について，線形ロボームスケーリングが用いられ，LHD における実験結果
８００

Project Review

2.2 Transport Simulation of Core Plasmas

R. Kanno et al.

と解析結果との比較が行われている．輸送コードの中で時行われてきた．プラズマサイズに比べ十分小さなジャイロ

間ステップごとに，ジャイロ運動論線形解析を行うことは半径で特徴付けられる乱流揺動のシミュレーションには，非常に計算コストがかかる．輸送障壁のような改善モード

磁力線に沿った微小領域を扱うフラックス・チューブ・モ

を得ることができる実験シナリオを提案するには，ある程デル［１６］が有効であり，これを用いてジャイロ運動論的シ度の領域でのプラズマパラメータサーベイが必要である．

ミュレーション・コード GKV が開発された［１７］．本研究

輸送コードで乱流輸送係数にどのモデルを選ぶべきか決めでは，従来のフラックス・チューブ・モデルを拡張した新るために，多くの微視的不安定性からどのモードが不安定

し い モ デ ル，フ ラ ッ ク ス・チ ュ ー ブ 連 結（flux tube化しているか，そしてその乱流レベルを評価する必要があtrain）モデル，を開発した［１８］．新しいモデルは，磁力線る．最初のステップとして，LHD での高イオン温度放電%以上）をもつ乱流方向に長い相関長（ポロイダル角 &で #（ショットナンバー８８３４３）における ITG モードについて研揺動を扱う場合に有効であり，時間ステップ幅の節約およ究を行った．混合長概算の値を大半径のイオン温度勾配長びバルーニング表示に関連した対称性の保存において優位に対する比に比例する形でモデル化する．径方向に安全係性をもつ．

数，密度，温度が変わる．したがって，径方向に比例係数

フラックス・チューブ・モデルでは，磁力線に沿って座や，ITG モードが不安定化するイオン温度勾配長を求め標の一つ（これをポロイダル角 &で測る）を取り，磁力線

る．また帯状流崩壊時間の径方向依存性を調べる．このとを横切る方向（トーラス小半径方向および磁力線ラベル方き，磁場配位を時間的に固定する．さらなるモデル化に向）には Fourier モード展開を行う．これはいわゆるバ

よって，高速に簡約化モデルの値を十分な精度で再現する

ルーニング表示&ことができた．さらなるモデル化したイオン乱流輸送係数''!& (
!'(
'!
&)
'!
&"#
%()
$ ! !

を 用 い て，イ オ ン 温 度 の ダ イ ナ ミ ク ス を 輸 送 コ ー ド

!

($!&

TASK3D に よ り 解 析 す る．イ オ ン の 加 熱 入 力 分 布 は，')(')&!&!"#%()#$&%!#(.
（１）

TASK3D から求める．初期条件として，LHDでの高イオン
温度放電（ショットナンバー８８３４３）で定常と考えられるイと密接に関連している．従来のモデルでは，&方向のシオン温度分布を用いる．解析の結果，多重解から正の電場（"& は正の整数）
ミュレーション領域を，!"&%%&%"&%

を選んだ時に，実験結果と矛盾のない結果が得ることがでとし，磁力線平行方向への相関長が長い場合には "& を大きた（図１０）．シミュレーションの際，ハッチした領域ではきくしてシミュレーション領域を拡大し，計算結果の収束

実験値を境界条件として固定している．このときの計算コ性が検証される［１６］．このとき，磁束関数 'でラベル付け

ストは，輸送コードの時間ステップごとにジャイロ運動論

される同心円磁気面上において，乱流揺動 !は，
$
$
!" &

線形解析を行う場合と比べ，非常に少ない［１５］．

'!
&)
&"#
%("&)
!'(
$! ! !$!'($$!($!&

２．
２．
９ フラックス・チューブ連結モデルの開発

（２）
$
!
')
%("&.
')(
')
&!&!(
"#
+
#$&%*
!#
-

磁場閉じ込めプラズマにおける乱流輸送過程を解析するとして表される．ここで，'はトロイダルモード数，)(')

有用な手法として，ジャイロ運動論的シミュレーションがは安全 係 数 で あ り，&!($$#'"&#)，#)
$
!
')
%"
'#)，$
$$
は動径方向のシミュレーション領域にわたる )の変化量として定義されている．この方法では，バルーニング構造をもつ乱流揺動を良く表現できるが，有限磁気シアのもとでは，,&,の増加とともに磁力線垂直方向の波数 %' が大きくなるという性質をもつ．このことは，時間ステップ幅へより厳しい制限を課すとともに，バルーニング表示がもつ（整数 (についての）対称性を破るという数値的な困難を引き起こす．この問題を解決するために，"& $"としたま
ま，複数（"*個）のフラックス・チューブを磁力線方向に連結するという新しい数値モデルを考案した．そこでは，
乱流揺動は次式のように表される（以下で &は整数）．
&

&)
!'(
'!
&)$ ! !$"&$$"" !'(
*

$$!&

')(
')
&!&!(
$
$
$
"
"*!
')
#$&%*
!#
"&
.
+
&

&

（３）

&!#
%&
)
$ ! !$!'(
$$!&

')(
')
&!&!(
$
!
')
%&
#$&%*
!#
!#
.
+
．

図１
０ LHD 実験結果とジャイロ運動論解析に基づく乱流輸送モデ
ルを用いた輸送シミュレーション結果との比較．

８０１

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016
$

%!"
$$
%&であり，またポロイダ&#!#"$##"" !%%ここで !#!%%の結果，多粒子種衝突モデルのゾーナルフローに対する衝ル角は "$$$%$"$%$"!&の範囲に制限される．こうして突減衰効果は，従来の簡約化モデルに比べて相対的に弱い構築されたモデルを用いることで，%の増加による時間ス

寄与を与えることが明らかとなった．また，電子−イオン

テップ幅の制限が緩和され，またバルーニング表示が持つあるいはイオン−電子間の衝突に代表されるような質量比

対称性も厳密に保たれる．

が著しく異なる粒子間衝突においては，速度空間上の衝突&

このフラックス・チューブ 連 結（flux tube train）
モデルを

周波数が急峻化し，数値微積分に関する誤差が増大する．

用いて，低磁気シアトカマク配位でのイオン温度勾配乱流

このような問題を解決するため，数値誤差を打ち消すよう

シミュレーションを行った．異なる"%および"&に対して，修正された衝突演算子を考案することにより GKV におけ

定常状態で得られた輸送係数を比較すると "%，"&$$にる多粒子種衝突モデルをさらに拡張した［２１］．修正された

おいて，従来の flux tube モデルの値と良く一致する結果がモデルは電子を含む多種イオンプラズマに対する熱緩和シ得られた（図１１）．一方，"% #!#の場合，従来のモデルでミュレーションに適用され，衝突計算に関する数値誤差がは解を得ることができなかったが，新しい手法では特段の丸め誤差程度の極めて小さいレベルに抑制され，粒子種間困難なく輸送係数の "&に対する収束性が確認できた．この温度緩和が正確に再現されていることが確認された．開の結果は，flux tube train モデルの優位性を示しており，他発された多粒子種衝突モデルは，今後の水素同位体やヘリの乱流輸送問題においてもその有効性の検証を進めていウム灰などの多種イオンで構成される燃焼プラズマ中の乱る．

流輸送解析に適用される．

２．
２．
１０ 高精度多粒子種衝突モデルの開発と乱流シミュレーションコード GKV への実装

２．
２．
１１ 国際ベンチマークの推進
ジャイロ運動論的乱流シミュレーション・コード GKV

水素同位体と電子が混合されたプラズマや不純物を含むの信頼性のさらなる向上をめざし，非軸対称系配位におけプラズマに対する乱流輸送シミュレーションへ拡張するたる乱流輸送シミュレーションのベンチマーク試験をプリン

めに，任意の異種粒子種間衝突の取り扱いを実現する多粒ストン・プラズマ物理研究所（PPPL）と共同で進めた．運子種衝突モデルを新たに開発し，ジャイロ運動論的シミュ

動論的手法に基づくプラズマ乱流シミュレーションは，

レーションコード GKV へ実装した．異なる平衡温度や平様々なコードによる研究が世界で行われているが，コード均流を持つ粒子種に対しても粒子数保存，運動量保存，エ間のベンチマーク試験が，特にステラレータでは遅れていネルギー保存を満足し，かつ，H 定理と呼ばれる衝突緩和

る．これまでに，PPPL やマックスプランク研究所（ドイ現象に関する物理的要請も満足する衝突モデル［１９］に基づツ）のグループと，NCSX のモデル配位におけるベンチ

き，GKVに従来実装されていた簡約化された同種粒子種衝

マークを行っていたが，今回は，LHD の実験配位を対象に突モデルを高精度の多粒子種衝突モ デ ル へ と 拡 張した

PPPL の GS2 コードと線形イオン温度勾配（ITG）不安定性［２０］．平衡分布への線形緩和に関するシミュレーションにの定量比較を進めた．ベンチマークでは，不安定性の線形より，粒子速度に依存した衝突周波数を含む多粒子種衝突

成長率や実周波数の波数スペクトルの比較に加え，固有関モデルの数値的健全性が確認された．また，線形ゾーナル

数の詳細比較も実施した［２２］．図１３は，GKV コードおよ

フロー減衰シミュレーションにより，残留ゾーナルフローび GS2 コードの双方から得られた ITG 不安定性の固有関数

レベルに対する衝突モデルの影響を検証した（図１２）．そをプロットしたものである．双方が，定性的な関数構造だ

図１
１ 従来のフラックス・チューブ・モデルと新しいフラック
ス・チューブ連結モデルを用いたイオン温度勾配乱流のシ
ミュレーションによって得られたイオン熱輸送係数の比較
［１
８］
．

図１
２ 従来の簡約化衝突モデル（LB）
および新たに開発した多粒子種衝突モデル（GK）を用いたゾーナルフロー減衰シミュ
レーション［２０］．

８０２Project Review

2.2 Transport Simulation of Core Plasmas

R. Kanno et al.

けでなく定量的にも非常に一致していることが分かる．ベ

が開始され，現在は九州大学・応用力学研究所で開発が進ンチマーク解析では，GS2 コードにおける時間積分法（陽められている乱流計測シミュレータ［２３］を適用し，実験で的解法または陰的解法）の違いにより，結果にわずかな差の乱流揺動計測を模擬することで，シミュレーションの異が生じ得ることも指摘され，今後のより詳細なベンチ

Validation 解析を試みた．これまでに，GKV コード独自のマーク研究への具体的な指針を得た．

座標系から実空間へのマッピング処理が完成し，乱流揺動分布の計測に威力を発揮する２次元位相コントラスト・イ

２．
２．
１２ 乱流計測シミュレータによる数値乱流計測

メージング（2D-PCI）による数値計測を開始した．本手法

磁場閉じ込めプラズマ研究において，乱流輸送機構の定により，乱流揺動のスペクトル構造や空間分布を実験と同量予測は重要な課題の一つである．ジャイロ運動論的シ様に調べることが可能となり，ジャイロ運動論的シミュ

ミュレーションによる数値解析は，物理の理解に向けた強レーションにおける定量的な乱流解析の手法開発が推進さ力な手法の一つである一方，計算結果に対する Validation

れた．本手法の基本部分は，現在，九州大学にて継続開発

解析も定量的な輸送評価には必須である．プラズマ乱流シされており，今後は，共同研究等を通じて，数値実験炉研

ミュレーション研究では，数例のみであるが，LHD プラズ

究プロジェクトにおいて，本手法の利用推進を図っていく．

マにおける熱輸送係数などの実験結果を定量的に再現できくことは，非局所的な効果等，更なるモデル改良の指針を

２．
２．
１３ ゾーナルフローに対する新しいクロー
ジャーモデルを用いたジャイロ流体 ITG
乱流シミュレーション

与えることが可能になる．ここでは，LHD 配位において５次元位相空間上の粒子分布関数を直接計算するジャイ

GKV コードを用いて実施した多くのジャイロ運動論的乱

ロ運動論的シミュレーションに比べ，３次元実空間上の流流シミュレーションの計算データに対して，NIFS で開発

体変数を求めるジャイロ流体シミュレーションは，計算容ている中，このような数値解析結果について，実際の実験で行われるような手法による乱流模擬計測法を確立してお量・計算時間が遙かに少ない点で魅力的であるが，その妥当性は，ランダウ減衰・有限ジャイロ半径・ジャイロ中心
軌道等の運動論的効果を記述するために用いられる流体変数のクロージャーモデルに大きく依存する．Dimits等［２４］によって示されているように，イオン温度勾配（ITG）乱流のジャイロ運動論的シミュレーションの結果に比べ，従来のジャイロ流体シミュレーションは，ITG 乱流により生成されるゾーナルフローの大きさを過小評価し，ITG 乱流によるイオンの乱流熱輸送を過大評価する傾向にある．
ゾーナルフローの生成過程を正確に記述するため，洲鎌等により新しいゾーナルフロークロージャーモデルが導出され［２５，２６］，それを用いたジャイロ流体 ITG 乱流シミュ
レーションが行われた［２７］．
図１４は，Cyclone DIII-D base case parameter set［２２］に対応するトカマクプラズマにおけるジャイロ流体ITG乱流シ
ミュレーションから得られた静電ポテンシャルのポロイダ
ル断面上の分布構造を表す．図１４（左）および（右）は，それぞれ，ゾーナルフローに対するクロージャーモデルを用いていない場合と用いた場合を示している．ゾーナルフ
ロークロージャーモデルを用いた場合の方が，大きなゾー
ナルフローが生成されることがわかる．
図１５（左）は，図１４の場合に対応するジャイロ流体 ITG
乱流シミュレーションから得られたイオン熱拡散係数の時間発展を表す．また，図１５（右）は，ジャイロ流体 ITG 乱流
シミュレーションから得られたイオン熱拡散係数のイオン
温度勾配依存性を表す．ゾーナルフロークロージャーモデ
ル［２５］を用いていない場合に比べ，用いた場合の方が，イオン熱拡散係数が少なくなり，ジャイロ運動論的シミュ
レーションの結果とよく一致している．また，図１
５（右）より，ITG 乱流輸送が現れてくるイオン温度勾配の値が ITG

図１
３ LHD 高イオン温度放電配位における線形イオン温度勾配不安定性における固有関数の比較（a）とその拡大図（b）．曲線が GS2 の結果，点（＋）
が GKV の結果を示す．

モードが線形不安定となるイオン温度勾配の値 "!
!"#!!
よりも上方にずれるという Dimits シフトの様子が，ゾーナ
８０３

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

く新古典輸送コードの比較により，従来から知られていたグローバル計算と局所近似の新古典フラックスの径電場依存性の差が，局所近似モデルで無視している
磁場ドリフト項の磁気面接線成分によって生じることを示した．この磁場ドリフト項の接線成分を残した新
しい局所近似モデルを構築し，グローバル計算で見られる傾向を少ない計算量で再現することが可能となった．また，新古典輸送シミュレーションをいくつかのヘリカル装置における実験解析に応用し，熱輸送や径
図１４ Cyclone DIII-D base case parameter set を用いたジャイロ
流体 ITG 乱流シミュレーション［２７］から得られた静電ポテ
ンシャルのポロイダル断面における分布構造．ゾーナルフ
ローに対するクロージャーモデル［２５］を用いていない場合（左）と用いた場合（右）［参考文献［２７］，Fig. 4（著作権：AIP Publishing）より許可を得て転載］．

電場の観測結果との比較を通じた Verification & Validation 研究が進められた．これらの研究を通じて，コ
アプラズマにおける輸送現象における新古典輸送部分の定量的予測性の向上に進展がみられた．
３）共鳴摂動磁場（RMP）を静的に加えた領域におけるプ

ルフロークロージャーモデルを用いることによりよく再現

ラズマの衝突輸送に対する理解のため，Rechester-

されていることがわかる．

Rosenbluth（R-R）理論と同様に電場を無視する仮定の下で実行した !
ドリフト運動論的シミュレーション!

２．
２．
１４ まとめにより，（１）径方向のイオン熱拡散係数は，R-R 理論値

数値実験炉研究プロジェクト新古典・乱流輸送シミュと比べて極めて小さくなることを確認し，（２）その物

レーショングループにおけるコアプラズマの輸送研究で理パラメータ依存性については，R-R 理論モデルとほは，コアプラズマの新古典および乱流輸送現象に対して，

ぼ同じであることが分かった．また，（３）RMP の効果

ドリフト運動論およびジャイロ運動論に基づく理論モデル

は，衝突頻度が大きくなるに従い弱められることが明の適用や大規模シミュレーション解析を進め，平成２５年度

らかになった．

から２７年度の期間においては，主に以下の成果を得た．

４）実際の LHD 実験や JT-60U 実験における有限ベータプ

１）LHD のような３次元磁場配位用に開発された新古典

ラズマに対する ITG および TEM 乱流輸送シミュレー

輸送コードを，トカマクの誤差磁場による軸対称性のション解析がジャイロ運動論的乱流シミュレーション

破れによって生じる新古典トロイダル粘性の評価に応

コード GKV によって実行され，コア領域におけるイ

用し，JT-60UのNBI実験において観測されたトロイダ

オンおよび電子熱輸送フラックスと温度勾配分布に関

ル回転分布を数値シミュレーションで再現性良く計算する実験結果を，計測誤差程度の範囲で再現する結果することが可能になった．

が確認された．さらに，GKV による ITG 乱流シミュ２）グローバル及び局所近似ドリフト運動論モデルに基づ

レーション結果から構築された簡約化輸送モデルを統

図１
５ Cyclone DIII-D base case parameter set を用いたジャイロ流体 ITG 乱流シミュレーション［２７］から得られたイオン熱拡散係数の時間
発展（左）とイオン温度勾配依存性（右）［参考文献［２７］，Fig. 3 および Fig. 5（著作権：AIP Publishing）より許可を得て転載］．左図の（w ZF）の曲線と右図の丸印で表したデータ点は，ゾーナルフローに対するクロージャーモデル［２５］を用いた場合を表し，左図の（w/o ZF）の曲線と右図の菱形で表したデータ点は，そのモデルを用いていない場合を表す［２
７］．白い丸印・菱形と黒い丸印・菱形では，有限ジャイロ半径効果を表すために少し異なるモデル係数を用いているが，その結果得られたイオン拡散係数にほとんど差はな
いことがわかる．

８０４

Project Review

2.2 Transport Simulation of Core Plasmas

R. Kanno et al.

［５］T.E. Evans et al., Nature Phys. 2, 419 (2006).
［６］R. Kanno et al., Plasma Phys. Control. Fusion 55, 065005(2013).
［７］A. Ishizawa et al., Nucl. Fusion 53, 053007 (2013).
［８］A. Ishizawa et al., Phys. Plasmas 21, 055905 (2014). アメリ
カ物理学会の招待講演論文，Physics of Plasmas 誌のEditor's Choice に選ばれた．
［９］H. Sugama et al., Phys. Plasmas 16, 112503 (2009).
［１０］A. Ishizawa et al., J. Plasma Phys. 81, 435810203 (2015)
レビュー論文，おすすめ論文に選ばれ Free download
可能．
［１１］A. Ishizawa et al., Nucl. Fusion 55, 043024 (2015)［１２］M. Nakata et al., Plasma Fusion Res. 9, 1403029 (2014) .
［１
３］M. Nakata et al., Proc. IAEA-FEC2014, TH/P7-38 (2014).
［１４］M. Nakata et al., Nucl. Fusion 56, 086010 (2016).
［１５］S. Toda et al., J. Phys.: Conf. Ser. 561, 012020 (2014).
［１６］M. A. Beer et al., Phys. Plasmas 2, 2687 (1995).
［１７］T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
［１８］T.-H. Watanabe et al., Phys. Plasmas 22, 022507 (2015).
［１９］H. Sugama et al., Phys. Plasmas 16, 112503 (2009)［２０］M. Nunami et al., Plasma Fusion Res. 10, 1403058 (2015) .
［２
１］M. Nakata et al., Phys. Commun. 197, 61 (2015) .
［２
２］D.R. Mikkelsen et al., Phys. Plasmas 21, 112305 (2014).
［２３］N. Kasuya et al., IAEA-FEC2010, THC/P4-15 (2010).
［２４］A.M. Dimits et al., Phys. Plasmas 7, 969 (2000).
［２５］H. Sugama et al ., Phys. Plasmas 14, 022502 (2007).
［２６］O. Yamagishi and H. Sugama, Phys. Plasmas 19, 092504(2012).
［２７］O. Yamagishi and H. Sugama, Phys. Plasmas 23, 032305(2016).

合コード TASK3D へ適用し，LHD 実験におけるイオ
ン温度分布が高い精度で再現された．
５）不純物を印加するプラズマ実験や ITER/DEMO における燃焼プラズマのように複数イオン種で構成されるプラズマに対する乱流輸送シミュレーションを構築するため，保存特性を満たす多粒子種衝突モデルの開発
や，乱流渦の沿磁力線構造をより正確に扱うフラック
ス・チューブ連結モデルの開発を行い，GKVに適用することでそれぞれの有効性の検証が成された．
６）国内および国際共同研究により，乱流計測シミュレー
タによる数値乱流計測の開発や，複数の非軸対称系配位を対象にしたクロスコードベンチマークが実施され，開発した輸送シミュレーションコードの verification や揺動計測などを活用したより詳細な validationの基盤構築が進展した．
７）新しいゾーナルフロークロージャーモデルを用いた
ジャイロ流体 ITG 乱流シミュレーションを行い，そのクロージャーモデルを用いない場合と比べ，ジャイロ
運動論的シミュレーションの結果とよく一致することを確認し，また，Dimits シフトの様子が，よく再現されることを明らかにした．

参考文献［１］S. Matsuoka et al., Phys. Plasmas 22, 072511 (2015).
［２］M. Honda et al., Nucl. Fusion 55, 073033 (2015).
［３］S. Kitajima et al., Nucl. Fusion 53, 073014 (2013).
［４］S. Satake et al., 20th International Stellarator- Heliotron
Workshop, Oct. 7 2015, Greifswald, Germany, S1-I5.

８０５

J. Plasma Fusion Res. Vol.92, No.11 (2016)８０６‐８０９

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

２．核融合プラズマシミュレーション研究
2. Simulation Researches in Fusion Plasmas

２．
３ 高エネルギー粒子・波動・加熱の物理
2.3 Physics of Energetic Particles, Waves and Heating藤堂泰１），WANG Hao１），西村征也２），横 山 雅 之１）

TODO Yasushi１）, WANG Hao１）, NISHIMURA Seiya２) and YOKOYAMA Masayuki１）
１）

自然科学研究機構核融合科学研究所，２）神戸市立工業高等専門学校（原稿受付：２
０
１
６年７月１
３日）

高エネルギー粒子・MHD 連結シミュレーションコード MEGA を用いて，環状プラズマにおける高エネル
ギー粒子駆動不安定性のシミュレーション研究を推進している．大型ヘリカル装置（LHD）
における高エネルギー
粒子駆動不安定性については，アルフベン固有モードによる高速イオン輸送のポアンカレ断面図解析と実形状３次元平衡データを用いた高エネルギー粒子駆動測地的音響モードのシミュレーションが進展した．実験で観測される高速イオン分布とアルフベン固有モードを再現するためには高速イオン減速時間スケールの長時間計算が必要である．そのような計算を可能にする高エネルギー粒子・MHD 連結マルチフェーズシミュレーションを開発
し，DIII-D トカマクでの実験に適用してその妥当性を実証した．ITER の２つの運転シナリオについて高エネル
ギー粒子駆動不安定性のシミュレーションを実行し，安定性と高エネルギー粒子輸送を明らかにした．以上の成果に加えて TASK3D コードを中心とした LHD における加熱解析の進展について報告する．
Keywords:
energetic particle, hybrid simulation, Alfvén eigenmode, geodesic acoustic mode, heating analysis

２．
３．
１ 高エネルギー粒子の物理（EGAM）も観測さ れ て い る．MEGA コ ー ド を 用 い て，

（１）LHD におけるアルフベン固有モードによる高速イオ

LHD と同等の物理条件をもつ軸対称トカマク平衡でのン輸送と高エネルギー粒子駆動測地的音響モード

EGAM の線形特性や非線形時間発展における周波数変調

LHD においてバースト的に繰り返し発生するアルフベを研究してきた［４，５］．最近では中性粒子との荷電交換を

ン固有モードによって高速イオンが規格化小半径の１０％程考慮した高速イオン分布関数を用いてシミュレーションを度輸送されたことが中性粒子解析器（NPA）の計測により実行し，LHD で観測された周波数がバルク温度に依存しな示唆されている［１］．この実験における２つのアルフベン

い EGAM を再現することができた［６］．さまざまな初期条

固有モードを AE3D コードを用いて解析し，それらのアル件を調べた結果，バルクプラズマの密度が低く，高速イオ

フベン固有モードを一定振幅で振動する擾乱とした際の高ン減速時間に対して荷電交換時間が十分に短い場合に，

速イオン軌道のポアンカレ断面図解析を行った［２］．ここEGAM の周波数がバルク温度に依存しないことがわかっでは HINT2 コードで解析した３次元平衡データを使用した．また，HINT2 コードによる LHD の３次元平衡データ

た．図１は LHD の横長断面における高速イオン軌道ポアを使用したEGAMのシミュレーションを開始した．このシンカレ断面図を示しており，アルフベン固有モードの振幅

ミュレーションで得られた EGAM の周波数は実験とほぼが増大すると軌道の乱れた領域が拡大することがわかる．

等しい 61 kHz であった．

図２は，この乱れた領域の半径方向の幅をアルフベン固有（２）高エネルギー粒子・MHD 連結マルチフェーズシミュ
レーションの開発と検証

モードの振幅に対して示している．実験結果から示唆された規格化小半径の１０％程度の輸送が発生するのは，アルフ

ITER などの核燃焼プラズマにおいて，アルファ粒子を

ベン固有モードの振幅が平衡磁場に対して!""!!##の場合

はじめとする高速イオンがアルフベン固有モード（AE）をであることがわかった．この結論は，LHD と同様の回転変

共鳴相互作用によって不安定化する可能性がある．核燃焼

換分布を持つ軸対称トカマク型平衡を用いた解析［３］と整プラズマにおいて高速イオンは粒子衝突を通して燃料プラ合している．

ズマを加熱する役割があり，高速イオンの分布変化と損失

LHD で は 高 エ ネ ル ギ ー 粒 子 駆 動 測 地 的 音 響 モ ー ドをもたらす AE は核燃焼プラズマに関する重要な研究課題
corresponding author’s e-mail: todo@nifs.ac.jp

８０６!2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchProject Review

2.3 Physics of Energetic Particles, Waves and Heating

Y. Todo et al.

の一つである．現在のトカマク，ヘリカルおよび球状トカ

測されている．DIII-D においては電子サイクロトロン放射

マク装置においても，高速アルファ粒子の代わりに中性粒（ECE）計測によりAEによる電子温度揺動が詳細に測定さ子ビーム入射（NBI）加熱や波動（ICRF）加熱による高速れるとともに，高速イオンの空間分布も測定されている．

イオンを用いた AE に関する実験が活発に行われている．

このような実験においては，高速イオン分布が AE によっ

米国の DIII-D トカマクでは，定常的に小振幅を保つ多数て強く影響されていることに注意しなければならない．実の AE とそれらによる高速イオン分布の顕著な平坦化が観験で計測される高速イオン分布は，AE によって速度空間と実空間の双方において緩和されたものである．したがって AE と高速イオンの両者の時間発展を高速イオンの衝突
減速時間程度の長時間追跡するシミュレーションが必要である．
このために，高速イオンと MHD を連結して両者の時間
発展を自己整合的に追跡するハイブリッドシミュレーショ
ンコード MEGA に現実的なビームイオン入射分布および粒子衝突と損失を導入した．さらに，MHD 揺動を含む計算と含まない計算を交互に行うマルチフェーズ計算により高速イオン減速時間程度の長時間の計算を可能にし［７］，
実験との詳細な比較のため熱イオン反磁性ドリフトとプラ
ズマのトロイダル回転を取り入れて MHD モデルを拡張した［８］．図３にシミュレーションにおいて多数のアルフベ
ン固有モードと相互作用して形成された高速イオン圧力分布を実験と比較する．両者が実験の誤差程度で一致していることと，MHD 揺動を含まずビームイオン入射と粒子衝
突だけの古典計算結果と比較すると両者が強く平坦化していることがわかる．図４はシミュレーションにおいてトロ!"のアルフベン固有モードがもたら
イダルモード数 !"!
す電子温度揺動のポロイダル面分布を示している．図４の電子温度揺動は捻れた構造を示している．この揺動は時計回りに回転しているので，プラズマ外縁部は振動を先導していることになる．図５はこの電子温度揺動の振幅と位相

図１

図２

アルフベン固有モードが存在する条件下での高速イオン軌
道ポアンカレ断面図；（a）
トロイダル磁場で規格化した磁場揺動振幅 8.9 × 10−4，（b）
7.1 × 10−3［２］．

図３ （a）
高速イオン蓄積エネルギー時間発展と（b）
高速イオン圧力 分 布 に 関 す る 高 エ ネ ル ギ ー 粒 子・MHD 連 結 マ ル チ
フェーズシミュレーションと MHD 揺動を含まない古典計
算および実験計測の比較．（b）
ではマルチフェーズシミュ
レーションが実験の高速イオン圧力分布をよく再現している［８］．

アルフベン固有モード振幅に対する乱れた軌道領域幅［２］．

８０７Journal of Plasma and Fusion Research Vol.92, No.11 November 2016の分布を実験と比較したものである．図５では，振幅分布用した．有限ラーモア半径効果がアルファ粒子とビームイが絶対値を含めてよく一致していることに加えて，位相分オンの双方について考慮されている．定常運転シナリオの布についてもよい一致が得られた．このように高速イオン

シ ミ ュ レ ー シ ョ ン の 結 果，ト ロ イ ダ ル モ ー ド 数

分布とアルヴェン固有モード分布について，実験結果の再"!!""!#のトロイダルアルフベン固有モードの線形成

現に成功したのは世界初の成果である．

長率が最も高いが，非線形局面では "!"のベータアルフ（３）ITER における高エネルギー粒子駆動不安定性と高エベン固有モードが支配的となることがわかった．図６に線

ネルギー粒子輸送

形局面と非線形局面でのプラズマ半径方向の MHD 速度揺ASTRA コードと EFIT コードによって計算されたプラ動分布を示す．これらのアルフベン固有モードの磁場揺動

ズマ電流９MAの定常運転シナリオと15 MAの標準運転シ

飽和振幅はトロイダル磁場で規格化して 0.4％ であった．

ナリオについて，高エネルギー粒子駆動不安定性と高エネ

図７に示すように，アルフベン固有モードがもたらすアル

ルギーアルファ粒子輸送を調べた［９］．高エネルギー粒子

ファ粒子と重水素ビームイオンの圧力変化は，それぞれのとしてはアルファ粒子だけでなく，入射エネルギー 1 MeV

プラズマ中心値と比較して６％および８％にとどまった．

の重水素ビームイオンを考慮した．高エネルギーアルファ

プラズマ電流 15 MA の標準運転シナリオについても粒子の速度空間分布は等方的であるが，NBI による高エネ

MEGA コードを用いたシミュレーションを実行した．ここルギーイオンは非等方的である．円柱座標系 #
!!
"!
$
$に対

で 用 い た ASTRA コ ー ド と EFIT コ ー ド に よ る 平 衡 でして そ れ ぞ れ
（２５６，
２５６，
５１２）個 の 格 子点 を 使 用 し，ア ル

は，安全係数 #!!の磁気面が３箇所に存在している．シ6

ファ粒子とビームイオンにそれぞれ 8.4×10 個の粒子を適

ミュレーションの結果，プラズマ電流とバルクプラズマ圧

図４

マルチフェーズシミュレーションにおける電子温度揺動のポロイダル面分布．各パネルはそれぞれ（a）
n = 3，（b）
n = 4，（c）
n =5のアルフベン固有モードによる揺動を示している．カラーバーの単位は電子ボルトである［８］．

図５電子温度揺動分布に関するマルチフェーズシミュレーションと実験データの比較．上段（a）‐（c）
は振幅分布，下段（d）
‐
（f）
は振動の位
相分布を示している．最も振幅の大きい n = 3 モードは振幅（a）と位相（d）
がともに実験とよく一致している．その他の n = 4（b および e），５（c および f）モードについてもよい一致が見られる［８］．

８０８Project Review

2.3 Physics of Energetic Particles, Waves and Heating

図７

図６

ITER 定常運転シナリオにおける（a）
線形局面および（b）
非線形局面でのアルフベン固有モードの半径方向速度揺動分
布［９］．

Y. Todo et al.

ITER 定常運転シナリオにおける（a）
高エネルギーアルファ
粒子とビーム重水素イオンの初期圧力分布（破線）と最終圧力分布（実線）の比較，および（b）
圧力変化の分布［９］．

01 では，NBI 加熱プラズマの解析に限定されていたが，TASK3D-a02 によって，ECH 単独，あるいは，重畳プラズ

力が駆動するトロイダルモード数!!"のMHD不安定性がマにも適用範囲が格段に拡がった．

発生し，アルファ粒子の強い分布変化が発生した．アル

ICRF加熱に関しては，TASK3D-a開発の一環として，京ファ粒子圧力の減少は最大でプラズマ中心値の２５％程度に

都大学：福山淳教授との共同研究によって，波動伝搬解析達した．安全係数を一様に 0.1 引き上げて "!!の磁気面を

コード TASK/WM の誘電テンソルの改良，吸収計算の高除去した平衡を初期条件としてシミュレーションを実行し精度化のための電場成分記述法の変更，径方向の有限差分たところ，発生したのは弱い MHD 不安性のみであり，高解法への変更などの改良を行っている．核融合研：關良輔

エネルギー粒子の分布変化は無視できるほど小さいことが助教によって，２０１５年の土岐コンファレンスにおいてそのわかった．従ってこれらのシミュレーションからは，"!!

初期成果について報告が行われた［１２］．LHD 実験への適の磁気面が存在しなければ，標準運転シナリオでは高エネ用をめざして，コード改良を加速する．

ルギー粒子駆動不安定性および高エネルギー粒子輸送の問

参考文献

題はないと判断できる．

［１］M. Osakabe et al., Nucl. Fusion 46, S911 (2006).
［２］S. Nishimura et al., Plasma Fusion Res. 8, 2403090 (2013).
［３］Y. Todo et al., Plasma Fusion Res. 3, S1074 (2008).
［４］H. Wang and Y. Todo, Phys. Plasmas 20, 012506 (2013).
［５］H. Wang et al., Phys. Rev. Lett. 110, 155006 (2013).
［６］H. Wang et al., Phys. Plasmas 22, 092507 (2015).
［７］Y. Todo et al., Nucl. Fusion 54, 104012 (2014).
［８］Y. Todo et al., Nucl. Fusion 55, 073020 (2015).
［９］Y. Todo and A. Bierwage, Plasma Fusion Res. 9, 3403068(2014).
［１０］T. Ii Tsujimura et al., Nucl. Fusion 55, 123019 (2015).
［１１］N. Marushchenko et al., Plasma Fusion Res. 2, S1129(2007).
［１２］R. Seki et al., "Development of TASK3D/WM module for
evaluation the ICRF wave propagation and absorption in
LHD", 25th International TOKI Conference (Nov. 3-6,
2015, Toki, Japan) P2-45.

２．
３．
２ LHD における加熱解析の進展
ECH 加 熱 分 布 を 評 価 す る LHDGauss コ ー ド お よ び
TRAVIS コードの統合輸送解析スイート TASK3D-a02 への導入作業が完了した（２．
５節参照）．LHDGauss コードは
核融合研オリジナルのコードであり，最近，ほぼショット
間での ECH 加熱解析が可能となり，ECH 加熱の最適化に活用されている［１０］．一方，TRAVIS コード［１１］は，マックスプランクプラズマ物理研究所：N. Marushchenko 博士
によって開発され，核融合研：吉村泰夫准教授によって，
単 体 と し て LHD 実 験 解 析 に 用 い ら れ て き て い た．
TASK3D-a02 は，LHD解析データサーバに登録されている
LHDGauss の解析結果の参照，および，TRAVIS のその場実行を行い，電子温度分布計測がある全時間タイミングでの ECH 吸収パワー計算結果を取り込んだエネルギーバラ
ンス解析が行えるように組み上げられている．TASK3D-a

８０９

J. Plasma Fusion Res. Vol.92, No.11 (2016)８
１０‐８１３

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

２．核融合プラズマシミュレーション研究
2. Simulation Researches in Fusion Plasmas

２．
４ 周辺プラズマの輸送とプラズマ・壁相互作用
2.4 Peripheral Plasma Transport and Plasma-Wall Interaction
鈴 木 康 浩１，２），河 村 学 思１，２），中 村 浩 章１，３），斎 藤 誠 紀４），伊 藤 篤 史１，２），高 山 有 道１）
SUZUKI Yasuhiro１，２）, KAWAMURA Gakushi１，２）, NAKAMURA Hiroaki１，３）, SAITO Seiki４）,
ITO M. Atsuhi１，２）and TAKAYAMA Arimichi１）
１）

自然科学研究機構核融合科学研究所，２）総合研究大学院大学，３）名古屋大学，４）釧路工業高等専門学校（原稿受付：２
０
１
６年７月１
３日）

本章では，粒子や熱が周辺プラズマ中でどのように輸送され，第一壁やダイバータ材料をどのように損傷させるかを議論する．周辺プラズマ部での輸送に関しては，最近の EMC3-EIRENE コードを用いた解析結果を中心に紹介する．プラズマ・壁相互作用については分子動力学を用いたシミュレーションの最新の成果について紹介する．
Keywords:
peripheral plasma transport, EMC3-EIRENE code, closed divertor, impurity, stochasticity,chaos coordinates, anisotropic heat transport, plasma wall interaction, tungsten, fuzz structure,
molecular dynamics, binary collision approximation, density functional theory, kinetic Monte Carlo method,
divertor tile, plasma facing material

２．
４．
１ 周辺プラズマの輸送シミュレーション

ものであり，粒子のリサイクリングを含めた解析へのバリ

周辺プラズマは，大型装置における熱負荷予測やその低デーションの一つと位置づけられる結果である．この成果

減手法，炉材料の損耗や再堆積など，工学要素と関わりがは第１４回 Plasma Edge Theory 国際ワークショップで報告

深く，重要度が増している．不純物を含めたプラズマおよを行った．

び中性粒子の輸送モデリングが必要とされているが，十分

２．
４．
１．
２ ガス排気と燃料供給を含めた粒子輸送解析に解析手法や計算モデルが確立されていないため，そのモ

粒子輸送は，能動的な排気と供給以外にダイバータ板でデリングは喫緊の課題である．このような状況から，周辺の表面再結合過程があるため，自己無撞着な取り扱いが必

輸送モデリングの基盤として，EMC３‐EIRENEコードを中

要である．解析のため，EMC３‐EIRENE に排気ポンプと２心にモデル開発とそれを用いた解析を進めてきた．進展の種類の供給源を組み込んだ．ポンプはドーム構造裏側に粒あった代表的な成果について以下に述べる．

子吸収条件を設定し，その排気と釣り合う供給として，ガ

２．
４．
１．
１ ダイバータ配位の中性粒子輸送への影響解析

スパフを模擬してダイバータ板でのリサイクリングを増加

LHD の開ダイバータから閉ダイバータへの変更の影響させる方法と，炉心へのペレット供給や NBI の粒子デポジを予測・説明する輸送計算を行った．この目的のため，エルゴディック領域とダイバータレグ領域を覆う計算メッ
シュを EMC３‐EIRENE コード専用に開発した［１］．図１に中性ガス圧の比較を示す．閉ダイバータのドーム構造下部で１０倍程度のガス圧上昇が見られ，実験で測定されたガス圧上昇と２倍以内の範囲で一致した．さらに，加熱入力を一定に保って電子密度スキャンを行ったところ，図２に示すように，計測とほぼ同じ依存性が得られた．圧力の絶対値については複数の要因で変化しうるため，物理的に妥
図１

当な範囲内で合わせているのが現状だが，依存性（スケーリング）の一致は物理機構を含めたモデルの妥当性を示す

水素分子の圧力比較．（a）
開ダイバータ，（b）
閉ダイバー
タ．気体温度は室温を仮定した．
corresponding author’s e-mail: suzuki.yasuhiro@LHD.nifs.ac.jp

８１０!2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchProject Review

2.4 Peripheral Plasma Transport and Plasma-Wall Interaction

図３
図２

水素分子圧力の計測との比較．赤＋：開ダイバータ計測
値，青×：閉ダイバータ計測値．実線はそれぞれの計算値．

Y. Suzuki et al.

水素分子の圧力比較．（a）
ガス供給，（b）
炉心供給．配位はともに閉ダイバータ．

あることがわかった．また，他の線状のピークについては，カメラから遠い位置にある同様のネオン密度ピークの寄ションを模擬して炉心側境界面からイオンとして供給する

与，およびプラズマの縁に沿った視線積分の寄与によって方法を用いた．ガスの排気も供給もない場合は，前項の生じていることがわかった．ネオンは炭素に比べてイオン

図１（b）が相当する．一方，排気と２種類の供給を行った

化エネルギーが高いため，より深くプラズマに進入するこ場合は図３の結果となった．ガス供給の場合，ポンプのあとが確認された．その結果として，プラズマの流れによるるドーム下部を除き，排気なしとほぼ同じ結果となる一摩擦力よりも温度勾配による熱力の影響を強く受け，ネオ方，炉心供給を行った場合，全体的にガス圧が減少する結ンは炭素に比べてエルゴディック領域に蓄積しやすく，よ果となった．これは，炉心密度を一定に保つ条件のためにり広い発光分布を持つことがわかった．

リサイクリング粒子束が減少し，ダイバータでの電子密度

２．
４．
１．
４ 乱れた磁力線領域での熱輸送に関する理論研究の減少とガス圧の減少が起きたためである．周辺の粒子輸

大型ヘリカル装置（LHD）をはじめとするヘリオトロ

送には粒子供給方法を含めたモデリングが重要であること

ン・ステラレータ磁場配位では周辺部にストカスティックが示された．この成果は第２４回 Fusion Energy Conference

な磁力線構造が存在する．このストカスティックな磁力線で報告を行った［２］．

構造はプラズマ応答により増幅するものと考えられてお

２．
４．
１．
３ 不純物ガスパフ時の不純物輸送解析

り，プラズマ周辺部の輸送解析を行うためにはストカス

第１７サ イ ク ル で 行 わ れ た ネ オ ン ガ ス パ フ 実 験

ティックな磁力線構造を含む解析を行う必要がある．一（＃１
１７４７８）を対象に，トムソン散乱計測で得られたプラズ

方，トカマクの場合は軸対称性を仮定することにより，通マ分布を再現するように輸送係数等の入力パラメータの決

常，２次元で周辺輸送の解析が行われる．しかし，閉じ込定を行い，スパッタリングで生じる炭素と，ガスパフで導め改善モード（H モード）に発生する周辺局在化モード

入するネオンの２成分不純物の輸送を水素プラズマと同時（ELM）を低減・抑制するために共鳴磁場摂動（RMP）をに解いた．図４に計測との比較を示す．縦長断面を上部のプラズマに重畳する実験が盛んに行われている．このこと

6.5 U ポートから少しトロイダル方向に傾く視線で観測しは，これまで２次元で解析が行われてきたトカマクにおいたものである．加熱入力 11.6 MW に対して，水素・炭素・ても３次元の輸送解析が必要であることを意味しており，ネオンの合計放射パワーが計測値の 2.5 MW と一致するよ

研究の新しい展開として注目されている．

うにネオンの導入量を調整した．発光分布の広がりは異な

これまで，無衝突の極限ではストカスティックな磁力線

るが，ひときわ強いピーク位置および線状の構造は非常に領域ではプラズマの閉じ込め性能はゼロとされてきた．こ

よく一致している．解析から，強いピークは X 点近くにあれは，ストカスティックな磁力線領域では磁力線に沿った

る比較的低価数のネオンイオン密度のピークの寄与である

熱輸送が，温度を直ちに拡散させるためである．しかし，ことがわかった．さらに，このピークを持つ領域は最も粒

衝突プラズマでは有限な磁力線平行方向の熱拡散係数によ

子束の大きいトーラス内側のダイバータ板と磁力線でつなりある程度の閉じ込め性能を有するが，それでも十分低いがっており，磁場に平行な輸送で直接もたらされたものでと考えられてきた．一方，LHD 実験やトカマクの RMP

図４

発光強度の比較．（a）
ボロメータ計測，（b）
EMC3-EIRENE 結果．

８１１

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

実験では，磁力線構造がストカスティックと考えられるプ
ラズマ周辺領域で有限な温度勾配が観測されている．このことは，ストカスティックな磁力線構造が有限な温度勾配を維持できる可能性を示している．このことを理論的に考察するためにカオス座標と呼ばれる，ストカスティックな磁力線構造領域で構築可能な新しい座標系の研究を，プリンストンプラズマ物理研究所と共同で行っている．この座標系を用いて，磁力線の統計的性質を考察すると共に，ス
トカスティックな磁力線構造中の熱輸送の解析を行った．
カオス座標系とは，これまで一般に用いられてきた磁気座標系と異なり，磁気島やストカスティックな磁力線構造
中から周期的な磁力線軌道を探し出して，座標系として利用するものである．このことにより，磁気島のセパラト
リックス中やストカスティックな磁力線構造中から KAM
面やカントール面をロバストに探索することが可能になった．このカオス座標系を用いて磁力線構造の詳細な解析を行うと，ストカスティックな磁力線構造中に KAM 面やカントール面が存在する場合にはプラズマ小半径方向に変位するストカスティックな磁力線に対してバリヤーとして作用することが明らかになった．このことは，たとえ磁力線
図５構造がストカスティックになったとしても，KAM 面のバリヤーが存在すれば，そこで有限な温度勾配が維持できる

乱れた磁力線領域をもつトカマク配位での熱輸送解析結
果．実線は温度の等高線を表す．

可能性を示すものである．このことを確認するために現実的な値を持つ磁力線に平行方向と垂直方向の熱拡散係数をダイバータ材料であるタングステンにヘリウムの照射を仮定した熱輸送解析を，RMP を重畳したトカマク配位に続けると数十ナノメートルスケールのヘリウムバブルがタ

対 し て 行 っ た．磁 力 線 に 平 行 方 向 と 垂 直 方 向 の 比 をングステン内部に発生する．照射量がさらに上がりフルーκ$〜1010 に設定し，以下の非等方熱輸送方程式の定常
κ(!

エンスが 1025 m−2 以上に達すると繊維状ナノ構造が形成さ解を数値シミュレーションにより計算した．

れる［３］．これらのナノ構造が形成されたタングステン表面は熱耐久性等の性能劣化や，表面でのアークの発生など""
κ(%("!κ$%$"'
!!
#%"
&
"
#

が危惧される．また別の視点から，ヘリウムバブルは固体金属中での気泡の形成現象という点で，繊維状ナノ構造は

図５に結果を示す．図中には磁力線のポワンカレ図，セパ

スパッタリングでも堆積でもない新しいナノ物質形成現象

ラトリックス上の安定・不安定多様体（stable/unstableという点で科学的にも未知の現象であった．このことか

manifolds）を赤と青の太線，実線で温度の等高線を示し

ら，ヘリウムとタングステンの相互作用問題は，プラズマ

た．ポワンカレ図から見ると磁力線構造がストカスティッ−壁相互作用研究においても重要課題と認識されている．

クになっているにもかかわらず，有限な温度勾配が存在す

シミュレーションでヘリウムバブルと繊維状ナノ構造のることがわかる．このことは，ストカスティックな磁力線問題を取り扱うにあたり我々が直面した問題は，これらの構造中に KAM 面から構成されるバリヤーが存在し，有限構造の発現は異なる物理階層の現象が複合的に組み合わな温度勾配を維持していると推測される．このことは，プ

さって起こるものであり，何か単一のシミュレーション手

ラズマ周辺部の輸送を理解する上で，重要な知見であると法で計算できるものではないということであった．そこ考えられる．この結果は第２０回国際ステラレータ・ヘリオ

で，まずは現象を分類することから検討を始めた．幸い国

トロンワークショップで報告された．

内には繊維状ナノ構造を発見した愛知工業大学：高村秀一
教授［３］を始めとして，実験室プラズマ実験から顕微鏡に

２．
４．
２ 数値シミュレーションを用いたタングス
テン挙動の解析よる観察まで，繊維状ナノ構造に纏わる実験研究者が揃っている．これらの実験研究者との議論を経て，ヘリウムバ

数値実験炉プロジェクトのプラズマ・壁相互作用研究グ

ブルと繊維状ナノ構造の形成現象は，（１）ヘリウムプラズ

ループでは，磁場閉じ込め型核融合装置のダイバータならマの入射過程，（２）タングステン内部のヘリウムの拡散・びに第一壁におけるプラズマ−壁相互作用（PWI）現象を

凝集過程，（３）ナノバブルの成長過程，（４）繊維状構造の成理解する為，分子動力学（MD）を中心としたシミュレー

長過程の四段階の過程に分類した［４］．

ション研究を行ってきた．特に，タングステンがヘリウム第一段階となるヘリウムプラズマの入射過程では，ヘリプラズマ照射を受けた場合に発現するヘリウムバブルおよ

ウムの侵入とスパッタリングによる金属表面損耗の競争がび繊維状ナノ構造の形成過程の解明に取り組んだ．

問題となる．二体衝突近似法（BCA）を用いた解析の結果
８１２

Project Review

2.4 Peripheral Plasma Transport and Plasma-Wall Interaction

Y. Suzuki et al.

［５］，ヘリウムの場合には 100 eV 以下の入射エネルギーでた（図６参照）［７］．これにより，繊維状ナノ構造の成長過はスパッタリングが殆ど起きずに侵入だけが起こることが程を再現することに成功した．繊維状ナノ構造の成長速度解った．一方で 100 eV 以下の入射エネルギーで侵入できるを見ると，高さがヘリウム照射時間の平方根に比例するこ

深さは 10 nm 以下であり，ヘリウムバブルの分布に比べてとが知られているが，この点についても MD-MC ハイブ

もかなり浅いことが解った．

リッドシミュレーション［１１］と実験［３０，３１］で一致した．

第二段階となるタングステン内部のヘリウムの拡散・凝以上のように，材料の数値計算による評価としてヘリウ

集過程では，ヘリウムの凝集可能性について調べた．ヘリ

ムバブルと繊維状ナノ構造の形成過程を中心課題に据えた

ウム原子がタングステン中でバラバラに孤立している状態

が，この課題を進めることで，複数の計算手法の利用ならと，一か所に凝集している状態のエネルギーを比較するこ

びにそれらのハイブリッド手法の開発によるマルチスケーとで，どちらがエネルギー的に安定で存在しやすい状態でル解析を，PWI 現象においても実用レベルで導入することあるかを見出す．エネルギーの算出には金属内の電子状態ができた．これらの計算手法を，さらに発展させ，PWIを 量 子 力 学 的 に 求 め る 必 要 が あ り，密 度 汎 関 数 法

課題の解決の為に有効に活用できることが期待できる．特（DFT）に 基 づ い た 計 算 を OpenMX コ ー ド［６］を 用 い てに MD-MC ハイブリッド手法は繊維状ナノ構造の形成だけ行った．その結果，水素原子はタングステン中の空孔におでなく，従来の MD では計算時間の制限により実験に比べいても６個程度までしか凝集しないが，ヘリウムだけでなて 106 倍以上の照射フラックスを設定せざるを得ないといくネオン，アルゴンといった希ガス原子はタングステン中う問題点を解決し，実験同等の照射フラックスを扱えるにで際限なく凝集できる可能性が示された［４，７］．希ガス原ように工夫している．この手法は，この現象以外にも適応

子は閉殻電子構造の為に金属中の自由電子からは斥力のみできることが期待できる．

を受ける為，凝集してタングステン原子を押し広げた方が参考文献

エネルギー的に安定になると考えられる［８］．
第三段階となるナノバブルの成長過程では，扱うべき空［１］G. Kawamura et al., Contrib. Plasma Phys. 54, 437 (2014).
［２］G. Kawamura et al., Proc. 25th IAEA-FEC 2014, TH/P639.
［３］S. Takamura et al., Plasma Fusion Res. 1, 051 (2006).
［４］A.M. Ito et al., J. Nucl. Mater. 463, 109 (2015).
［５］S. Saito et al., J. Nucl. Mat. Suppl. 438, S895 (2013).
［６］OpenMX http://www.openmx-square.org/.
［７］A. Takayama et al., Jpn. J. Appl. Phys. 52, 01AL03 (2013).
［８］T. Tamura et al., Mater. Sci. Eng., 22, 015002 (2014).
［９］A.M. Ito et al., Phys. Scr. T159, 014062 (2014).
［１０］R. Kobayashi et al., J. Nucl. Mater. 463, 1071 (2015).
［１１］A.M. Ito et al., Nuclear Fusion 55, 073013 (2015).
［１２］M.J. Baldwin and R.P. Doerner, J. Nucl. Mater. 404, 165(2010).
［１３］Y. Noiri et al., J. Nucl. Mater. 463, 285 (2015).

間スケールが数ナノメートル以上に及ぶため DFT による計算は難しく，精度は落ちるが MD による古典力学計算を行う必要がある．しかし，プロジェクト発足時の段階ではタングステンとヘリウムの相互作用を精度良く扱えるポテ
ンシャルモデルが存在せず，MD 計算を実行できる状況ではなかった．そこで，Downfolding 法を用いた MD 用のポ
テンシャルモデルの開発を行った．１００原子程度からなる
サンプル構造を数百種類用意し，その構造に対するエネル
ギーを DFT と暫定ポテンシャルモデルで算出する．それらの差の二乗平均が最小になるように暫定ポテンシャルモ
デルのパラメータを修正していくことで，精度の良いポテ
ンシャルモデルを得ることができる．これにより，DFTとの差が 0.01 eV/atom 以下となる高精度なタングステン−ヘリウム系のポテンシャルモデルを得た［９］．本ポテン
シャルモデルを用いてヘリウムバブルの表面付近での破裂を調べた［９］．また名古屋工業大学：小林亮助教らは本研
究で開発したポテンシャルを用いた MD シミュレーション
により，バブル成長に伴うループパンチング機構を明らかにした［１０］．
第四段階となる繊維状ナノ構造の成長過程は，空間スケールで言えば MD の適用範囲内であるが，MD シミュ
レーションで扱える時間スケールはせいぜい１００ナノ秒程度である．その為，長時間スケールに及ぶタングステン内部のヘリウム拡散を MD では再現できず，浅い領域のヘリウム濃度が上昇しすぎて非現実的な破壊が起こってしまう．そこで，MD と動的モンテカルロ法（kMC）を組み合

図６

わせた MD-MC ハイブリドコードを開発し，ヘリウムバブ
ルからの圧力によるタン グ ス テ ン の 変 形 を MD で表現
し，ヘリウム原子の拡散をランダムウォークとして表現し

８１３

タングステンにヘリウムを照射した時の時間発展．紫がタングステン，水色がヘリウムを表す．タングステン表面に繊維状ナノ構造が生成される様子を，MD-MC ハイブリッ
ドシミュレーションを使って再現することができた．

J. Plasma Fusion Res. Vol.92, No.11 (2016)８
１４‐８２０

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

２．核融合プラズマシミュレーション研究
2. Simulation Researches in Fusion Plasmas

２．
５ 統合輸送解析コードの構築
2.5 Development of the Integrated Transport Analysis Code横山雅之
YOKOYAMA Masayuki
自然科学研究機構核融合科学研究所，総合研究大学院大学（原稿受付：２
０
１
６年７月１
３日）

大型ヘリカル装置（LHD）
をはじめとした環状磁場閉じ込めプラズマの物理機構解明と体系化に資するため，統合輸送解析コード体系 TASK3D の開発と運用を行っている．TASK3D は，トカマクを対象とした統合解析コード TASK［１］を基盤として，ヘリカルプラズマにおける３次元磁場配位に起因する特有の物理機構や，各種物理要素モジュールの３次元平衡データとの連携機能などを追加する形で開発を進めてきている．前回のプロ
ジェクトレビュー以降，「放電の全（多数）時間解析・予測」，「時間変化追随解析・予測」へと着実な進展を遂げている．これは，数値試験炉構築の基盤形成となる進展である．また，比較的簡便かつ短時間の要素コード（モジュール）統合の段階から，「プロセスされた実験データ提供を通じての大規模シミュレーションとの連携，および，それらの Verification & Validation 支援」の側面での位置づけも国内外において確立している．本章では，LHD 実験解析型，および，到達パラメータ予測型の開発・運用状況について報告する．
Keywords:
integrated transport code, TASK3D, LHD, experimental analysis, prediction

２．
５．
１ LHD 実験データ解析型 TASK3D-a の進展

データセットを蓄積でき，その系統的理解にもつながる機

LHD 実験適用型の Analysis Version TASK3D-a は，第一能拡張である．

版 a０１のリリース（平成２４年９月）から，第二版 a02 へと機ECH 加熱分布を評価する LHDGauss コード［３］および

能拡張を進めてきている．a01 と a02 を比較する形で，図２

TRAVIS コ ー ド［４］の 導 入 作 業 が 完 了 し た．LHDGaussに TASK3D-a02 のモジュール構成を示す．TASK3D-a02

コードは核融合研オリジナルのコードであり，最近，LHD

における新機能を以下に概略する．

解析データサーバとの連携確立（ジャイロトロンの情報：

新古典熱・粒 子 拡 散 フ ラ ッ ク ス を 評 価 す る GSRAKE

入射パワー，アンテナの向き，偏波の情報などを，実験

コードの導入作業が完了した．GSRAKE コード［２］はバウ

データベースから自動的に入手）によって，ほぼショット

ンス平均化法に基づく新古典輸送解析コードであり，LHD

間での ECH 加熱解析が可能となり，ECH 加熱の最適化にのように比較的単純な磁場リップル構造を有する磁場配位活用されている．一方，TRAVIS コードは，マックスプラ

に対しては，短時間に新古典拡散フラックス評価を行うこ

ンクプラズマ物理研究所のN. Marushchenko博士によってとができる．GSRAKE コードは，従来，単体で LHD 実験開発され，単体として LHD 実験解析に用いられてきてい解析に用いら れ てい た が，入 出 力 の 整 備 を 行 っ た上でた．TASK3D-a02 では，LHD 解析データサーバに登録され

TASK3D-a02 に 組 み 込 ん だ．新 古 典 粒 子 フ ラ ッ ク ス のている LHDGauss の解析結果の参照，および，TRAVIS

両 極 性 条 件 か ら，径 電 場 の 評 価 も 行 う こ と が で き る．

のその場実行を行い，電子温度分布計測がある全時間タイ

TASK3D-a01段階から可能であった実験的なエネルギーバ

ミングでの ECH 吸収パワー（TRAVIS では ECCD 電流評

ランス評価と，この新古典エネルギーフラックス評価が同価も）計算結果を取り込んだエネルギーバランス解析が行時に行われることになった．
（実験的なエネルギーフラッ

えるように組み上げられている．これにより，これまで，

クス）‐（新古典エネルギーフラックス）を「乱流輸送によるNBI加熱プラズマに限定されてきたTASK3D-aの適用範囲寄与」と考えると，TASK3D-a02 の実行によって，新古典が格段に拡がった．高イオン温度プラズマへのECH重畳に輸送と乱流輸送が径方向エネルギー輸送に占める割合のよる電子・イオン同時高温化プラズマの輸送解析などへのNational Institute for Fusion Sceience, SOKENDAI (The Graduate University for Advanced Studies), Toki, GIFU 509-5292, Japanauthor’s e-mail: yokoyama@lhd.nifs.ac.jp

８１４!2016 The Japan Society of Plasma
Science and Nuclear Fusion ResearchProject Review

図１

2.5 Development of the Integrated Transport Analysis Code

M. Yokoyama

統合輸送解析コード TASK3D 構築に関する取り組み全体の概要．LHD プロジェクト，核融合工学研究プロジェクト，および，京都大学との強力な連携によって大きな進展を遂げている．

適用がすでに始まっている．３次元平衡，計測温度・密度

している．アルフベン固有モードの解析や，それに伴う高分布などを揃えた形での LHDGauss，TRAVIS 間の系統的エネルギー粒子の挙動解析，新古典プラズマ流の解析ななベンチマーク作業も格段の進展が期待できる．

ど，現在，個別に計算や実験解析が行われている物理課題

また，国際的なヘリカル系実験の共同作業である，プラについても，TASK3D-a からプラズマの平衡や，プラズマ

ズマ閉じ込め時間や分布計測結果のデータベースに対しの温度・密度分布などを提供することで，連携を拡張させて，フォーマットを合わせた形でデータ提供する数値解析ている．これらの大規模シミュレーションコードはコンプログラムの組み込みを完了した（図２の左下部）．これ

ピュータ資源を多く必要とするものが多く，現時点では，により，LHD からのデータ供給がルーチン化することとな

そのままの形態で TASK3D-a に取り込むことは現実的でり，LHD の高性能プラズマのデータや，昨年から実験を開はないが，実験データの提供という「ゆるい連携」から今始した Wendelstein 7-X 装置（ドイツ：マックスプランク

後の展望を拓いていくことを検討している．これまでに，プラズマ物理研究所）の実験データなどを含んだ形で，ヘTASK3D-aからLHD平衡，計測分布データを提供した大規

リカル系プラズマ閉じ込めに関する新たなスケーリング則模シミュレーションコードの例を，以下にいくつか挙げの導出や，分布を用いたコードのベンチマーク作業などのる．
!FORTEC-3D コード［５］による国内外ヘリカル実験装

国際協力研究を活発に展開することができる．
TASK3D-a02における機能拡張を活きたものにするため

置での新古典熱フラックス評価の実験検証共同研究，には，上記のような，新古典両極性径電場，径方向エネル!GKV-X コード［６］による LHD 高イオン温度放電のイ

ギー輸送に対する乱流輸送寄与の系統データなど，実験計

オン温度勾配駆動微視的不安定性が引き起こすエネル

測との比較や多数データの取り扱いによる系統理解を進め

ギーフラックスの評価，!米国 MIT の M. Landremann 氏によって開発が進んでることが必要である．
さ ら に，TASK3D-a は，様 々 な 物 理 課 題 に 関 し て，いる新古典輸送解析コード SFINCS［７］による LHD「LHD 平衡・分布データ提供」という役割を通じて，大規プラズマの具体的計算（新古典輸送 解 析 コ ー ド の模シミュレーションコードと連携を進めるフェイズに到達

Verification 国際共同研究展開へのテスト ケ ー ス：
８１５

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

図２

TASK3D-a01（上部）と a02 での機能拡張．

TASK3D-aから供給したデータで計算が可能であるこ衡，同じ計測分布データが円滑に供給されることとなり，とを確認），

Verification & Validation（V&V）国際共同研究の加速に貢・アルフベン固有モード解析のための AE3D コード［８］，

献している．

アルフベン固有モード存在下での粒子軌道追跡計算も国際学会やワークショップ，CWGM（ヘリカル系調整作可能である DELTA5D［９］（いずれも，米国オークリッ

業会合）などでの継続的な研究発表を通じて，TASK3D-a

ジ国立研究所：D.A. Spong 氏との共同研究で LHD 解の国際的認知度も着実に高まっている．TASK3D-a の国際析が行われていたもの）．TASK3D-a で生成される結展開の例を以下に挙げる．

果などを用いて，AE3D，DELTA5D コードへの入力!N. Pablant 氏（米国：プリンストンプラズマ物理研究を自動的に生成するルーチンが作成された．両コード所）は，X 線イメージング結晶分光器（XICS）を LHD

とも TASK3D-a のサーバー群では計算が非常に困難に設置し，NBI 入射のない状況でのイオン温度やフであるが，計算実行のための入力ファイルを自動的に

ロー計測を可能とした．特に，ECH 入射による LHD

生成できるようになったことは，損失高速イオンプの高電子温度プラズマのエネルギー輸送特性評価のた

ローブ計測［１０］との比較に関して，解析数の増加によ

め，TASK3D-a01 を用いた輸送解析を行った．イオンる系統的知見の導出にも役立つ成果である．

温度分布計測データも取得できるため，そのようなプ

また，D.A. Spong 氏から，PENTA コード［１１］一式の提ラズマにおけるイオンのエネルギー輸送評価も行える供を受けた．これもサーバー群では計算が非常に困難であ特長があり，２０１３年度，２０１５年度の国際ステラレータ‐

るため，プラズマシミュレータで計算実行することになる

ヘリオトロンワークショップいずれにおいても本人が予定であるが，TASK3D-a の実行で，必要な平衡や入力

招待講演を行う［１２］など，TASK3D-a も含めた国際的ファイルが自動的に準備できるようになり，それらを用い評価の高い国際共同研究となっている．

て，LHD プラズマにおける新古典プラズマ流の磁気面平均!P. Vincenzi氏（イタリア：RFXコンソーシアム博士課

諸量を求めることができる段階となっている．さらに，磁程学生）が，ヨーロッパの EUROfusion 経費にて核融

気面上での分布を評価するルーチンを組み込むことで，荷合科学研究所に滞在（２０１４年，２０１５年）した．LHD

電交換分光計測の視線に沿った計測結果との比較を行うべの NBI グループと RFX コンソーシアムの共同研究のく，作業を進めている．

一環として，「NBI 加熱モジュールの高精度化と実験

新 古 典 プ ラ ズ マ 流 の 解 析 に 関 し て，FORTEC-3D，

データによる検証」に取り組むことが先方の指導教官

PENTA，SFINCS の３者に，TASK3D-a によって，同じ平と合意された．その合意を受けて，TASK3D-a におけ
８１６

Project Review

2.5 Development of the Integrated Transport Analysis Code

M. Yokoyama

る NBI 加熱評価の高精度化（ターゲットプラズマに複を TASK3D-a が有するようになったことは，今後の比数のイオン種が存在する場合へのルーチン的適用），

較研究にも大いに役立つ．この成果は，２０１４年のアメおよび，核融合反応による中性子発生量評価のための

リカ物理学会において，A.Y. Pankin 氏によって「Use

コード開発（反応断面積の最新の知見への更新）を

of Uncertainty Quantification Techniques for Interpre-

行った．LHD 第１
８サイクル実験において，水素，ヘリ

tive and Predictive Transport Analysis of Burning Plas-

ウム比を系統的に変えた高イオン温度プラズマの生成

mas」として口頭発表された［１４］．

が行われたが，それらを題材として，開発したコードTASK3D-aの適用によって，複数のLHD高イオン温度放の適用を行い，その成果について，２０１５年６月のヨー電の複数タイミング（合計およそ２
００タイミング）の解析

ロッパ物理会議，同１１月の第２５回国際土岐コンファレ

データベースを作成することができた．このデータベース

ンスにおいて本人が発表を行った．コード改良の詳細をに関して，「Modification of fit3d code for the analyses

乱流輸送の区分や，乱流輸送の原因となる揺動の種類など

of LHD deuterium experiments」としてまとめ，NIFS-

を全く意識することなく，径方向分布も含めたイオン熱輸

Memo として発刊した［１３］．

送係数の簡便な回帰式の導出に成功した．２
０１４年の IAEA

ビッグデータの視点で捉えることで，新古典輸送／

・２０１４年に，米国 Tech-X 社の A.Y. Pankin 博士が，核融核融合エネルギー会議におけるポスター発表などで成果発

合研客員教授（JIFT 客員）として，３か月滞在され表を行った［１５］．LHD 実験データベースと統合輸送解析

た．この滞在によって，実験観測分布の誤差範囲を考スイートとの連携によって初めて可能となった成果であ慮した熱輸送係数の値（計算結果）の幅をルーチン的り，将来の核融合炉におけるコアプラズマ制御における到に計算す る ツ ー ル（Uncertainty Quantification: UQ）

達温度予測の一手法として検討を進めていく予定である．

図３がTASK3D-aに導入された．LHDプラズマの熱輸送係

TASK3D-a による LHD プラズマ解析の主な具体例を紹数の多 数 ケ ー ス 解 析 に お い て，温 度・密 度 分 布 の介する．個別の解析に関する詳細や物理的知見についてのフィッティング形状などに結果がどのように依存する

説明は各文献に譲り，ここでは，TASK3D-a がどのようにかについて，系統的な知見を簡便に得ることができる

使われ，図が作成されているかについて概説する．紹介さようになった．実験解析と理論モデルとの比較などをれている図の一覧を図３に示す．

行っていく際に，双方の不確定性（データの幅）を意

１．参考文献［１６］（NUCLEAR FUSION 誌 Most Cited

識しておくことが必要であり，それを可能とする機能（from 2011）論文）の Fig. 2：プラズマ周辺部の径電場

TASK3D-a の解析結果を含む主な LHD 論文の図の一覧．TASK3D-a のどのモジュールの結果が用いられているかを示している．図の解釈等の詳細は各文献を参照されたい．

８１７

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

計測と新古典両極性径電場を，衝突周波数に対してプ成された図である．

ロットした図（Cal が計算値）．複数のケースに対し

６．参考文献［２１］（２０１４年 IAEA 核融合エネルギー会議のて，TASK3D-a を実行し，GSRAKE コードの結果を用口頭発表）の Fig. 4：高イオン温度プラズマのイオン

いて作成された図である．

温度勾配増進とイオン密度で規格化されたイオン熱フ

２．参考文献［１７］の Fig. 1：Co-，Ctr-NBI 入射プラズマに

ラックスの増え方（weak ITB と strong ITB の比較）のおけるトロイダル回転計測結果を，その他のプラズマ

時系列データを示した図．ビーム減衰も考慮した NBI

パラメータと共に示した図．NBI によって持ち込まれ吸収分布と，密度・温度の時間変化をも加味した熱フる外部トルクの図（最下枠）は，NBI 加熱計算の時系

ラックスの時系列データ（dytrans モジュールの結果）

列データ（FIT3D モジュールの結果）から作成されたから作成された図である．

図である［１７］．

２．
５．
２ 到達パラメータ予測型 TASK3D の進展

３．参考文献［１８］（九州大学との共同プレスリリース論文）の Fig. 1：プラズマ回転減速過程に お け るイオ

予測版については，京都大学との共同研究を通じて，ン・電子熱輸送係数の時間変化を示す図（最下枠）．

LHD の高イオン温度放電の熱輸送モデリング，時間変化を

ビーム減衰も考慮した NBI 吸収分布と，密度・温度の含めた予測シミュレーション，さらには，その実験検証が時間変化をも加味した熱フラックスから求めた熱輸送

進んでいる．粒子輸送や，水素ヘリウム比の変化に伴う達係数の時系列データ（dytrans モジュールの結果）か

成イオン温度の変化（LHD 第１
８サイクル実験で得られた系

ら作成された図である．

統的実験結果）に関する輸送モデリング，重水素プラズマ

４．参考文献［１９］（２０１２年 IAEA 核融合エネルギー会議のの到達温度予測などにも大きな進展があるが，ここでは，

口頭発表）の Fig. 8：高電子温度プラズマにおけるイ

現在の LHD 高イオン温度放電を対象とした輸送モデリン

オン，電子熱拡散係数を新古典の値と比較した図．過グに関して概説する．

渡的振る舞いがほとんどない状態であるので，trsnap

径方向の熱・粒子輸送拡散方程式を解く TASK/TR モ

モジュールと GSRAKE モジュールの結果から作成さジュールに，３次元平衡 VMEC コード，新古典拡散係数

れた図である．

データベースである DGN/LHD，両極性径電場を求めるモ

５．参考文献［２０］（２０１４年の EPS 招待講演）の Fig. 10：ジュール ER，さらに，ジャイロ平均化されたドリフト運動（a）電子，（b）イオン温度勾配に対する，電子密度で規論方程式を解く GNET
（その時間発展解析版 GNET-TD）

格化された（a）電子，（b）イオン熱フラックスの時間

コードを統合した体系での解析であ る（図４）．複数の変化を示す図．ビーム減衰も考慮した NBI 吸収分布

LHD 高イオン温度放電から，熱拡散係数として，電子はと，密度・温度の時間変化をも加味した熱フラックス

ジャイロボーム，イオンはジャイロボームに規格化イオンの時系列データ（dytrans モジュールの結果）から作

温度勾配をかけた形式（係数は，計測温度分布を最も再現

図４

京都大学で開発が進んでいる予測型統合輸送コードの計算フロー図．GNET（-TD）
が組み込まれている．輸送モデルを用いて，到達パ
ラメータ予測を行うように組み上げられている．

８１８

Project Review

2.5 Development of the Integrated Transport Analysis Code

するように決定されている）のものを利用している．LHD

M. Yokoyama

盤整備の観点からの研究連携である．

プラズマの高イオン温度状態は，炭素ペレット入射後の密核融合工学研究プロジェクトでは，LHD での密度・温度度減衰フェイズで実現されており，この状況に対応するた

分 布 を 直 接 外 挿 す る 手 法（Direct Profile Extrapolation:

め，複数イオン種，時間発展を加味した加熱分布計算を行DPE）［２３］によって，径方向の熱・粒子輸送拡散方程式をう GNET-TD コード［２２］を適用している．しかし，上記の解くコード体系の構築を開始していた．ペレット入射によ

イオン熱拡散係数での計算では，イオン温度の時間変化をる粒子供給に関して NGS モデルを適用して，ペレット入射全く再現することができなかった．そこで，炭素ペレット

LHD プラズマの密度・電子温度の時間変化を再現するな

入射（4.75 s）後に実効電荷数 Zeff が図５（右上図）のよう

ど，実験検証を行いながらの構築が進んでいた．このコーに時間変化することに着目して，Zeff に依存してイオン熱ド体系を FFHR-d1 の立ち上げシナリオ検討に適用するに拡散 係 数 が 低 減（乱 流 輸 送 の 抑 制 に 対 応）す る モ デ ル

当たり，新古典エネルギー輸送フラックスがプラズマのエ（図５中に記されている Zeff 依存モデル）を模索したとこ

ネルギーバランスを満たしているかどうかを定量的に評価

ろ，イオン温度の時間変化を比較的精度よく再現することするために，TASK3D-a のモジュール（VMEC コ ード，ができた．この熱輸送モデルを，図５の例よりは低いイオ

Boozer 座標への変換コード，さらに GSRAKE コード）をン温度（4 keV 程度）のプラズマ（炭素ペレット入射）に適移植した．これにより，密度・温度の時間発展に伴う新古用したところ，イオン温度分布の時間変化を再現できた．

典エネルギーフラックスの逐次計算が可能となり，エネル

従って，この輸送モデルを導出したプラズマ以外への適用

ギーバランス成立性を探る複数のシナリオ比較が容易に行性も確認が進んでいる．この輸送モデルの理論的根拠を探えるようになった（図６）．さらに，広範囲の密度・温度パると同時に，未踏イオン温度領域への予測計算への進展がラメータ範囲を有する POPCON（Plasma Operation Con-

期待されている．

tour）ダイアグラム上に，新古典エネルギーフラックスがまた，核融合研における到達パラメータ予測型統合計算

エネルギーバランスを満たしうる境界線を明示するなど，システムの構築も大きな進展をみせている．核融合工学研

単体の新古典輸送計算では困難である検討も行われ，その究プロジェクトが主体となって設計活動が進められている

成 果 も 含 む FFHR-d1 コ ア プ ラ ズ マ 設 計 の 最 新 の 知 見

LHD 型原型炉 FFHR-d1 設計グループと統合輸送シミュが，２０１５年の第２０回国際ステラレータ‐ヘリオトロンワー

レーショングループとの連携による，FFHR-d1 立ち上げシ

クショップにて，招待講演として発表された．このプロ

ナリオの定量解析である．まさに，数値試験炉としての基

ジェクト間連携をさらに進め，核融合研における予測型統

図５

LHD 高イオン温度プラズマにおけるイオン温度変化を再現するイオン熱輸送モデリングの概要．（右上）
炭素ペレット入射による高イオン温度放電でのペレット入射後の Zeff の時間変化．左上に記載されているイオン熱輸送係数を用いた計算を行ったところ，イオン温度の時間変化（左下：実験計測）をほぼ再現する結果（右下）を得ることができた．

８１９

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

図６

FFHR-d1 立ち上げシナリオ検討に当たり，TASK3D-a のモジュールを移植し，逐次，新古典エネルギーフラックスの計算が行える環境を整備した．このことにより，時間発展を伴う立ち上げ時のエネルギーバランス成立性のシナリオ比較などが容易に行えるようになった．

［１２］N. Pablant et al., Plasma Phys. Control. Fusion 58, 045004(2016).
［１３］P. Vincenzi et al., NIFS-MEMO-74, Nov. 19 (2015).
［１４］A. Pankin et al., Proc. of 56th Annual Meeting of the APS Division of Plasma Physics, 59 N15 (Oct. 27-31, 2014, New Orleans, Louisiana) JO3.00009.
［１５］M. Yokoyama, Plasma Fusion Res. 9, 1302137 (2014).
［１６］K. Kamiya et al., Nucl. Fusion 53, 013003 (2013).
［１７］K. Ida et al., Phys. Rev. Lett. 111, 055001 (2013).
８］K. Ida et al., Nature Communications, 08 January (2015).
［１［１９］H. Takahashi et al., Phys. Plasmas 21, 061506 (2014).
［２０］K. Ida et al., Plasma Phys. Control. Fusion 57, 014036(2015).
［２１］K. Nagaoka et al., Nucl. Fusion 55, 113020 (2015).
［２２］H. Yamaguchi et al., Plasma Fusion Res. 9, 3403127 (2014).
［２３］J. Miyazawa et al., Fusion Eng. Des. 86, 2879 (2011).

合輸送コード（数値試験炉の基盤）として整備・拡張していく検討も行っている．

参考文献［１］http://bpsi.nucleng.kyoto-u.ac.jp/task/［２］C.D. Beidler and W.D. D'haeseleer, Plasma Phys. Control.
Fusion 37, 463 (1995).
［３］T. Ii Tsujimura et al., Nucl. Fusion 55, 123019 (2015).
［４］N. Marschenko et al., Plasma Fusion Res. 2, S1129 (2007).
［５］S. Satake et al., Plasma Fusion Res. 3, S1062 (2008).
［６］M. Nunami et al., Plasma Fusion Res. 5, 016 (2010).
［７］M. Landremann et al., Phys. Plasmas 21, 042503 (2014).
［８］D.A. Spong et al., Phys. Plasmas 17, 022106 (2010).
［９］D.A. Spong et al., Bull. Am. Phys. Soc. 44, 215 (1999).
［１０］K. Ogawa et al., Nucl. Fusion 52, 094013 (2012).
［１
１］D.A. Spong, Phys. Plasmas 12, 056114 (2005).

８２０

J. Plasma Fusion Res. Vol.92, No.11 (2016)８２１‐８３０

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

３．計算科学の高度化と基礎物理シミュレーション
3. Sophistication of Computational Science and Fundamental Physics Simulations
石 黒 静 児，伊 藤 淳，宇 佐 見 俊 介，大 谷 寛 明，坂 上 仁 志，
樋 田 美 栄 子，長 谷 川 裕 記，堀 内 利 得，三 浦 英 昭
ISHIGURO Seiji, ITO Atsushi, USAMI Shunsuke, OHTANI Hiroaki, SAKAGAMI Hitoshi,
TOIDA Mieko, HASEGAWA Hiroki, HORIUCHI Ritoku and MIURA Hideaki
自然科学研究機構核融合科学研究所（原稿受付：２
０
１
６年７月１
３日）

バーチャルリアリティ（VR）可視化，多階層シミュレーション手法などの計算科学の高度化，ブロブの運動論的挙動，波と粒子の相互作用に伴う粒子加速などの基礎プラズマ過程，拡張 MHD 理論及びシミュレーション，
レーザー核融合の関連研究など，中長期的視野で展開している幅広い研究について紹介する．
Keywords:
simulation, VR, simulation science, multi-hierarchy, blob, extended MHD, laser plasma, shock, turbulence

数値実験炉研究プロジェクトでは，特定の装置の特定のを崩壊させていると考えられている．そのため，ダスト粒

物理現象に着目した核融合シミュレーション研究に加え子の生成・輸送機構を調べるため，LHD 実験では高速カメて，シミュレーションによる予測性を高めたり，その適用

ラを２台取り付けて立体画像を撮影することで，ダスト粒範囲を拡大する先進的なシミュレーション手法に関する研子の３次元座標の時系列データ取得に成功した［３］．この究，シミュレーション研究の基盤となる可視化，並列計算

ダスト粒子軌道と磁力線の解析を２次元平面で行う場合，技術などの計算科学高度化のための研究，特定の装置に限

奥行き情報が失われているため，その立体的構造を把握す

定されない基礎的な物理現象の解明をめざした研究を包含ることが困難である．そこで，VR空間でそれらの立体構造している．本章ではこれらの研究について概観する．

を調べるため，本手法の開発に着手した．実験観測データ
をデータ変換を経ずに読み込むためのインターフェースの３．
１ バーチャルリアリティ空間における同時可視化

開発，大量のダスト粒子をスムーズに VR 空間に表示するための可視化機能の開発を行った．また，その軌道の観測

バーチャルリアリティ（VR）
装置"CompleXcope"は，観からダスト粒子の輸送方向を把握するための工夫も行っ

測者を VR 空間に投入させ，あたかも目の前に物体が存在た．その VR 可視化の例を図１に示す．その解析の結果，ダしているかのように感じさせることができる装置である．

スト粒子は周辺領域に存在していること，多くのダスト粒

また，観測者の動きに合わせてその物体が動いたり，コ子が磁力線に沿って輸送されている一方，いくつかのダスントローラーで動かしたりすることができるので，プラズ
マ現象を中から見たり外から見たり，いろいろな角度から観測することができる．まさに，コンピュータが作る様々な VR の世界に観測者を入り込ませ，物理現象を「目の当りに」することで，人間の空間を認識する能力を総動員して観測することができる．このような VR 装置の特徴を生かして，シミュレーションデータと実験装置データを VR
空間で同時に可視化する手法を活用してきた［１］．
この可視化手法の研究をさらに発展させて実験観測デー
タを同時に可視化する手法を開発した［２］．
図１

大型ヘリカル装置（LHD）の長時間放電実験では炭素放射の急増とともにプラズマの閉じ込めが壊れることが観測されている．この急増はダイバータ領域から大量のダストが放出される時点と同期しており，ダスト粒子が閉じ込め

ダスト粒子データとシミュレーション，実験装置データの同時 VR 可視化．左図は観察している様子を表し，右図は3-O ポートから観測されたダスト粒子データをトロイダル
方向から見た図．緑線は磁力線，マゼンタ色はプラズマの等圧面を表している．

National Institute for Fusion Science, Toki, GIFU 509-5292, Japancorresponding author’s e-mail: ishiguro.seiji@nifs.ac.jp

８２１!2016 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

ト粒子が磁力線を横切るように輸送されていること，その存在領域によって輸送の方向が異なっていることを容易に観測することができた．このように，ダスト粒子の分布する領域や軌道と磁力線の比較，その立体的な構造の理解を，訓練を積んだ研究者だけではなく，初学者にも容易に行うことができるようになった．
本研究は，兵庫県立大学：大野暢亮教授，神戸大学：
陰山 聡教授，甲南大学：田村祐一教授との共同研究であ
る．
図２

３．
２ 無衝突磁気リコネクションの多階層シミュ
レーション

多階層モデルによる磁気リコネクション結果．（a）
では単一の X 点，（b）
では複数の X 点を持つ磁気リコネクションが駆動されている．

磁気リコネクションは，核融合をはじめとした様々なプラズマで見られる普遍的なエネルギー解放過程であり，ま関する国際シンポジウム（ISSS）といった国際会議におけた典型的な階層横断現象でもある．リコネクションが起こる招待講演として選ばれた．

ると磁場のトポロジーが系全体に及ぶスケールで変化し，

３．
３ プラズマコヒーレント構造の微視的ダイナ
ミクスのシミュレーション

大規模な輸送も発生する．その一方，リコネクションが駆動されるために必要な電気抵抗は，運動論的な過程が基になって生成される．このマクロとミクロの過程は独立して近年，磁場閉じ込めプラズマ装置の実験において，そのいるのではなく，複雑に絡み合って発展している．我々は，装置形状によらず，周辺領域で間歇的に発生するフィラメ

階層横断現象として磁気リコネクションの全体像をつかむ

ント状のプラズマコヒーレント構造が観測されている．こことをめざし，マクロとミクロの物理を同時にかつ自己無の構造は，プラズマブロブとも呼ばれ，従来の拡散的描像

撞着に解く多階層シミュレーションモデルの開発を行ってでは説明できない動径方向のプラズマ輸送を担っているといる．我々の多階層モデルでは，シミュレーションの実空考えられている．このブロブ現象に関するこれまでの理論間を分割するという大きな特色がある．すなわち，領域を的・数値的研究では，その多くで，流体的な描像に基づく

マクロ階層とミクロ階層に分け，マクロ階層の物理は磁気２次元簡約化モデルを用いた議論がなされてきたが，その流体（MHD），ミクロ階層は粒子（PIC）シミュレーション

ような巨視的モデルでは，磁力線方向の電流を評価する

により計算する．また，２つの階層間には互いのデータを際，固体壁とプラズマとの間の電位構造やイオンと電子のスムーズに交換するため，有限の幅を持ったインターフェ

速度差の効果などを簡易的なモデルで仮定するなどしておイス領域を挿入する．この手法により，物理的根拠が曖昧

り，その妥当性の評価が問題となっている．そこで，上記な仮定・モデルを用いず，第一原理計算に基づいたミクロのような粒子的効果がプラズマコヒーレント構造の挙動にの物理をマクロ階層に取り入れることができる．

与える影響を明らかにし，今後のより精緻な多階層連結型

この多階層モデルを用いて，マクロのダイナミクスが磁シミュレーションの開発［５］，そして，周辺プラズマ輸送

気リコネクションのミクロな物理に及ぼす影響について初の予測に資することをめざ し，３次元静電粒子コードの解析を行った［４］．我々は，MHD 領域（マクロ階層）か（p3bd コード）によるプラズマコヒーレント構造ダイナミら流入させるプラズマインフローのパターンが異なる多階

クスの研究を行っている［６］．これまでの p3bd コードに層シミュレーションを実行した．インフローは上流境界のよる計算では，図３（a）のように，構造内に自己無撞着に形

中央（#!!
）から流入し始め，流入する領域は速度 "+ で #成されるらせん状の電流系が観測された［７］．このことは，

方向に広がっていくものとした．図２は時刻 $)*!!"'&!
!
"(，（b）は%における 磁 力 線 を示 し て お り，（a）は "+ !!
"+ !$!
"(（"( はアルフベン速度）の場合の結果である．
!
（a）においては，PIC 領域（ミクロ階層）で駆動されたリコ
ネクションは単一の X 点を持っており，定常なリコネク
ションに緩和していくと見られる．一方（b）では，複数のX 点を持ち，間欠的なリコネクションが継続する．多階層シミュレーションにより詳しく調べた結果，"+ ""( の場合!
"( では複数の X 点を持つリコネク!
は単一の X 点，"+ ##
図３ （a）
ブロブ構造内部に自己無撞着に形成される電流系（左手前側の流線．右奥側の等値面は電子密度分布を示している．
）と（b）
ブロブ内部の高ポテンシャル側，および，低ポ
テンシャル側における磁力線方向の電子速度分布（太線が高ポテンシャル側，細線が低ポテンシャル側を示している）．

ションが駆動される傾向があることがわかった．
この多階層シミュレーションによる成果は，「核融合プ
ラ ズ マ 理 論」に 関 す る Varenna-Lausanne 国 際 ワ ー ク
ショップ，アジア太平洋プラズマ理論会議（APPTC），計算物理国際会議（CCP），「宇宙空間シミュレーション」に
８２２

Project Review

3. Sophistication of Computational Science and Fundamental Physics Simulations

S. Ishiguro et al.

粒子シミュレーションを用いることにより，ブロブ構造における電流クロージャの問題を解決できる可能性を示している．また，ブロブ内部に形成される温度構造も明らかにされた［７］．ブロブの内部には，その断面において，ダイ
ポール状の電位構造が形成されるが，その高ポテンシャル
側，および，低ポテンシャル側における磁力線方向の電子速度分布を調べたところ，図３（b）に示したように，高ポ
テンシャル側の電子温度のほうが高いことがわかった．

３．
４ LHD の Hall MHD シミュレーション
図４

LHD における短波長圧力駆動型（交換型，バルーニング）不安定性の飽和に関するシミュレーション研究を行っ

LHD 中のバルーニングモードの３次元シミュレーションによる，圧力分布．

ている．バルーニングモードに関するこれまでの完全３次元圧縮性 MHD シミュレーション研究から，以下の点がわ

ロイダル波数 !!"!程度）以上の波数の不安定モードを打かっている．

ち切るような数値フィルターを作用させた．この結果，数１．粘性が十分に大きい場合にはバルーニングモードの成値シミュレーションはバルーニングモードの飽和にいたったが，その飽和プロファイルにおける !値は，Hall 効果な長は飽和し，かつ，良好な磁気面をもつ状態に自律的しで磁力線平行方向熱伝導効果を入れた場合よりも低く，に回復する［８］．

また，密度分布もホロー分布になるなど，むしろ不安定効

２．粘性を，実験で想定されるような小さな値にするほど果に作用する結果になった．このシミュレーションからわにバルーニングモードの成長は増大し，プラズマコア

かったことは以下の点である．

を大きく壊すまで，その成長は飽和しない．飽和した段階では，密度分布が往々にして初期の圧力勾配とは１．短波長成分の打ち切りは，飽和プロファイルに対して逆方向に勾配をもつホロー分布になるなど，通常想定

悪影響を与える．これは，短波長成分の打ち切りが長されるような飽和状態とはまったく異なった状態（崩波長成分との相互作用を失わせ，結果として長波長成壊状態）となる．また，良好な磁気面への自律的な回分の線形成長段階を延長させる作用を及ぼすためであ復は今のところ見られない．

る．

３．この崩壊過程において，圧力分布の平坦化，圧縮性に２．Hall 効果は不安定性を強化しえる．これは，Raleigh-よる成長率の低下などは発生している（あるいは考慮

Taylor 型不安定性の２次元拡張 MHD シミュレーショ

されている）が，飽和に対して十分な影響を及ぼしてンで得られた結果と同様に，２流体効果が平衡プロいるようには見られない．

ファイルの性質によって，不安定性の強化にも抑制に

４．３次元シミュレーションにおいて，磁力線平行方向へも作用し得ることと関連していると考えられる．

の熱伝導効果を高精度で実装することにより，この崩

２次元 Rayleigh-Taylor 不安定性のシミュレーションから壊をある程度緩和することがわかった．また，平行方

も，LHD の３次元短波長不安定性の抑制には，２流体効果

向熱伝導の影響により，平行方向に対して熱エネル

だけでなく有限ラーモア半径効果も同時に必要であると考ギーが運動エネルギーに変換され，結果として運動エえられる．また，２流体効果も有限ラーモア半径効果は短

ネルギーは大きくなりえることがわかった［９］．

波長成分も，磁場および速度場の短波長成分を著しく励起

以上の研究を踏まえて，LHD の短波長圧力駆動型不安定性

することがわかっているため，適度にこれらを打ち切りつについて，Hall 効果などを取り入れた拡張 MHD シミュ

つ，非線形段階においては，短波長成分が非線形項を通じレーション研究をはじめた．この研究は，LHD での３次元て長波長成分に及ぼす影響を適宜代替し得るようなモデル

シミュレーションだけでは困難が予想されるため，２次元の開発を行う必要がある．このような考え方は，流体工学

矩形形状での Rayleigh-Taylor 型不安定性の拡張 MHD シ

において Large Eddy Simulation（LES）と呼ばれる．Hall

ミュレーションによる短波長成分生成の基礎研究，Hall効果による LES の最初の取り組みは，国際会議において発効果が短波長成分の励起と空間構造の変化に及ぼす影響に表され［１１］（招待講演），その後の進展を含めて論文［１２］ついての大規模一様 Hall MHD 乱流シミュレーションを平として出版された．

行して行っている．これら２件の研究については他の節を３．
５ MHD 平衡・安定性理論の拡張のための数値解析参照されたい．
LHD での３次元シミュレーションに Hall 効果を入れた・強いポロイダル流をもつ平衡に対する磁気座標

シ ミュレ ー シ ョ ン を 行 っ た［１０］．Hall 項 は そ の 性 質 から，きわめて短波長の分散性波動を生成し，時間および空

強いポロイダル流をもつ有限ベータトロイダル平衡の安間について高い数値解像度を要求する．このため，比較的定性解析に必要な磁気座標系を構築した．流れをもつ軸対低波数のモードに焦点を当てることを念頭に，中程度（ト称平衡を表す一般化された Grad-Shafranov 方程式を逆ア
８２３

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

スペクト比で展開することで，ポロイダル音速程度の流れの両方を含む場合（FLR+TF）においては，有限ラーモアをもつ高ベータトカマク平衡に対する方程式系が得られ半径効果のみの場合とは異なるベータ値で完全安定化の消［１３］，円形断面の固定境界と線形の磁気面量を仮定するこ失が現れ，どちらか一方を含む場合と比べ最も強い安定化とで磁束関数の解析解が得られる［１４］．１次の磁束関数のを示さない場合があることを明らかにした．また，固有関

解は流れのない場合と等しく，２次の解に流れの効果が含数と局所解析における成長率の空間依存性とを比較し，拡まれる．この解析解を用いた座標変換式を数値的に解き，張 MHD シミュレーションで見られた有限ラーモア半径効

ポロイダル断面における磁気面量 +とポロイダル角$によ果のみを含む場合に短波長で成長率が局所解析に漸近しなる磁気座標 %+!
$&を構築した．図５は各磁気面における圧い現象［１５］が，有限ラーモア半径効果による固有関数の非力のポロイダル分布を示したもので，亜音速流の場合（a）

局在化と，局所成長率が周辺領域で中心よりも大きくなっにはトーラス外側（$#"付近）で圧力が高くなり，超音速ていることの組み合わせに起因していることを示した．

流の場合（b）にはトーラス内側で（$#$)付近）で圧力が数値的に求まり，径方向分布のポロイダルマッハ数に対す

３．
６ ２次元 RT 不安定性の拡張 MHD シミュレーション

る依存性を調べた結果，亜音速領域では超音速領域に比２次元矩形領域中の Rayleigh-Taylor 不安定性の線形・高くなる．この座標系を用いることで圧力の磁気面平均がべ，より中心にピークした分布になることが明らかとなっ非線形不安定性に対する２流体および有限 Larmor 半径

た．また，座標変換式を逆アスペクト比で展開することで（FLR）効果を数値シミュレーションによって調べた．この

解析的に磁気座標が得られ，ポロイダル流速が小さい場合

際の注目点は，低 %領域における，高波数成分のダイナミに数値的に得られた座標と良い一致を示した．

クスおよびこれに伴う空間構造形成である．本研究では，・有限ベータプラズマにおける Rayleigh-Taylor 不安定性

Roberts と Taylor（１９６２）［１６］，P. Zhu（２００８）［１７］らによっに対する二流体及び有限ラーモア半径効果のパラメータて用いられたものと同様の，２流体効果およびジャイロ粘依存性性項（FLR 項）を含む拡張 MHD 方程式を数値的に解くこ

非一様な磁場をもつ有限ベータプラズマ平衡におけるとで，上記の問題に取り組んだ．この方程式は以下のような形をしている．

Rayleigh-Taylor（RT）不安定性に対する二流体及び有限ラーモア半径効果のパラメータ依存性を，拡張 MHD モデ,*, *#! %
，
&&
,
/
,
,
&

ルを用いた線形解析により調べた．短波長モードに対する局所解析の結果を図６に示す．有限ラーモア半径効果によ,, *#*-'-&!#
% &
#! #
'
&
,
/
,
, &
&

り臨界波数以上で完全に安定化する現象が，非一様な磁場で実現する圧力勾配の小さい平衡では低ベータで見られな

##
%
.
&'
"!+" ( ("&&'"&#&'$
!*%
!
$，
$

くなる．二流体（TF）
効果と有限ラーモア半径（FLR）効果

,
!%
%
.
, !
-'"-(&#('(#! '% %"+&
,
/
,
,
&
%
-$，
".('
)#
*!*
(
)*"
+
!% # #*-&-&"
，
(!#
$
,
#
, -#
&
#! %
% ' &!#
'&&
,
/
,
,
&
!'
&
'
(

+$
, ' ' "# ,．
! ()* ) *!
"&
,
/
,
/
(
'*

ここで #
#
,
/
&はそれぞれ磁場ベクト
("
&，"
&#'
&
'
(,
'および 図５ル，電流密度ベクトル，速度ベクトルの第&成分である．記各磁気面（+= const.）における圧力のポロイダル角 $ 方向
分布．（a）
亜音速流（b）
超音速流の場合．

図６

号 &&'，'，*および +はそれぞれ Kronecker のデルタ，&'(，(

成長率の波数依存性の流体モデル間比較（a）%= 0.01（b）%= 0.1（c）%= 0.3．

８２４

Project Review

3. Sophistication of Computational Science and Fundamental Physics Simulations

S. Ishiguro et al.

Levi-Civita の反対称テンソル，比熱比，質量密度，および!
$成分
圧力である．記号 ""# はジャイロ粘性テンソルの第"!
#である．式中で繰り返される添え字については１および２の和をとる．数値シミュレーションは２次元矩形領域を直交格子上の４次精度中心差分法と４次精度の Runge-KuttaGill 法を用いて行った．
この研究については，２本の論文が発表された［１５，１８］．
最初の論文［１８］では，二流体効果，FLR 効果を独立したパラメータと見なし，それぞれの効果の有無に応じて不安定性の成長がどのように変化するかを調べたものである．いくつかのパラメータの組み合わせから，一流体効果，FLR
効果はそれぞれ単独では線形成長に与える影響が小さいのに対して，両者が加わった場合には，線形成長率を著しく低減する効果を発揮することがわかった．また，この両者が加わることで線形段階の成長率は低下するものの，非線形成長段階における飽和レベルには大きな影響を与えないことが示された．さらに，混合層幅の成長（これは非線形成長の指標の一つでもある）においても，二流体効果とFLR 効果の組み合わせは初期成長を抑止するものの，最終的には一流体 MHD の場合と同程度の大きさまで成長する事が示された．
第二の論文［１５］は，第一の論文とは異なる平衡についての研究が行われた．#値１％程度の領域において，二流体効果および FLR 効果の２つの効果が，高波数の RayleighTaylor 不安定モードを安定化することが，数値シミュレー
ションおよび線形（WKB）解析によって示された．二流体効果，FLR 効果はそれぞれ単独では十分な安定化効果を発揮しないが，両者が加わった場合の線形解析では，有限波数において成長率がゼロとなる．非線形シミュレーション
では，線形段階における成長率が WKB 解析と良好な一致

図７

Rayleigh-Taylor 不安定性の２次元拡張 MHD シミュレーション．（a）
線形段階初期，（b）
十分に発達した段階．

図８

Rayleigh-Taylor 不安定性の２次元拡張 MHD シミュレー
ション．二流体項およびFLR項を伴う拡張MHDシミュレー
ションによって得られた密度等高線（下段）．上段のパネ
ルは水平（x）
方向の密度，速度場の変化を表す．

を見せた．また，二流体効果と FLR 効果の両方を取り入れらたシミュレーションでは，線形段階において低波数モードのみが有限の成長率をもつため，Rayleigh-Taylor不安定性特有のマッシュルーム構造が大規模に成長する一方，
FLR 効果のために密度界面が急峻化し，これに合わせて密度界面に沿って急峻な速度せん断層が形成される．この結果，密度線断層に沿った Kelvin-Helmholtz 的二次不安定性が発生 し，短 波 長 成分 は 二 流 体 効 果，FLR 効 果 の な い
MHD モデルによる Rayleigh-Taylor 不安定性に比べて高波数成分が卓越することがわかった．図７には論文［１８］より，一流体 MHD シミュレーションにおける密度場の時間変化（a）線形段階初期，（b）十分に発達した段階を，図８には，拡張 MHD シミュレーションにおける密度場（下段）及び速度場水平方向成分（上段）を表す図を示す．なお，この論文は Phys. Plasmas 第２２巻３号において"Featured article''に選ばれるとともに，論文中の図が表紙を飾った．

３．
７ レーザープラズマのシミュレーション
コーン付きターゲットを用いた高速点火レーザー核融合
で重要な役割を果たす高速電子は，大きな発散角を持っており，かつエネルギーが高すぎるためコア加熱効率の著しい低下が懸念されている［１９］．そこで，高エネルギーイオ
８２５Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

ンを用いて追加熱を行う手法について統合シミュレーショ

効率の原因を調べると，Au コーンの側壁と CH 薄膜が接地ンを行い，イオンの寄与および加熱効率を評価した．イオ

しているため，シース電場が生成されると Au コーンからン加速の機構は，主に二つある．ターゲット裏面におけるCH 薄膜に背景電子が流れ込み，シース電場を弱めるため

シース電場によるイオン加速［２０，２１］では，コーンチップ

であることがわかった．そこで，Au コーンと CH 薄膜の間の前方に薄膜を置きイオン源とする．また，ターゲット表に 0.5 μm の真空ギャップを設けて，同様のシミュレーショ

面におけるポンドラモーティブ力を使うイオン加速［２２］ンを行った．その結果を図９に赤色で示す．イオンエネル

では，コーンチップ前面に低密度のフォームを置きイオン

ギーには大きな差はないが，+'! の最大イオンビーム強度源とする．この両機構を用いたイオン加速について，２次は２倍となり，変換効率も２％になった．

元粒子コードを用いたシミュレーションによりイオンビー

イオンを超高強度短パルスレーザーで加速する別の機構

ムの特性を明らかにすると共に，イオン源ターゲットの構として，ポンドラモーティブ力による直接加速がある．こ造を工夫して，加熱効率の向上をめざした［２３‐２５］．

の機構でも，シース電場による加速機構と同様にイオン源

イオンを超高強度短パルスレーザーで加速する機構にのターゲット材を CH とし，コーンチップを低密度 CH

は，先行する高速電子によって生成されるイオン源となるフォームで 8 μm コーティングした．このターゲットの電ターゲット裏面に垂直なシース電場を用いるものがある．

子密度を図１０（a）に示す．ポンドラモーティブ力により加＋イオン源として不純物から生成される H イオンを用いる

速されるイオンのエネルギーは，レーザー強度と平均電離と，加速されるイオン量等を制御できないので，イオン源度に比例し，フォームの電子密度に反比例する．レーザー＋のターゲット材を CH とした．このため，H イオンだけで強度は外部要因として一定であると考え，平均電離度も完はなく +'! イオンも補助加熱に利用できる．また，イオン

全電離を仮定すると一定となるので，イオンエネルギーをエネルギーの違いによるコア到達時刻の差を小さくするた

大きくするためには，フォームの電子密度を小さくしなけめ，そのイオン源ターゲットをコーン内のできるだけコアればならない．しかし，フォームの電子密度を相対論的な近傍に設置することを考えた．そこで，チップ幅10 μm，開臨界密度以下にすると，レーザーはフォーム中に異常侵入#-2，!""*(，""$!）
口角６０度のコーンを Au プラズマ（#!

するため，ポンドラモーティブ力によるイオン加速は抑制

!
#-2，
(
"%
とし，チップの 5 μm 手前に 4.5 μm の CH 薄膜（"されてしまう．このため，ポンドラモーティブ力による最!
#-2，!""，"""）を置いた．つまり，)'!""#，""'と #大イオンエネルギーは，フォームの電子密度を相対論的なコーンチップと CH 薄膜の間に 5 μm の真空領域があり，こ臨界密度にしたときに得られる．相対論的な臨界密度は，

こにシース電場が生成される．このターゲットの電子密度

レーザー強度の 1/2 乗に比例するので，結局，イオンの最を 図９（a）に 示 す．そ し て，加 熱 レ ー ザ ー と し て 強 度

大エネルギーは，レーザー強度の 1/2 乗と平均電離度に比10 W/cm ，波長 1.06 μm，半値全幅 10 μm のスーパーガウ

例してスケールされることになる．そこで，加熱レーザー

を照射
シアンビーム（""&，#
2
0
3
."#
/
,
1
1"50 fs，#/
1,
4"400 fs）の強度を 5×1019，1×1020，5×1020 W/cm2 とし，それぞれし，コーンチップの後方 1 μm において高速電子および加のレーザー強度のときの CH フォームの電子密度が相対論速されたイオンを観測した．

的臨界密度，つまり 6.12，8.60，19.13 #-2になるよう H＋と

202

イオンビーム強度の時間発展を図９（b）に，時間平均し+'! の密度を決定した．そして，加熱用レーザーを照射し，たイオンエネルギー分布を図９（c）に青色 で 示 す．ここ同様に高速電子および加速されたイオンを観測し，前述ので，実線は H＋，破線は +'! である．最大イオンビーム強度

スケーリングについて調べた．なお，レーザーのパラメー'!

は，H で 4×10 W/cm 程度，+ で 10 W/cm 程度であ
＋

17

2

18

2

タは，強度以外はすべて同じである．

り，加熱用レーザーの強度が 1020 W/cm2 であることから，

イオンビーム強度の時間発展を図１０（b）に，時間平均しレーザーからイオンビームへの変換効率は，１％以下と低たイオンエネルギー分布を 図１０（c）に示す．ここで，赤いことがわかる．また，H＋は 10 MeV 程度，+'! は 50 MeV色，青 色，緑 色 は，そ れ ぞ れ レ ー ザ ー 強 度 が 5×1019，

程度のエネルギーまで加速されたこともわかる．低い変換

1×1020，5×1020 W/cm2 であり，実線は H＋，破線は +'!

（a）ターゲットの電子密度（b）
イオンビーム強度の時間発展
図９

シース電場によるイオン加速．

８２６

（c）
時間平均したイオンのエネルギー分布

Project Review

3. Sophistication of Computational Science and Fundamental Physics Simulations

S. Ishiguro et al.

である．イオンの最大エネルギーは，概ねスケーリング通撃波は，衝撃波中に形成される強い電磁場を介した様々なりになっており，イオンビーム強度はレーザー強度に，ほ機構で粒子を強く加速できることが明らかにされてきたぼ比例している．

［２７］．例えば，イオンは衝撃波面での電磁場の急激な立ち

２次元粒子コードのシミュレーションで得られた高速電上がりにより，一部が反射され加速される．加速イオンは，子および高エネルギーイオンを１次元フォッカープランク

電子のジャイロ周波数 ##がプラズマ周波数%&#より大きい

コードに対する時間依存性がある粒子源のデータとして，

とき，相対論的エネルギーになり得る．そして，この機構

統合シミュレーションを行った［２６］．その結果を表１に示す．

によって太陽フレアにおけるイオンの急激な加速を説明できると提案された．

シース電場によるイオン加速の場合，真空ギャップがな

一方，電子は磁場に対して斜めに伝播する衝撃波によっ

いとイオンによる補助加熱の平均コア電子温度の上昇率は２７％に留まったが，真空ギャップがあると１００％，つまり，て，一部が捕捉され加速される．##!
%&#"!で，衝撃波の伝平均コア電子温度は高速電子だけの場合に比べて２倍と（ここで，!は光速，$は伝播角）に近
播速度 "'$が !"%'$なった．しかし，イオンを加速するシース電場は高速電子

いとき，加速電子のエネルギーは超相対論的となる［２８］．

流を阻害するため，真空ギャップがある場合の高速電子

また，捕捉電子は，衝撃波の主パルス付近の電磁場に大き

ビーム強度は低下し，平均コア電子温度の最大温度は逆にな影響を与える．空間１次元の場合は，磁場に平行方向の悪くなった．このため，高速電子と高エネルギーイオンは電場が強くなり，電子の捕捉が安定化する．それに対し，

トレードオフの関係にあり，イオン源の構造や素材につい

空間２次元のシミュレーションでは，捕捉電子はホイッスて，最適設計する必要がある．

ラー波を励起し，その非線形発展によって主パルス中の多一方，ポンドラモーティブ力によるイオン加速の場合，次元電磁擾乱が大振幅となる［２９］．そして，この多次元電レーザー強度が大きくなる程，高速電子のエネルギーが高磁擾乱により，捕捉電子は主パルスから解放され，その一くなりすぎてコア加熱には適さなくなるため，イオンによ部は旋回運動に伴って衝撃波を出入りすることにより，更る補助加熱の割合が大きくなることがわかった．しかし，に高いエネルギーになること示された［３０，３１］．

平均コア電子温度の最大温度は，投入したレーザーエネル

このように捕捉電子が引き起こす電磁場の変化は，電子ギーには比例しておらず，この場合も，レーザー強度と生の運動に重要な影響を及ぼすことが明らかにされたが，捕成される高速電子おおび高エネルギーイオンの特性を勘案

捉電子がイオンに及ぼす効果は未解明であった．そこで，しつつ，コア加熱を最適化するターゲット設計が必要であイオンの衝撃波面での反射による加速が，捕捉電子によっる．

てどのような影響を受けるかを調べた［３２］．まず，物理的考察を行い，イオンの反射が多次元電磁擾乱によって促進３．
８ 衝撃波と粒子加速のシミュレーション

されることを示した．次に，空間２次元・速度３次元の相

太陽フレアや超新星爆発など，宇宙では様々な場所で無対論的電磁粒子シミュレーションを用いて，斜め衝撃波中

衝突衝撃波による粒子加速が起こっていると考えられていの電子の運動と電磁場の発展，イオンの運動を詳細に調べる．しかし，無衝突衝撃波の形成過程や粒子加速の機構は，た．また，イオン反射における２次元電磁擾乱の効果を明

未だ完全な解明には至っていない．一方，相対論的電磁粒確にするため，１次元的に平均化した電磁場中のテストイ

子コードを用いたシミュレーションによって，磁気音波衝

オンの軌道を計算した（テストイオンは２次元電磁擾乱を（a）ターゲットの電子密度（b）
イオンビーム強度の時間発展（c）
時間平均したイオンのエネルギー分布

図１０ ポンドラモーティブ力によるイオン加速．
表１

統合シミュレーションによるコア加熱の評価．
シース電場
gap なしgap 有り電子のみ
電子＋イオン
イオン向上率％

温度上昇
keV

0.15
0.19
27

0.08
0.16
100

８２７

ポンドラモーティブ力
5×10191×10205×1020
0.05
0.06
20

0.11
0.15
36

0.35
0.48
37

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

ルの屈曲や乱流磁気リコネクション，乱流ダイナモの影響が広く調べられている．ここでは，Hall 効果が乱流に与える影響を調べた２編の論文［３３，３４］について報告する．
第一の論文では，一様等方 MHD 乱流におけるスケール
間エネルギー伝達について，直接数値シミュレーションに
よって調べた．この直接数値シミュレーションでは，非圧縮性 Hall MHD 方程式
(
-(
( -(
" '$
#! #
$ ( )!'"!+" '"' &()，
"$
('
)%
(
/
(
/
(
,
# * * ()
)
)
（１）
(
-*
#!,
(
/
*

（２）

(
-( (
" (，
(
/
/
) (
(

（３）

(
!(
(
"*
#!%
，
(
)
*
(
,
(
/
)

（４）

(
!*
#!，
(
/
*

（５）

"(#!%
-(!%
!)"&
$
#$
(
)
*$
(%
(.

（６）

&()#

を，擬スペクトル法と Runge-Kutta-Gill 法で数値的に解いた．ここで記号 !(および -(はそれぞれ，磁場ベクトルおよび速度場ベクトルの第(成分を，記号'，&，%# は粘性率，
電気抵抗率，Hall パラメータである．記号 $()および %()*
は Cronecker のデルタ記号と Levi-Civita の反対称テンソル
であり，添え字が繰り返されたものについては和をとる．
図１
１ 衝撃波中の電磁場のプロファイルとイオンの速度分布．２次元電磁擾乱を感じている 2Ds イオンの方が，感じていな
い 1Dt イオンに比べて，より多くの粒子が衝撃波面で反射され加速されている．

まず，格子点数 %$#&"#$で外力を課さない減衰性乱流のシミュレーションを，同一の初期条件，同一の粘性率，電気抵抗率，時間刻みで行うことにより，Hall MHD 乱流!
!&）および MHD 乱流（%（%）の２つの場合を初期
# #!
# #!
値問題として厳密な比較対象可能な形で計算し，エネル

感じない）．図１１は，２次元電磁粒子シミュレーションのギースペクトル及びスケール間のエネルギー伝達を調べた

イオン（2Ds イオン）とテストイオン（1Dt イオン）の速度［３３］．この論文では，一様等方性Hall MHD乱流の磁気スペ.+*#"のイオンが反射されたイオンであ分布を表す．.,"

クトルが波数 *について *!("$となる傾向を確認した．まるが，2Ds イオンの反射粒子が 1Dt イオンのそれよりも多た，ローパスフィルターを用いて高波数成分を除去したエいことを示している．すなわち，反射イオンの数が，２次

ネルギー伝達関数（これは各 *でのをエネルギーの増減を元電磁擾乱によって増加することが確認できた．また，捕表す）を構築したところ，カットオフ波数が高波数である

捉電子が増えると主パルス中の２次元電磁擾乱は大振幅と場合には磁気エネルギー伝達関数の低波数成分側にはあまなり，2Ds イオンと 1Dt イオンの差が増加することも示し

り影響がない一方，カットオフ波数が低くなり過ぎると最た．

低波数（最も大きなスケール）のエネルギーの増減まで著

３．
９ 一様等方 MHD 乱流のシミュレーション（a）
運動エネルギー伝達関数および（b）磁気エネルギー伝

しく影響が及ぶことを明らかにした．図１２は文献［３３］より，
MHD 乱流は，宇宙プラズマ，太陽風や液体金属の流れ

達関数のそれぞれについて，カットオフ波数を *)#'%おなどの様々な流れに見られる．LHD のバルーニングモード

よび *)#$#とした場合の結果である．この図及び他の解シミュレーションにおいても，不安定モードの非線形飽和析から，カットオフ波数にある閾値があり，この閾値より段階においては，乱流と言ってよい状態が見られる．この低いカットオフ波数のローパスフィルターを作用させるこ中で特に一様等方乱流は，流体力学の乱流との対比も念頭とは低波数側のエネルギー増減を（したがってプラズマのに調べられてきた，古典的な乱流の一つである．近年では，

運動そのものも）著しく変化させる可能性がある事がわ

微視的効果が乱流に与える影響が重要視されるようになっ

かった．このような閾値の存在は，場の方程式を粗視化しており，特にイオン表皮長（Hall）効果を取り入れた Hallて，微細構造を現象論的モデルで置き換えるラージ・エ

MHD モデルは，太陽風における磁気エネルギースペクト

ディ・シミュレーションなどを考える際に重要な指標とな
８２８

Project Review

3. Sophistication of Computational Science and Fundamental Physics Simulations

S. Ishiguro et al.

図１
２ ローパスフィルターを作用させた磁場および速度場から構築した（a）運動エネルギー伝達関数および（b）磁気エネルギー伝達関数．

る．
次に，大規模な数値シミュレーション（"$!"!
#%$）の直接数値シミュレーションにより，Hall 効果がもたらす空間構造の変化について研究を行った［３４］．Hall パラメータ以外のパラメータを上記と同様に固定し，初期値問題として相互対象比較が可能なシミュレーションをいくつかのパラメータについて実施した．この中には，一流体 MHD（#&!!）シミュレーションも含まれる．このシミュレーショ
ンにおいて，Hall パラメータが大きいほど磁気エネルギー
スペクトルの高波数部が励起されるものの，運動エネル
ギースペクトルについては明確な変化がなく，Hall パラ
メータが違っていても相互に区別がつきにくいものであった．他方，空間構造を渦度（速度場の回転）の強さ（エンストロフィー密度）及び電流密度（磁場の回転）の等値面で可視化した と こ ろ，電 流 密 度 よ り も 寧 ろ エ ン ストロ
フィー密度に明確な違いが発見された．図１３（a）および（b）
はそれぞれ，MHD 乱流および Hall MHD 乱流におけるエン
ストロフィー密度（色付きの等値面）および電流密度（灰色の等値面）を，格子点数 5123 の領域について描画したものである．図１３（b）を，全格子点数 10243 を使用し，拡大表示したものが図１３（c）である．MHD 乱流においては，磁場と速度場の凍り付きの条件から，エンストロフィー密度と電流密度の等値面はそれぞれほぼ一体化し，層状構造を形成して い る．他 方，Hall MHD 乱 流 の 場 合，エ ン ス トロ
フィー密度の等値面は電流密度の等値面との一体化を失

図１
３ エンストロフィー密度（色付きの等値面）および電流密度，（b）（灰 色 の 等 値 面）．（a）
MHD 乱 流（格 子 点 数 5123）
Hall MHD 乱流（格子点数 5123），（c）
Hall MHD 乱流（全格．
子点数 10243）

い，流体乱流に良く見られる管状渦構造となる．また，電流密度は電流層の名残はあるものの，高波数の励起にともなって，短冊状の小片となって等値面が描かれている．これは，Hall 項の導入によって凍り付き条件が変わり，エンストロフィー密度と電流密度の等値面が単純に重なり合う

皮長スケールよりも高波数の領域（#""!
#!）にのみ発生理由がない事に加えて，磁場の高波数成分が Hall 効果にするわけではなく，より低波数領域にまで変化が発生してよって励起されるため，磁場の散逸効果と，これに付随しいることが確かめられた．これは，高波数部における Hallて発生し得る磁気リコネクションによる磁場のトポロジー

効果が，低波数領域まで含めた大規模な構造遷移を引き起変化が顕著に現れる一方，渦渡場（速度場の主に高波数成こしている証左である．

分）は磁場との凍り付きから解放され，Kelvin-Helmholtz

更に，これらの構造変化が高波数成分の性質に及ぼす影

不安定性などのプロセスを経て管状渦構造が形成されるた響の一端を，高波数成分が作る場を Smagorinsky 型のサブめと考えられる．更に，ローパスフィルターを使った解析

グリッドスケール（SGS）
モデル［３５］に射影することで評価から，この構造変化は Hall パラメータが指定するイオン表

した［３４］．この解析から，Hall MHD 方程式
（１）
‐
（６）の右
８２９

Journal of Plasma and Fusion Research Vol.92, No.11 November 2016

［１
１］H. Miura et al., 24th International Conference on Numerical
Simulation of Plasmas (Aug. 12-14, 2015, Golden, Colorado,
U.S.A.) (Invited talk).
［１２］H. Miura et al., J. Comput. Phys. 316, 385 (2016).
［１３］A. Ito et al., Plasma Fusion Res. 3, 034 (2008).
［１４］A. Ito and N. Nakajima, Plasma Phys. Control. Fusion 51,
035007 (2009).
［１５］R. Goto et al., Phys. Plasmas 22, 032115 (2015).
［１６］K.V. Roberts and J.B. Taylor, Phys. Rev. Lett. 8, 197 (1962).
７］P. Zhu et al., Phys. Rev. Lett. 101, 085005 (2008).
［１［１８］R. Goto et al., Plasma Fus. Res. 9, 1403076 (2014).
［１９］T. Johzaki et al., Nucl. Fusion 51, 073022 (2011).
［２０］S.P. Hatchett et al., Phys. Plasmas 7, 2076 (2000).
［２１］R.A. Snavely et al., Phys. Rev. Lett. 85, 2945 (2000).
［２２］S.C. Wilks et al., Phys. Rev. Lett. 69, 1383 (1992).
［２３］H. Sakagami et al., Nucl. Fusion 49, 075026 (2009).
［２４］H. Sakagami et al., Proc. 41st EPS Conf. on Plasma Physics,
P1.102 (Berlin, Germany, June 23-27, 2014).
［２５］H. Sakagami et al., J. Phys.: Conf. Series 688, 012096 (2016).
［２６］T. Johzaki et al., Laser Part. Beams 25, 621 (2007).
［２７］Y. Ohsawa, Phys. Rep. 536, 147 (2014).
［２８］N. Bessho and Y. Ohsawa, Phys. Plasmas 6, 3076 (1999).
［２９］K. Shikii and M. Toida, Phys. Plasmas 17, 082316 (2010).
［３０］M. Toida and J. Joho, J. Phys. Soc. Jpn 81, 084502 (2012).
［３１］M. Toida and J. Inagaki, Plasma and Fusion Res. 9, 341026(2014).
［３２］M. Toida and J. Inagaki, Phys. Plasmas 22, 062305 (2015).
［３３］H. Miura and K. Araki, Plasma Phys. Contrl. Fus. 55,
014012 (2013).
［３４］H. Miura and K. Araki, Phys. Plasmas 21, 072313 (2014);
http://dx.doi.org/10.1088/0741-3335/55/1/014012.
［３５］F. Hanba and M. Tsuchiya, Phys. Plasmas 17, 012301(2010).
［３６］H. Miura et al., The 24th International Conference on Numerical Simulation of Plasmas (Sep. 12-14 2015, Golden, Colorado, U.S.A.) (Invited talk).
［３７］H. Miura et al., J. Comput. Phys. 316, 385 (2016).

辺に現れる非線形項の高波数成分は，Smagorinsky 型 SGSモデルで効率よく表現可能であるものの，Hall 項についてはカットオフ周波数が比較的高波数の場合にのみ有効であるという結果が得られた．
これらの一様等方性 Hall MHD 乱流の性質，特に Hall
項が引き起こす構造遷移などは，LHD などトーラスプラズ
マの短波長不安定性の飽和などにも深く関係していると考えられる．一様等方性はトーラスプラズマとの比較の上で条件が異なりすぎるため，強い背景磁場を課した一様非等
方 Hall MHD 乱流の研究を行い，上記と同様の結果を得た［３６，
３７］．また，拡張 MHD モデルによる Rayleigh-Taylor型の不安定性成長と飽和のシミュレーション研究への応用も進行中である．

参考文献［１］石黒静児，大谷寛明：プラズマ・核融合学会誌 90, 325(2014).
［２］H. Ohtani et al., Contrib. Plasma Phys. 56, 692 (2016).
［３］M. Shoji et al., J. Nucl. Matter 463, 861 (2015).
［４］S. Usami et al., J. Phys. Conference Series 561, 012021
(２
０
１
４).
［５］S. Ishiguro and H. Hasegawa, J. Plasma Phys. 72, 1233
(2006).
［６］H. Hasegawa and S. Ishiguro, Plasma Fusion Res. 7,2401060 (2012).
［７］H. Hasegawa and S. Ishiguro, Phys. Plasmas 22, 102113(2015).
［８］H. Miura et al., Theory of Fusion Plasmas Joint VarennaLausanne International Workshop (Varenna, Italy, August28-September 1, 2006); AIP Conference Proceedings 871,
157 (2006).
［９］H. Miura and N. Nakajima, Nucl. Fusion 50, 054006 (2010).
［１０］H. Miura et al., The 25th IAEA Fusion Energy Conference
TH/P5-17 (Oct. 13-18, 2014, St. Petersberg, Rusia).

８３０

J. Plasma Fusion Res. Vol.92, No.11 (2016)８
３１

プロジェクトレビュー 核融合科学研究所における数値実験炉研究プロジェクト

４．まとめと展望
堀内利得
自然科学研究機構核融合科学研究所（原稿受付：２
０
１
６年７月１
３日）

第２期中期計画が始まった平成２２（２０１０）年に，数値実験

LHD プラズマ実験をはじめとする様々な実験結果との比炉研究プロジェクトの将来計画の研究目標として，大きな較による信頼性の検証（validation）を通じて，より高精度

３つの柱を設け段階的に達成することを計画した（図１）．

で信頼性のあるシミュレーションコードへの開発と改良を第１段階では，MHD 平衡・安定性，コアプラズマ輸送，乱進め，さらに（３）開発検証されたコードや明らかにされた流，高エネルギー粒子，加熱，周辺プラズマ輸送，プラズ

研究成果を数値実験炉に組み込み，LHD プラズマのさらなマ壁相互作用等の核融合プラズマを構成する要素物理の解る最高性能化研究や現在進められている原型炉設計に精度明とそのためのシミュレーションコード整備である．このの高い予測や重要な知見を与えるためのモデルやインター

研究は，図１からわかるように，第２期中期計画の主目標

フェースの構築，あるいは必要とされる周辺基盤技術・研となっていた．第２段階としては，得られた物理要素の統

究環境の整備等を行うことを，その基本戦略としてきた．

合化に向けた研究と核燃焼プラズマへの応用である．この同時に，これらの研究目標を実現するための研究体制の研究は，第３期中期計画期間における主目標となる．統合構築・整備を進めてきた．この間，国内外の研究者の英知化に向けた研究は，２．
５章に記載したように，既に，第２期を結集して数値実験炉研究プロジェクトに関する共同研究

中期計画の中でも着実に進めてきた．

を推進することはもとより，シミュレーション研究者の複さらに，第３段階として，「高精度・高信頼化された数数のタスクグループへの参加，タスクグループ間の研究交

値実験炉の構築」に向けた研究をプロジェクト全体で強力流，LHD 実験グループや炉工学研究グループとの研究交流に推し進めることにより，ヘリカル型核融合原型炉（以下を積極的に取り入れ，研究活動の活性化を図ることによでは，「原型炉」と略記）の設計へ大きく貢献することをり，研究目標の達成へ向けたプロジェクト研究を着実に進研究の柱として掲げている．既に，核融合工学研究プロ展させた．さらに，この進展を支えるものとして，プラズ

ジェクトを中心として推進している原型炉設計の開発研究

マシミュレータの性能向上をはじめとする研究環境の整備

に対し，これまでに数値実験炉研究プロジェクトにおいてやプラズマシミュレータ等を共同研究で利用するための支開発されてきた幾つかのシミュレーションコードを供与す

援体制の整備・拡充を進めた．特に，共同利用施設であるるなどの貢献をしている（２．
５章）．第３期中期計画の終盤

プラズマシミュレータは，平成２４年度および平成２７年度以降の主研究計画となっている．

に，それまでのそれぞれ４倍及び８倍の性能を有するシス

この研究計画を具体化するに当たっては，核融合プラズ

テムへと性能向上し，常に最先端の利用環境を全国の共同マを構成する多くの個別要素物理現象の解明を基本としな研究者に提供してきた．

がら，（１）世界各地で開発されている類似コードとのベン

これらの研究活動の成果は，多くの学術研究論文や国際

チマークテストによる理論検証（verification）および（２）

会議での発表として報告されており，その内の代表的なものに関しては２〜３章にて紹介した．その結果，数値実験
炉研究プロジェクト関連の論文の総数は，平成２２（２０１０）年
から２７（２０１５）年までの第２期中期計画期間の６年間で４８０編に達している．更に，これらの研究成果により，主たる国際会議での６１件の招待講演の推薦や２３件の学会賞等の受賞を受けるなど，プロジェクト研究体制の下，国際的に高いレベルの研究をプロジェクト全体として進めてきた．第３期中期計画の６年間（平成２８（２０１６）年〜平成３３（２０２１）年）
では，引き続き，国内外の研究者の英知を結集し，LHD
実験グループや炉工学研究グループとの連携の強化を図る
ことにより，「核融合プラズマ閉じ込めの物理機構解明とその体系化」および，ヘリカル型核融合原型炉の設計へ貢献する「高精度・高信頼化された数値実験炉の構築」に向

図１

第３期中期計画のローマップにおける位置づけ．

けたプロジェクト研究を強力に推進する．

4. Summary and Prospect HORIUCHI Ritoku
National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail: horiuchi.ritoku@nifs.ac.jp

８３１!2016 The Japan Society of Plasma
Science and Nuclear Fusion Research