J. Plasma Fusion Res. Vol.95, No.2 (2019)５
７‐９８

プロジェクトレビュー

IFERC 計算機シミュレーションセンターの高性能計算機システム Helios を利用したシミュレーション研究プロジェクト
Simulation Research Projects Using the High-Performance
Computer System Helios at IFERC-CSC

編集：量子科学技術研究開発機構 六ヶ所核融合研究所 核融合炉システム研究開発部
BA 計画調整グループ，IFERC 事業チーム，IFERC-CSC ホームチーム
Edited by QST BA Project Planning and Coordination Group, IFERC Project Team, and IFERC-CSC Home Team
（原稿受付：２
０
１
８年１
０月１
２日）

幅広いアプローチ協定に基づいて青森県六ヶ所村の国際核融合エネルギー研究センター（IFERC）の計算機
シミュレーションセンター（CSC）に設置された高性能計算機システム Helios は，２０１２年１月に運用を開始
し，２０１６年１２月末に計画通り運用を終えた．日欧の磁場閉じ込め核融合シミュレーション研究に供された５年の間，Helios は高い利用率を維持し，炉心プラズマ物理から炉材料・炉工学にわたる広い分野で多くの研究成果に貢献した．本プロジェクトレビューでは，５年間の運用，研究成果への貢献について概観したのち，先のプロジェ
クトレビュー以後に国内の大学や研究機関において Helios を利用して進められたシミュレーション研究プロジェ
クトとその成果を，各研究プロジェクト毎にその目的，用いられる計算手法，これまでの研究成果，そして今後必要とされる計算を紹介する．

はじめにであり，核融合研究専用機としては世界最速であった．増

青森県六ヶ所村に位置す る 国 際 核 融 合 エ ネ ル ギー研

強システムについては２０１６年８月号掲載の講座［２］に利用究センター（International Fusion Energy Research Centre:

例も含めて紹介されているので参照されたい．図１は，IFERC）の計算機シミュレーションセンター（Computa-

Helios の通常運用が始まった２
０１２年４月から運用を終了しtional Simulation Centre: CSC）は日欧間のより広範な取組

た２０１６年１２月までの Helios 本体システムの月平均稼働率

を通じた活動（Broader Approach: BA）協定に基づいて実（Availability）と月平均利用率（Utilisation）を示してい施されている事業で，磁場閉じ込め核融合研究専用計算機

る．２０１２年７月までは５０％程度の利用率にとどまっていたとして最新鋭の計算機を提供し，核融合研究への大規模・高効率シミュレーションの活用を促進することを目的としている．この目的のために設置された IFERC-CSC 高性能
計算機 Helios（愛称：六ちゃん）は２
０１２年１月に運用を開始し，計画通り２
０１６年１２月末に５年間の運用を終えた．
Helios 設置の経緯，機器構成，日欧の責任分担などの詳細については２
０１５年１１月号掲載の解説記事［１］に譲り，ここでは Helios の概要を述べるに留める．Helios は，CPUのみ搭載する本体システムと，CPU の 他 に メ ニ ーコア（Many Integrated Core; MIC）アーキテクチャをベースにしたコプロセッサも搭載する増強システムで構成されていた．Helios 本体システムの２
０１２年の運用開始時におけるLinpack 性能は 1.237 PFLOPS で，これは当時の TOP500

図１ ２
０
１
２年４月から２０１
６年１２月までの Helios の月平均稼働率（Availability）と月平均利用率（Utilisation）．

ランキングで世界１２位，国内２位（１位は京コンピュータ）
５７!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

ものの，２０１２年８月からは８０％を超え，その後も緩やかな数も第１サイクルの１８８名から順調に増加し，第３サイク

上昇傾向となり，２０１５年４月以降はほぼ９０％以上の高い利ルでは４００名に達し，第５サイクルでは最多の４５４名となっ用率を維持した．

た．特に欧州のプロジェクト数，利用者数がともに第３サHelios を利用する国際選考枠の研究プロジェクトの選

イクルで大きく増加しているのは，欧州域内の大型計算機

定，計算資源の割当は，以下のサイクルごとに日欧各５名が運用を停止し利用者が移ってきたためと考えられる．

の 委 員 か ら 構 成 さ れ る 常 設 委 員 会（Standing Commit-

このように Helios が安定に運用され，多くの研究者が利tee）によって決定されていた（国際選考枠とは別に日本の用し，高い利用率を維持した結果，図４に示されるように会計年度で切り替わる日本枠もあった．日本枠の各サイク

Helios を利用して得られた研究成果を発表した査読付き論

ルの期間についての詳細は２０１６年３月号に掲載のプロジェ

文数も毎年，堅調に増加した．Helios を利用した後の論文

クトレビュー［３］を参照）．

の作成・掲載には時間を要するため，第１サイクルが終了
１２年１１月までに報告のあった査読付き論文は３１報でした２０

第１サイクル ２０１２年４月〜２０１２年１１月あったものの，その後は１年で約１５０報づつ増え，２
０１７年２

第２サイクル ２０１２年１１月〜２０１３年１１月

月には累計で６３９報に達している．日欧の比較では，利用者

第３サイクル ２０１３年１１月〜２０１４年１１月

数の多い欧州側の論文数が日本側論文数の２倍以上となっ第４サイクル ２０１４年１１月〜２０１５年１１月ている．なお，日欧の論文数の比（２０１７年２月までの累計

第５サイクル ２０１５年１１月〜２０１６年１２月

で４４８/１９１）は，プ ロ ジ ェ ク ト 数 の 比（第５サ イ ク ル で

図２は国際選考枠のプロジェクト数の推移を示してい

８３/３５）や利用者数の比（第５サイクルで３０３/１２６）に近い．

る．赤色で示される日本側プロジェクトの数は日本の研究

これはプロジェクトあたりの論文数では，日欧で同程度で機関に所属する研究者が代表を務めるもの，青色で示されあることを示唆するが，この点は後でさらに分析する．査ているのは欧州の研究機関に所属する研究者が代表を務め読付き論文のうち日欧共著の論文数の推移を示したものがるプロジェクトの数である．第１サイクルでは日欧合わせ

図５である．およそ１０％程度が日欧共著論文になっていて５８だったのが，第２サイクルでは７３となり，第３サイクる．興味深いことに日欧共著論文は日本側プロジェクトにル以降は１００を超え第１サイクルの２倍程度まで増加した．

よるものが多く，欧州側プロジェクトによるものの２倍程

日欧の比較では，研究者数が多い欧州が第１サイクルの時度となっている．次に掲載誌別の論文数を図６に示す．上からプロジェクト数で日本を上回っていたが，第３サイク

位には，Physics of Plasmas，Nuclear Fusion，Plasma Phys-

ル以降その差は顕著になり，欧州のプロジェクト数が日本

ics and Controlled Fusion などプラズマ・核融合分野で著名のプロジェクト数の倍以上となった．図３は利用者数の推移を示している．プロジェクト数の増加に対応して利用者

図２

国際選考枠のプロジェクト数の推移．各サイクルで採択されたプロジェクト数を表す．

図３

Helios 利用者数の推移．各サイクルにおける登録利用者数
を表す．

５８

図４

査読付き論文数の推移．各集計時点までの累計数を表す．

図５

査読付き論文数（日欧共著）の推移．各集計時点までの累計数を表す．

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

図６

掲載誌別の論文数．

図７

BA Project Planning and Coordination Group et al.

研究分野ごとの論文数．

な論文誌が並んでいる．上位６位までの論文誌の順位は解説記事［１］にある２０１５年９月までのものと変わっていない．
Physics of Plasmas と Nuclear Fusion は日欧共に 多 いが，Plasma Physics and Controlled Fusion と Fusion Engineering and Design は欧州が多く，Plasma and Fusion Researchは日本が多いなど地域性が現れている論文誌があることがわかる．
さて，プロジェクトの研究分野は以下のように分類されている．
!プラズマ乱流および関連する輸送過程

図８!高速粒子に関する物理!線形，非線形および拡張電磁流体力学

研究分野ごとの１プロジェクトあたり論文数．

ジェクトの成果については，前回のプロジェクトレビュー!周辺プラズマの物理［３］において第３サイクルまでの国内の研究プロジェクト!加熱と電流駆動で得られた成果を中心にまとめられた．今回の報告では，!核融合プラズマの統合モデリング

その後の第４サイクル，第５サイクルに実施された国内の!核融合炉材料

研究プロジェクトで得られた研究成果を中心にまとめられ!核融合炉工学

る．表１は，本プロジェクトレビューで報告される国内シ

これらの分野ごとの論文数（２０１７年２月までの集計）を示

ミュレーション研究プロジェクトの一覧である．前回のプしたものが図７である．統合モデリング分野の論文数で日

ロジェクトレビューと同様の形式で，各研究カテゴリー毎

本側が上回る他は，全体の論文数が多い欧州が日本を上に，研究プロジェクトを研究代表者（Principal Investiga-

回っている．特に加熱・電流駆動，核融合炉工学の分野で

tor: PI）の五十音順に示している．同一の研究プロジェク

欧州の論文が多いことがわかる．次に分野ごとの論文数をト名でサイクルを継続している場合もあれば，サイクル毎

１プロジェクトあたりで示したものが図８になる．なお，に異なる研究プロジェクトを提案している場合もある．次

ここでは複数のサイクルで継続されたプロジェクトは１つ章以降に各プロジェクトの内容が表１の順に，研究プロのプロジェクトとして扱っている．これを見ると，乱流，ジェクトの目的，用いられている計算手法，これまでに得高速粒子，MHD，周辺プラズマの分野では，論文の絶対数られた成果，そして今後必要とされる計算の記述に絞ってでは欧州が上回っていたものの，１プロジェクトあたりの記述される．なお，表１に記載の研究プロジェクト名は日

論文数では日欧とも同程度（３〜４）であることがわかる．

本語であるが，次章以降の各プロジェクトの報告にある英MHD に関しては日本は欧州を少し上回っている．核融合

語のプロジェクト名が本来のもので，括弧書きされている

炉材料については，日欧共に他の分野より明らかに多く，ものは各プロジェクトの短縮名になっている．本プロジェ

特に日本の論文数は７に達している．一方，統合モデリン

クトレビュー及び前回のプロジェクトレビューにおいて，

グについては日欧共に前述の４分野より少ないなど分野に

欧州側の研究プロジェクトの内容については触れていない

よる違いが表れている．平成３０年版の科学技術白書［４］にが，欧州のプロジェクトも含めたオーバービューが第２６回

おいて我が国の論文数が減少傾向にあり，国際ランキング

IAEA 核融合エネルギー会議（京都，２
０１６年）
［５］
において

が低下しているとの指摘があるが，ここでの結果を見る

報告されている．また IFERC-CSC 事業の全体及び科学的と，核融合分野において日本の研究者の研究能力は欧州の成果についてはIFERC-CSC最終報告書［６］にまとめられて研究者に比べて劣ってはいないと考えられる．

いる．

以上が Helios の概要と Helios を利用したシミュレーショ

現在，IFERC-CSCではHeliosの後継となる新しい計算機

ン研究プロジェクトの全体概要となる．個々の研究プロ

シ ス テ ム Japan Fusion Reactor Simulator 1（JFRS‐１，愛
５９Journal of Plasma and Fusion Research Vol.95, No.2 February 2019称：六ちゃん‐II）が稼働している．IFERC-CSC 事業が当び研究成果は JFRS‐１に引き継がれ，さらに発展していく

初 の 計 画 通 り に 進 展 し，大 き な 成 果 を 上 げ た こ と が，ことが期待される．

JFRS‐１の予算が認められた大きな要因であることは間違（BA 計画調整グループ

いない．今後，Helios で開発された数値計算プログラム及

６０

宮戸直亮）

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

表１
番号
PI 氏名
１．乱流と輸送
１．
１ 石澤明宏１．
２ 井戸村泰宏
１．
３ 井戸村泰宏
１．
４ 今寺賢志１．
５ 糟谷直宏１．
６ 佐竹真介１．
７ 龍野智哉PI 現所属シミュレーション研究プロジェクト一覧．
研究プロジェクト名

実施サイクル頁

高ベータプラズマにおける乱流輸送
４，５Full-f ジャイロ運動論コードのベンチマーク
４，
５
Full-fジャイロ運動論シミュレーションによる乱流輸送の同位体効果に関する研究 ４，５
ジャイロ運動論コードによるトカマクプラズマの数値実験
４，５
乱流計測シミュレータプロジェクト
４，
５
ITER-BA のための先進的新古典輸送シミュレーション研究
４，
５

62
63
63
64
65
66
67
67
68
69
70

１．
１
２ 前山伸也
１．
１
３ 松岡清吉１．
１
４ 宮戸直亮１．
１
５ 矢木雅敏１．
１
６ 矢木雅敏１．
１
７ 渡邉智彦２．高速粒子２．
１ 藤堂 泰電気通信大学 スラブ配位における平行，垂直位相混合の数値的な調査
１
‐
５
核融合研
３次元プラズマにおける乱流輸送モデリングの研究
４，５
兵庫県立大学 磁気リコネクションの多階層シミュレーション
４，
５
核融合研
周辺プラズマの運動論的シミュレーションコードの開発
４
量研機構
JT-60SA，ITER，DEMO プラズマのシミュレーションのための統合モデリング ４，
５
コード開発
名古屋大学 高ベータジャイロ運動論シミュレーション
４，
５
原子力機構 ジャイロ運動論シミュレーションによる新古典トロイダル粘性研究
４，
５
量研機構
プラズマの炉心−エッジ結合の研究
４，
５
量研機構
磁気核融合プラズマにおけるマルチスケール乱流シミュレーション
１
‐
５
量研機構
トカマクプラズマにおける乱流拡散と帯状流抑制
５名古屋大学 複数のフラックスチューブを用いたジャイロ運動論的シミュレーション
４，
５

71
72
72
73
74
74

核融合研

75

２．
２
２．
３

核融合研
量研機構

１．
８
１．
９
１．
１
０
１．
１
１

沼波政倫
沼田龍介
長谷川裕記
林 伸彦

奴賀秀男
ビアワーゲ
アンドレアス

２．
４ 松山顕之３．MHD 現象３．
１ 相羽信行

京都大学
原子力機構
原子力機構
京都大学
九州大学
核融合研

BA Project Planning and Coordination Group et al.

ITER１
５MA 標準シナリオにおけるアルヴェン固有モードと LHD における高エネ ４，５
ルギー粒子駆動測地的音響モード
トカマクにおけるディスラプション時の逃走電子生成シミュレーション
１
‐
５

76

トカマクプラズマにおける高エネルギー粒子輸送の大局的非線形シミュレーション ４，
５

77

量研機構

垂直移動現象と逃走電子ビームのシミュレーション研究

４，
５

77

量研機構

トカマクプラズマにおけるプラズマ回転・３次元性・運動論効果が MHD 安定性 ２‐５に与える影響に関するシミュレーション研究
三次元 MHD 解析の拡張
４，５
三次元平衡の線型及び非線型 MHD 解析
４，５
高ベータトロイダルプラズマの非線形 MHD 現象
１‐
４
高ベータトカマクプラズマにおけるフィッシュボーンモードのようなシミュレー ４，
５
ションのための拡張ハイブリッドモデル開発
磁化プラズマの圧力駆動型モードの拡張 MHD シミュレーション
４，
５
統合コードを用いた ITER ディスラプションのシミュレーション研究
１‐
５

78

81
82

４，
５
１
‐
５
４，
５

82
83
84

３．
２
３．
３
３．
４
３．
５

市口勝治
核融合研
市口勝治
核融合研
佐藤雅彦
核融合研
ビアワーゲ
量研機構
アンドレアス３．
６ 三浦英昭核融合研
３．
７ 矢木雅敏量研機構
４．周辺プラズマ物理
４．
１ 菅野龍太郎
核融合研
共鳴摂動磁場を加えたトロイダルプラズマの粒子・エネルギー輸送解析
４．
２ 星野一生量研機構
原型炉設定に向けたシミュレーション研究４．
３ 矢木雅敏量研機構
トカマクプラズマにおける L-H 遷移のシミュレーション研究
５．統合モデリング
５．
１ 畑山明聖慶応義塾大学 SOLPS-IMPGYRO 統合コードによる不純物輸送シミュレーション
５．
２ 矢木雅敏６．核融合炉材料
６．
１ 伊藤篤史６．
２ 大澤一人６．
３ 鬼塚貴志６．
４ 鈴土知明６．
５ !山有道６．
６ 谷川 尚６．
７ 土屋 文
６．
８ 星野 毅６．
９ 渡辺淑之７．核融合炉工学
７．
１ 近藤恵太郎
７．
２ 佐竹信一７．
３ 宮本賢治

78
79
79
80

量研機構

統合ディスラプションコードのためのモジュール開発

５
２
‐
５

85
86

核融合研
九州大学
福井大学
原子力機構
核融合研
量研機構
名城大学
量研機構
量研機構

プラズマ‐物質相互作用の分子動力学研究
第一原理計算によるタングステン中の空孔拡散の研究
第一原理計算によるウラン合金中水素の挙動解析
プラズマ対抗材料および構造材料の照射・熱時効下での微細構造発達
第一原理計算によるプラズマ対向材料の研究
マルチスケールモデリングによるトリチウム増殖材料の性能評価
ナノ微細粒子分散型フェライト鋼中のヘリウムおよび水素の相対的安定性
リチウムイオン伝導体を用いたリチウム同位体分離に関する研究開発
材料照射相関の解明に向けた計算機シミュレーション研究

１
‐
５
４，
５
２
‐
５
４，
５
２
‐
５
２
‐
５
４
４，
５
４，
５

87
88
89
89
90
91
92
92
93

４，
５
４，
５
５

94
95
96

量研機構
並列モンテカルロ計算による IFMIF 原型加速器（LIPAc）の中性子工学解析
東京理科大学 剪断乱流の大規模スケール MHD シミュレーション
鳴門教育大学 核融合用負イオン NBIシステムの開発に向けた負イオン源プラズマの数値シミュレータ

６１

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

１．乱流と輸送１．
１ 高ベータプラズマにおける乱流輸送
Turbulent transport in high-beta plasmas（THETA）
プロジェクト代表：石澤明宏［京都大学］
国際枠

第４−第５サイクル

核融合プラズマにおいて，ベータと呼ばれる規格化された圧力は，核融合反応率やトカマクの定常運転に必要なブートストラップ電流生成などに直接かかわる重要なパラメータである．そして，微視的乱流は，磁場閉じ込めプラ
ズマのエネルギーおよび粒子の異常輸送を引き起こすことが理解されており，乱流輸送のベータ値依存性は核融合プ
ラズマ研究における中心課題のひとつである．この依存性を理解するためには，有限ベータになって初めて現れる磁場揺動がどのように乱流と帯状流のバランスに影響するか理解する必要がある．さらに，ストキャスティックになった磁場構造における磁力線方向の輸送など，磁場揺動によって直接生じる乱流輸送の理解も必要である．本プロ
ジェクトの目的は，この電磁的な乱流輸送のベータ依存性の理解を進めることである．本プロジェクトで用いる手法
は，磁場閉じ込めプラズマにおける乱流輸送の解析および予測に有用な電磁的局所ジャイロ運動論シミュレーションである［７］．
はじめに，これまでのトーラスプラズマにおける電磁乱流研究をレビューした論文［８］の概要を以下に述べる．有限ベータプラズマでは，磁場揺動が不安定性の線形成長率
や帯状流生成率など乱流輸送を決定する量に直接影響する．線形解析により以下のことが理解されている．イオン

図９ （a）
Cyclone base DIII-D case における乱流によるイオンエ
ネルギー拡散係数と電子エネルギー拡散係数のイオンベー
タ依存性．（b）
静電ポテンシャル揺動によるエネルギー拡
散係数と磁場揺動によるエネルギー拡散係数のイオンベー
タ依存性（参照論文［８］より転載）．

温度勾配不安定性（ITG）は本質的に静電的であるが，有限ベータでは磁力線折り曲げ安定効果によりその成長率はベータ値とともに減少する．捕捉電子モード（TEM）と電子温度勾配不安定性は有限ベータ効果の影響をほとんど受けない．一方，高ベータでは運動論的バルーニングモード

いの渦を引きちぎることにより飽和する．

という運動論効果によって補正を受けた MHD 不安定性が有限ベータプラズマにおいて乱流輸送を計算するために現れる．そのほか有限ベータでは微視的テアリングモード

は，磁場揺動を生成する電流密度揺動 y の評価が必要となが現れる．これらの不安定性はバルーニングパリティとテるため運動論的電子の取り扱いが必須となる．この運動論

アリングパリティに分類され，パリティを保存しながら線的電子効果導入を利用し，非軸対称プラズマである大型ヘ形成長する．そして，プラズマの非線形性によりパリティ

リカル装置およびヘリオトロン J プラズマにおける乱流輸は混合され，有限ベータにおけるバルーニングパリティ

送への運動論的電子の影響を調べた．解析においては実験

モード（ITG，TEM など）の飽和状態では，このパリティ

磁場配位を用いて，実験結果との比較を行った．その結果，混合が乱流による磁気面破壊の物理機構となる．また，パ実験で観測されるイオンおよび電子のエネルギー輸送を説

リティ混合は帯状流生成に影響を及ぼし，有限ベータでは明するためには，運動論的電子効果によるドリフト波不安

乱流と帯状流の相互作用を複雑にする．そして，高ベータ定性の成長率増加が必須であることを明らかにしたでは帯状流が弱くなり，数値シミュレーションにおいて飽

［１０，
１１］
．非軸対称プラズマにおける乱流輸送に対して今

和状態が得られなくなる（図９（a））．これはランアウェイ後必要とされる計算は，従来の計算は最も不安定な領域をと呼ばれ，有限ベータプラズマにおける乱流の飽和機構に通る磁力線の周りにおける乱流輸送の解析に限られてきた

対するジャイロ運動論シミュレーションの未解決問題であので，乱流輸送の磁力線ラベル依存性を明らかにすることる．この問題は非軸対称プラズマにおいては解決しているである．

ことを注意する［９］．局所近似を用いた場合，通常，最も不安定になるモードはバルーニング角がゼロであるが，有限ベータ非軸対称プラズマにおいてはモード構造が有限バ
ルーニング角を持ち，その結果，傾いたモード同士がお互
６２Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

１．
２ Full-f ジャイロ運動論コードのベンチマーク

量 ! の依存性は再現できないことがわかった．この結果をBenchmark test of full-f gyrokinetic codes（BMTFF）

受けて，第４サイクルではfull- #
シミュレーションのためのプロジェクト代表：井戸村泰宏［日本原子力研究開発機構］

新しいハイブリッド電子モデル［１３］を開発し，TEM（捕捉国際枠

電子モード）乱流のシミュレーションを開始した．

第４−第５サイクル

まず，ITG-TEM 乱流の線形解析を実施し，電子衝突時間%!
"の比で特徴付けられる TEM#とイオン通過時間 $%$Full- #ジャイロ運動論モデルは，巨視的なプラズマ分布と微視的な乱流成分を第一原理によって同時に発展させるモデルであり，従来の ##モデル（分布を固定し乱流成分のみを計算する）では扱うことのできない自己無撞着なプラ
ズマ分布と乱流の相互作用を扱うことができる．これまで，##シミュレーションに関しては乱流輸送のコード間ベンチマークにより計算結果のロバースト性が確保されているが，full- #シミュレーションでは以下の困難さのため，系
統的なコード間ベンチマークは行われていない．まず，full- #シミュレーションにおいては，乱流輸送だけでなく，新古典輸送，プラズマ分布形成，及び，それらの相互作用を含むマルチスケール・マルチフィジックス問題のベンチ
マークとなる．これに加えて，大域的シミュレーションに
おいては，境界モデルや熱源モデルをどのように一致させるかも大きな問題となる．本プロジェクトでは，full- #シ
ミュレーションの計算結果のロバースト性向上を目的とし，GT5D（原子力機構）と GYSELA（CEA）による熱源
駆動乱流シミュレーションのコード間ベンチマークを行った．
図１０は，GT５D と GYSELA による熱源駆動乱流シミュ
レーションの比較を示す．図１０（b），（c）から，熱流束は
主に雪崩的輸送によって引き起こされること，雪崩的輸送の伝搬方向が径電場構造によって周辺領域では外向き，炉心領域では内向きに逆転することなど，コード間で輸送の時空間的な特徴が一致していることがわかる．また，図１０（a）の周波数スペクトルについても，自己組織化臨界現象"
への遷移に関しを特徴づける!
から乱流を特徴づける!
!
#!
#て良い一致を得た．これは，full- #シミュレーションにおける自己組織化臨界現象やそれによって発生する雪崩的輸送がロバーストに再現できることを示す．この成果を Physics of Plasmas 誌に投稿した［１２］．今後は，自発回転を形成する運動量輸送等，乱流輸送と新古典輸送の相互作用が重要な問題のベンチマークに取り組む予定である．
１．
３ Full-f ジャイロ運動論シミュレーションによる乱流輸
送の同位体効果に関する研究
Study of isotope effects on turbulent transport using full-f
gyrokinetic simulations（GT5DISO）
プロジェクト代表：井戸村泰宏［日本原子力研究開発機構］
国際枠

第４−第５サイクル

GT5DISO プロジェクトでは full- #ジャイロ運動論コード
GT5D を用いて閉じ込め性能の水素同位体効果を明らかにすることを目標としてシミュレーション研究に取り組ん
図１
０ （a）
イ オ ン 熱 流 束 "i の 周 波 数 ス ペ ク ト ル 及 び GYSELA
（b）
，GT5D
（c）
により得られたイオン熱流束 "i の時空間発展．

だ．第３サイクルまでの結果から，断熱的電子モデルを用いた ITG（イオン温度勾配駆動）乱流では同位体効果のう
ち規格化ラーマ半径 $!の依存性は再現できるが，規格化質６３Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

閉じ込め核融合プラズマを高性能化する上で重要な鍵を握っているが，第一原理に基づいた乱流輸送シミュレー
ションでは再現されていなかった．これは，大規模な計算機資源を必要とする大域的乱流輸送コードの開発・実行が容易でないこと，さらには径電場に代表される背景分布と乱流の相互作用が物理的に解明されていないことに起因している．
そこで本研究では，ITB が実験的に観測されている反転
磁気シアプラズマにおいて，熱源駆動型乱流輸送シミュ
レーションを実施し，内部輸送障壁の形成メカニズムに関する解析を行った．図１２は，各時刻におけるイオン温度分
布を示している．安全係数分布は黒線で示した反転磁気シ
ア配位を用いており，#!"&!!から#!$#!!の間，図に示す

図１
１ ITG-TEM 成長率 !vti/R の温度勾配 R /Lti，電子衝突度 "＊e依存性．"＊
e ＝ 0 では粒子コード GT3D との良い一致を示し，"＊
e ＝ 0.03 では水素と重水素の TEM 成長率に同位体効果が示される．

領域に運動量の入射を行っている．この結果から，運動量を入射してもその領域で温度分布の急峻化がほとんど起きない一方で，入射を止め十分時間が経過した後，")(*近傍
で温度分布が急峻化することが新たに明らかとなった．この衝突安定化効果が水素（H）と重水素（D）で大きく異なる

れは磁気シアが正の領域からの運動量ピンチの効果で，回ことから ! 依存性が存在することを明らかにした（図１１）

転分布の急峻化が引き起こされ，その結果，力学的平衡を［１４］．次に，運動論的電子モデルの非線形シミュレーショ

介して径電場シアが形成されたためと考えられる．つまンの検証を目的として，同じパラメータの ITG 乱流シミュ

り，反転磁気シア配位の場合，内部輸送障壁の位置は運動

レーションを断熱的電子と運動論的電子で比較し，熱輸送量分布に関係はなく，")(*の位置で決定されると考えられ特性や温度分布は大きく変化しないが，乱流飽和をもたらる．

す径電場の形成機構や運動量輸送の機構が異なることを明次に内部輸送障壁の応答特性を検証するため，運動量入らかにした［１５］．最後に，運動論的電子モデルの実験的検射に関する変調数値実験を行った．図１３は，#!%!!!まで証を目的として，電子変調加熱数値実験を実施し，実験的運動量入射を行った後，入射を止め，#!'!!!から再度入に観測されている電子加熱による密度分布の尖頭化，トロ射を行った際の ")(*近傍における（A）温度勾配と乱流熱流

イダル回転トルクの反転を再現することに成功した．この束の関係，（B）
温度勾配と径電場シアの関係，および（C）

数値実験では，ITG 乱流から TEM 乱流への遷移によって，温度勾配の空間‐時間発展を示している．局所的な温度勾粒子輸送の向きが逆転し，それがトロイダル角運動量保存

配や径電場シアの変化よりも早い時間スケールで乱流輸送

則を介してトロイダル回転トルクを反転させる機構を明らが変化した結果，ヒステリシス性が確認された．これは熱かにした［１６］．

アバランチや乱流拡散などの非局所輸送がコア部から ")(*

今後，JT-60U で観測されている同位体効果の再現に向に流入することで，内部輸送障壁が崩れることによるもの

けて，水素プラズマと重水素プラズマにおける ITG-TEM

であり，非局所性が現象を支配していることが明らかと乱流の系統的なパラメータスキャンに取り組む予定である．

なった．

１．
４ ジャイロ運動論コードによるトカマクプラズマの数
値実験
Gyrokinetic Numerical Experiment of Tokamak（GKNET）
プロジェクト代表：今寺賢志［京都大学大学院］
国際枠

第４−第５サイクル

本 プ ロ ジ ェ ク ト で は，大 域 的 ジ ャ イ ロ 運 動 論 コ ード
GKNET を用いて，（A）
反転磁気シアプラズマにおける内部輸 送 障 壁（ITB）形 成 に 関 す る シ ミ ュ レ ー シ ョ ン［１７］，（B）熱源駆動型乱流の特性に関する統計解析［１８］，（C）帯状流，および測地線音響モードに対する磁気面形状効果の解析［１９］，（D）
GKNET への運動論的電子効果の実装，および大域的捕捉電子モード不安定性に対する線形解析を行った．本稿では，（A）の成果について報告する．
乱流を抑制することで輸送係数を低下させ，断熱層／粒図１
２ イオン温度分布（t ＝ 1600，t ＝ 2400，t ＝ 3200，t ＝ 4000）と熱／運動量ソース，シンク，安全係数分布．

子遮蔽層として作用する「内部輸送障壁」の形成は，磁場６４

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

は実験で行われているプラズマに摂動を加えることで非局所応答の様相を探る研究に対応したものである．周期的に加えた摂動に対して条件付き平均法によりプラズマ応答を抽出している．ヘリカルプラズマでドリフト交換型モード
は共鳴面近傍の有限幅に存在するので，複数のモードが存在する場合，磁場構造に依存した分布を持つ．その不安定性が分布する中で圧力ソース変調シミュレーションを行
い，プラズマ応答が径方向に伝播する様相を調べた［２３］．
不安定モードが駆動する熱流束の連鎖的応答が径方向の圧力伝播を加速する．ソース変調直後の大変動とともにその後の飽和状態に近づく過程での圧力小振幅変動を解析した．非線形飽和状態は多くのモードからの寄与の総和で飽和レベルを維持している．一方で揺動スケールの時間変動についてエネルギー時系列の相関解析を行った．図１４に各半径位置での圧力平均成分と代表的モードからの非線形寄

図１
３ qmin 近傍における（A）
温度勾配と乱流熱流束の関係，（B）
温度勾配と径電場シアの関係，および（C）
温度勾配の空間‐
時間発展．

与の時間変動の相関係数を示す．圧力変動はそれぞれの場
所で異なる，共鳴モードからの非線形寄与成分と相関関係にあることがわかった．非線形性による自己維持機構が働

１．
５ 乱流計測シミュレータプロジェクト

く中で，磁場分布に依存する不安定な単一モードがプラズ

Turbulence Diagnostic Simulator Project（TDS）

マ応答を主導し，径方向に変動を伝播する役割を果たしてプロジェクト代表：糟谷直宏［九州大学応用力学研究所］

いる様相が示された．加えるソース摂動の分布を変え，応日本枠

答の違いを比較した．プラズマ中心部に局在させた摂動に

第４−第５サイクル

対して動的応答のポロイダル断面２次元分布を求め，過渡
ITER のような核燃焼プラズマにおいては，実験計測が過程におけるバルーニング性を持つ熱流束の渦巻き構造を制約されるので，得られた実験データから最大限の情報を示した．さらに摂動を加える径方向位置を変えたシミュ

引き出すことが不可欠である．そのためシミュレーション

レーションを行い，摂動の伝搬可能域の違いと分布硬直性

データを用いてプラズマ乱流の理解を深め，そしてその結を見た．

果を実験観測に適用することで新たな発見をめざす手法はヘリカル型装置における実配位３次元シミュレーション有効である．計算機技術の発展に伴い可能となったプラズ

による定量的解釈を目的として数値診断を行った．LHD

マ乱流の大域的様相を模擬した大規模シミュレーション，

装置参照平衡磁場配位を用いた MHD 不安定性シミュレー

および実験計測技術の発展に伴って可能となった時間的，ションでは非線形飽和過程における不安定モード構造の数

空間的に分解能の高い揺動計測の両者を統合して，プラズ

値診断を行った．圧力駆動型バルーニング不安定性の駆動

マ乱流の理解を加速することができる．そこで我々は乱流と圧力分布緩和に伴う非線形飽和において，熱拡散係数の計測シミュレータ［２０］という乱流コードと実験計測を模擬

大きさにより励起モード数とエネルギー移送過程に違いがする数値診断モジュールから構成されるコード群の開発を現れることが明らかになった．また，ヘリカルプラズマでし，計測信号の解釈，定量的評価及び物理機構の解明を

揺動は３次元的分布を持ち，観測との対応が問題となるの行っている．本プロジェクトではトロイダルプラズマにおで，実配位を用いたジャイロ運動論シミュレーションによける非線形シミュレーションを行い，揺動場の３次元時系

り得られた揺動場の時系列データを対象とし，観測視線方

列データの取得，およびそのデータの数値診断により乱流向に積分した密度揺動信号と局所揺動の関係を比較した．

輸送ダイナミクスの理解を進めた．

実験的に観測されるのは線積分値のみであるが，磁力線方本研 究 の た め に乱 流 コ ー ド を Helios へ 移 植 し 実行した．研究対象によって複数種類のコードを用いている．乱流コード DISA では３場の簡約化磁気流体モデルを用いて，ヘリカル状磁場配位中での流れ関数，ベクトルポテン
シャル，圧力の時間発展を計算し，ドリフト交換型モードを解析する［２１］．方位角，軸方向に周期境界条件を課してフーリエモード展開をするスペクトル法を用いている．半
径方向には離散化した非線形偏微分方程式を解いている．
時間発展は予測子修正子法を用いており，線形部分は陰的に解いている．MPI+OpenMP のハイブリッド並列実行が

図１
４ 平均圧力成分とそれに対する各モードからの非線形寄与の関係［２３］．時間変動の相関の径方向分布をプロットしている．

可能である．このコードを用いて，プラズマ中心部に加熱摂動を加える大域的シミュレーションを行った［２２］．これ
６５

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

向の空間分布を利用した局所値の抽出法（磁気シア法）を)$''は左辺の !(全体，あるいはその径方向成分 %(
#
(+ を無
)
+

適用し，空間分解能を評価した．
このように ITER プラズマにおける乱流輸送の物理的理'&視していたが，FORTEC-3D は )%('項を無視せ"%(*#($解を深めるためのデータ解析技法の基盤研究を着実に進展ず，有限軌道幅効果を含めた正確な案内中心軌道を考慮しさせることができた．今後，トカマク配位およびヘリカルて新古典輸送計算が行える点が特徴である．

配位における不安定性の駆動力の違いによる多次元的な乱（１）新古典トロイダル粘性と摂動磁場の空間構造について

流輸送ダイナミクスの同定や実配位シミュレーションにお

NTVは+
，$"%)
$
##!$
$",
$
#-!$
#'*
#' !のように(#(#.

ける輸送の定量化を図る．

.

.

(-#
#
表される．
ここで，
，
，
(-$
'
$
#-$%&$$#' $%&$-! '"$'#+0,は磁気面平均を表す．これまでの研究では，磁気面平１．
６ ITER-BA のための先進的新古典輸送シミュレーショ
ン研究

均量のみが評価されてきたが，本研究で初めて NTV の空

Advanced Neoclassical Transport Simulation for ITER-BA間分布構造も評価してみた．NTV トルクの空間密度は，（ANTSI）
プロジェクト代表：佐竹真介［核融合科学研究所・総合研

"
( )
$
# "$
#-*という量から評価できる．ただし， "
)
(
$/)
( '

究大学院大学］

は Boozer 座標系のヤコビアン，$/$+$#+'，$はトロイダ

国際枠

ルフラックス 'の磁気面に囲まれる体積を表す［２９］．

第４−第５サイクル

ITER で計画されているように，３か所 TBM を挿入した
ITER や JT60-SA トカマクにおける核融合研究において状況での MHD 平衡磁場を計算したところ，挿入したトロ重要となる新古典輸送現象として，軸対称性の破れが生イダル位置に局在化した摂動磁場が，図１５のように形成さむ，新古典トロイダル粘性（NTV）がトロイダル回転に与れることがわかった．この摂動磁場は，内 側 の 磁 気 面

える影響がある．本課題では第１〜３サイクルにおいて，!
%）まで深く影響していることがわかる．一方，最$%（+$!

ヘリカルプラズマ用に開発されてきた新古典輸送コード!
%）においては，)$"&のトロイダル磁'%
外殻近く（+$!

FORTEC‐
３D［２４］を拡張して NTV の数値計算法を開発し，場リップルが卓越している．この２つの 磁 気 面 上 で のトカマクにおける磁場のトロイダルリップル成分による

NTV トルク密度を評価すると，摂動磁場の空間構造に対NTV を評価し，その影響を含めたトロイダル回転分布の応して NTV トルクが発生していることがわかった．この予測計算と JT-60U の観測結果の比較を通じて
［２５，
２６］，定

)となる共鳴有理面を持たない，計算例のように，*$(#

量的予測性の検証を行ってきた．第４，５サイクルでは，)）はポロイダル及びトロイダ
非共鳴なフーリエ成分［（("ITER で 設 置 が 計 画 さ れ て い る テ ス ト ブ ラ ン ケ ッ ト モルのモード数］の磁場摂動を与えた場合，その摂動磁場はジュール（TBM）のように，真空容器内に構造体を入れた

磁気軸近くまで入りこみ，NTV を発生させる．一方，共鳴場 合 に 生 じ る 局 在 化 し た 磁 場 摂 動 や，共 鳴 磁 場 摂 動

有理面を持つ )$$の RMP を JT60-SA に与えたケースで（RMP）を掛けた場合の JT60-SA について，NTV の評価

は，摂動磁場の空間構造と NTV トルクの空間分布の関係と，そのトロイダル回転への影響について研究を進めた．

は非共鳴摂動の場合ほど明確ではなかった．詳細な報告はまた，トカマクに限らず一般的な３次元配位における新古［２９］にあるが，理想 MHD で評価した MHD 平衡は，プラズ

典輸送計算の評価において，従来から用いられてきた局所マの応答によって共鳴モードの摂動磁場分布が真空を仮定近似などの近似モデルと，より第一原理的な大域新古典輸
送計算法を比較し，ドリフト運動論を解く際の様々な近似が新古典輸送の定量的評価にどのように影響するかについての研究［２７，２８］も行った．NTV がトカマクのトロイダル
回転分布に与える影響の評価については，プロジェクト
JIMS の項において報告される．
FORTEC‐
３D は５次元位相空間 )
+
"
&
"
(
"
"
%*における案
内中心分布関数の局所マクスウェル分布 '+"-*からのず
*)
れ$
'
%'
!'
* に対するドリフト運動論方程式
) %
'
,)$$
'
"% "%&"%(
*
#
("$
# ")
)
, - ) (
)
'
,)"'
$
'
!!%(
#
(""!)
* （１）
)
- *

を 2-weight $'法と呼ばれる数値手法で解くモンテカルロ
コードである．ここで !)$'*は線形化衝突項を表す．また，

図１
５ ITER に TBM を設置した条件で評価した，非軸対称摂動磁場成分（上段）と NTV トルク密度（下段）の，２つの磁気面上における空間分布．（(，&）は Boozer 座標系におけるトロイダルおよびポロイダル角．

'&，%(
は磁気ドリフト項の磁気面垂直成分（１）式左辺の %(および，接線成分を意味する．従来の局所近似計算手法で
６６

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

した場合と異なっており，これが NTV の評価に重要であることが計算結果から示唆されている．
（２）局所近似モデルの新古典輸送計算への影響（１）式のドリフト運動論方程式を解く際に用いられる近似の 影 響 を 統 一 的 に 調 べ る た め，
（１）を そ の ま ま 解 く%
のみを落とFORTEC-3Dを基に，左辺の径方向ドリフト!$$を 残 す Zero-Orbit-Width（ZOW）モ デ
し，接 線 成 分 !$ル，磁 場 ド リ フ ト 項 全 体 を 落 と す Zero-Magnetic-Drift（ZMD）モデル，そこからさらに，DKES コード［３０］のように，%!!ドリフト速度 !% "%!$!
$#を，%!$!
$#'で近&似した DKES-limit モデル（%!$ドリフトの非圧縮近似に相 当）の３通 り を１つ の コ ー ド で 解 け る 局 所 近 似 版
FORTEC を構築した．図１６に，LHD 配位における新古典
図１
６ LHD における，新古典粒子フラックスの径電場依存性のモデル間比較．GSRAKE はバウンス平均コード（ZMD と同様のモデル）の計算結果であり，Global が本来の大域的なFORTEC-3D 計算の結果である．

輸送の計算例を示す．ZMD，DKES モデルでは，径電場!
"領域と呼ばれる，ヘリカル捕捉粒"&( !において所謂 "子の寄与による極端に大きな新古典フラックスが生じているが，磁場ドリフトの接線成分を残した ZOW モデルではそのピークがシフトし，かつ全体的に大幅に小さくなり，
磁場ドリフト項を全て残した global 計算の傾向を再現していることがわかる．この傾向は電子よりイオンにおいて特に顕著であった．この研究結果は，ヘリカルプラズマの新古典輸送において重要な，電子とイオンの新古典フラック
スが釣り合う両極性径電場を評価する際に，径電場がゼロ$に近く，!% #!$に な っ て い る 場 合 に は 従 来 の ZMD やDKES モデルの評価は正しくないということを，明確に示した初めてのシミュレーション例である．
１．
７ スラブ配位における平行，垂直位相混合の数値的な調査
図１
７ ２次元自由減衰乱流における無衝突保存量の減衰則［３２］．
低衝突な場合（D），衝突と非線形ダイナミクスがバランスする場合（B）
ともに，無衝突保存量 W，E はいずれも時間のべきで表される時間発展を示す．

Numerical Investigation of Parallel and Perpendicular Phase
Mixing in Slab Geometry（NIPPPMSG）
プロジェクト代表：龍野智哉［電気通信大学］
国際枠

第１−第５サイクル
また，磁場に平行な位相混合とのカップリングについて

ITER のような先進核融合装置では，低衝突のため運動調べる大規模シミュレーションを行うため，コードの並列論的乱流における速度空間の構造を理解することが非常に化 効 率 を 高 め る 高 速 化 を 行 っ た．rmhdper を 用 い て重要になる．本プロジェクトでは，簡単化したスラブ配位

OpenMP による 自 動 チューニングシステムを 開 発 し［３４］，における，磁力線に平行，垂直な方向のそれぞれについて，

コードの部分ごとに並列度を変化させることで，全体を通位相混合の基礎物理を理解するための数値的な調査を行う

して同じ並列度を用いたシミュレーションよりも高速化でことを目指した．数値コードは，トーラス配位において世きることを示した（図１８）．同様のスキームを AstroGK に界的に広く使われているジャイロ運動論コード GS2 を簡単も移植中であり，今後はさらに大型のシミュレーションを化した AstroGK と，数値スキームのテストのために２次元行いたい．

の簡約化電磁流体コード rmhdper を用いた．
まず磁力線方向に一様性を仮定した２次元の静電揺動に

１．
８ ３次元プラズマにおける乱流輸送モデリングの研究

ついて，自由減衰乱流の理論・数値シミュレーション研究

Gyrokinetic Turbulent Transport Modeling in three-を行った［３１‐３３］．２次元配位では，無衝突保存量が２つ存

dimensional plasmas（GKMDLST）

在するため，磁場に垂直な位相混合のダイナミクスは双カプロジェクト代表：沼波政倫［核融合科学研究所］

スケードを示すと考えられる．この場合非線形相互作用の国際枠

第４−第５サイクル

保存性とカスケードの性質から簡単な現象論を構築することに成功し，無衝突保存量の減衰則を理論的に導くことが本プロジェクトでは，３次元配位における磁場閉じ込め

できた．数値シミュレーションは理論的な減衰則とよく一プラズマの異常輸送現象に対して，ジャイロ運動論的シ

致し，理論が正しいことが示された（図１７）．

ミュレーションを用いた数値計算を実施し，乱流輸送フ６７

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

スの温度勾配依存性を様々な径方向位置で調べた［３５］．同時に，LHD 実験における温度分布の計測誤差デ ータから，許容される温度勾配長の誤差幅を見積もることで，温度勾配の誤差内での輸送フラックスの変化量を特定した．
温度分布は，計測データから，径方向のイオン分布を規格化小半径の６次から１０次までの関数形と し て 仮 定 し て
フィッティングし，その誤差分布から温度勾配長の誤差を評価した．図１９で示すように，関数形に応じて誤差幅が大きく異なることがわかる．また，今回実施した広範な温度勾配長領域での数値解析から，実験で観測される輸送フ

図１
８ 自動チューニングシステムを用いたrmhdperの高速化効率［３
４］．B は高速化前，数字は固定されたスレッド数，A は自動チューニング時の時間ステップ数あたりの計算時間を示す．

ラックスを実現し得る温度勾配長を見積もった．図１９から，シミュレーションから予測される温度勾配長が，ほぼ全領域で，上で評価した誤差範囲内に留まっていることがわかる．近年のトカマクでの研究で，実験での温度分布のラックスのプラズマプロファイル依存性の詳細な解析と，計測誤差に伴う温度勾配の誤差幅の範囲内で，ジャイロ運

輸送予測やプロファイル予測を実施すると同時に，実験結動論によるシミュレーション結果と実験との輸送フラック

果との比較から，シミュレーション結果の検証および妥当スの差異に関する原因を説明し得るとの議論［３６］もある性の確認と，より高精度の乱流輸送モデルの確立を目的とが，LHD における本解析では温度勾配の誤差幅自体が大きした．特に本課題では，大型ヘリカル装置（LHD）におけ

くないことがわかり，温度勾配の誤差幅に対して，より正る乱流解析を，５次元位相空間上のプラズマ分布関数を解

確な評価手法の確立が必要であることが示唆された．

く局所乱流輸送コードである GKV を用いて行い，磁場配位として，VMEC コードにより計算された LHD の実際の

１．
９ 磁気リコネクションの多階層シミュレーション

放電配位を正確に模擬された磁場配位を適用した．

Multiscale Simulation of Mangetic Reconnection（MSMR）

イオン温度勾配不安定生（ITG モード）が駆動する乱流

プロジェクト代表：沼田龍介［兵庫県立大学大学院］

は，イオン温度勾配分布に極めて敏感である．ここでは，

日本枠

第４−第５サイクル

温度勾配長に関する広範なパラメータ領域においてジャイ
ロ運動論的乱流シミュレーションを実施し，輸送フラック

磁気リコネクション現象は，プラズマのエネルギー変換
機構や輸送現象に本質的な役割を果たしており，その理解は，核融合プラズマの閉じ込め性能の向上に欠かすことはできない．磁気リコネクションは本質的に多階層現象である．粒子間衝突や電子慣性といった微視的（粒子的な）効
果によって磁力線のつなぎ変えが起こり，大域的な磁場構造の変化を引き起こす．本研究プロジェクトでは，運動論モデルを用いて，マクロからミクロまでを包括的に取り扱う多階層シミュレーションを行った．磁気リコネクションの爆発的な時間スケールや，磁気リコネクションの結果として生じる磁気島の成長，飽和過程，および加熱や輸送といった磁気リコネクションを通したエネルギー分配過程などにおける運動論効果を明らかにすることが目的である．
本プロジェクトでは，電磁的ジャイロ運動論コード
AstroGK［３７］を用いて，２次元スラブにおける磁気リコネ
クションの計算を行った．AstroGK は連続体近似に基づいており，背景磁場に直交する面内ではフーリエスペクトル法を用い，時間積分は，線形項に対しては陰解法を，非線形項に対しては３次の Adams-Bashforth 法を用いている．
第１〜３サイクルにおいて，磁気リコネクションの非線形シミュレーションを行い，位相混合によるプラズマ加熱機構を実証した［３８］．電子は沿磁力線方向の位相混合により，分布関数に !!方向の急峻な構造を形成し，これが衝突

図１
９ GKV シミュレーションの結果に基づき得られたイオン温度
勾配分布の予測（上図）と実験観測から見積もられる温度勾配の誤差分布（下図）．実験観測から見積もられた温度
勾配（実線）とその誤差幅に対して，シミュレーションからの予測値（丸印）は，よい一致を確認できる．

によって熱化され加熱が起こる．高ベータプラズマにおいては，イオンも同様に位相混合によって加熱されるが，イオンの場合は有限ラーマ半径効果に起因する，磁力線に直６８Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

１．
１０ 周辺プラズマの運動論的シミュレーションコードの交方向の非線形な位相混合が重要になる．さらに，電流開発

シート構造の２次的な不安定性によるプラズモイド形成によって，加熱が強く引き起こされることが新たに明らかにDevelopment of Kinetic Dynamics Simulation Codes for Pe-

なった．このことは，プラズモイドが活発に生成される乱

ripheral Plasma（KINDSP）

流的な状況において，より効率的なエネルギー変換が起こ

プロジェクト代表：長谷川裕記［核融合科学研究所］

ることを示唆している．

日本枠

第４サイクル

第４〜５サイクルにおいては，磁気リコネクション過程に乱流が及ぼす影響を明らかにするため，ジャイロ運動論

境界層のプラズマ輸送や非接触ダイバータプラズマ過程

モデルに外力による乱流駆動項を実装する方法を検討しについては，多くの理論的，数値的な研究がなされてきた

た．ジャイロ運動論モデルと簡約化 MHD モデルの結果を，が，それらの多くでは，実験装置サイズの現象を扱う必要

乱流強度をパラメタとして比較できるようにするため，入性から，流体モデルが用いられてきた．しかし，そのよう力エネルギーが一定となるランダム外力項［３９］を採用し，な巨視的手法では，これら諸現象における微視的，すなわ

簡約化 MHD コードに実装した．動作確認を行った結果，

ち，運動論的な物理過程（例えば，粒子ジャイロ運動の効一定の外力によるエネルギー入力と粘性散逸が釣り合う定果，ダイバータ板近傍のシース電位形成，速度分布が非マ常状態において，理論的に予測されるエネルギースペクト

クスウェル分布となる現象，放射，電離，荷電交換，再結ル（磁場がない場合 !""!"，磁場がある場合 !""!"!!）

合などの原子分子過程等）を，自己無撞着に取り入れるこが得られることを確認した．引き続き同手法の AstroGKとは不可能であり，微視的過程の影響を評価することがでへの実装を進めている．

きない．本プロジェクトでは，周辺プラズマ現象におけるまた，関連する MIC プロジェクトとして，AstroGK コー

運動論的な素過程を明らかにし，それが周辺プラズマ現象

ドを MIC ノードに移植し計算効率の測定を行った．フラッに与える影響を評価するための運動論的プラズマシミュ

ト MPI 並列化されていたコードに対し，OpenMP+MPI

レーションコードの開発を行う．物理現象としては，非接によるハイブリット並列化を実施し，計算効率の向上を測

触ダイバータプラズマや境界層の非拡散的な動径方向輸送定した．図２０に強スケーリング測定結果を示す．MIC 上でなどに着目し，各々の現象に対応したコードの開発を進め

単一コアあたり１MPI プロセスを起動した場合，最大の６０ている［４０‐４２］．

コアまで理想的なスケーリングが得られていること，また本サイクルでは，それらの運動論的コード群のうち，特物理コアあたり複数の MPI プロセス／OpenMP スレッドに，非接触ダイバータプラズマ過程を解明するためのを起動することにより最大２倍程度のスピードアップが得

PAMCADE コー ド（Particle-in-Cell Simulation Code with

られることがわかる．しかし，理論ピーク性能に対する計Monte Carlo collision for Detachment dynamics）の開発を

算効率はCPUに比べて劣っていることや，コアあたりのメ

進めた．本コードは，Particle-in-Cell（PIC）法に基づく空モリサイズに制約があることなどの問題点があり，MIC

間１次元速度３次元の静電粒子コードに，モンテカルロ法

ノードを用いて実用的な計算を行うにはさらに検討が必要

によって様々な衝突過程を取り入れている（Monte Carlo

である．

Collision（MCC））．本サイクルまでに，すでに，プラズマ

今後，新たに開発した乱流駆動手法を用いて乱流磁気リ粒子と中性粒子の素過程（水素イオンと水素原子の弾性衝

コネクションにおけるエネルギー分配過程を明らかにする

突，及び，荷電交換，電子と水素原子の弾性衝突，水素原計算を進め，実験室プラズマにおけるプラズマ加熱現象に子の電離，及び，励起の各過程）の null-collision 法［４３］を用

対して位相混合加熱機構の検証を行う．

いた実装がなされており，本サイクルでは，さらに，プラズマ粒子間のクーロン衝突過程を Nanbu 法［４４］により導入した．
また，シミュレーションでは，計算システムの一方の境界（下流側境界）をダイバータ板に相当する粒子吸収壁，
もう一方の境界（上流側境界）を反射境界とし，上流境界側の有限幅（システム長の 1/4）の領域を粒子源領域としている．粒子源領域では，上流側の密度が一定となるように，Emmert 型の速度分布［４５］でイオン及び電子を系内に放出し，プラズマ流を生成する．そして，ダイバータ板
側の有限幅（システム長の 1/4）の領域を中性粒子領域（中性粒子の空間分布は固定）とし，プラズマ密度，温度などの時空間変動を観測した．ダイバータ板近傍に中性粒子を存在させた場合，図２１に示したように，ダイバータ板近傍における電子温度低下，イオン温度低下，プラズマ密度増加が観測された［４０］．ここで，大部分のイオンのエネル

図２
０ MIC ノードにおける強スケーリング測定結果．

ギー損失はイオンと中性粒子の弾性衝突および荷電交換衝６９Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

国際枠

第４−第５サイクル

本プロジェクトの目的は，JT-60SA，ITER，DEMO のプ
ラズマ性能予測，制御手法や運転シナリオの開発のために
統合モデリングコードを開発することである．コードの中
心は 1.5 次元輸送コードで，それに様々なモジュールを開発・結合して，統合コード全体を TOPICS と呼ぶ．
TOPICS とそれに結合するモジュール毎に，様々な並列数の計算を行う．TOPICS の中心部分は実際の放電時間を模擬し，その中でモジュールによる計算を制御する比較的軽い計算を長時間行うので機構内の PC クラスタで実行する．一方，計算の重いモジュールは Helios で実行して，遠隔にある複数の計算機で連携した統合計算を行った．
!３次元摂動磁場効果によるトルクとトロイダル回転の評価
トロイダル回転は MHD 不安定性だけでなく径電場を介して乱流輸送の抑制に寄与するので，トカマクであっても磁場の３次元性を考慮した場合に生じるトルクの評価とそれにより駆動される回転のモデリングがプラズマの予測・制御に必要である．NBI や !粒子による衝突・"!!トルク
を OFMC，摂動磁場を含む３次元平衡を VMEC で構築した 上 で 新 古 典 ト ロ イ ダ ル 粘 性（NTV）
によるトルクをFORTEC-3D（NIFS の佐竹氏・鈴木氏との共同研究）で評
価し，TOPICS の運動量輸送ソルバと径電場ソルバと連携した計算ができる統合フレームワークを用いて，ITER におけるテストブランケットモジュール（TBM）に起因するトルクとトロイダル回転を評価した．燃焼開始前の L モー

図２
１ 時刻 t ＝ 24 μs における電子温度，及び，イオン温度の空間分布（上）と電子密度，及び，イオン密度の空間分布（下）［文献［４０］Fig. 4 及び Fig. 5 より］．0.00［m］＜ x ≦ 0.05［m］が中性粒子領域，0.15［m］［m］
が粒子源領域であ≦ x ≦ 0.20
る．

ドプラズマにおいてTBMによって生じるNTVトルクが炉心深くにおいても存在し（図２２）トロイダル回転を減速させるが，GKWによる評価の結果，熱輸送に大きな影響は与えないことがわかった［４６］．

突により引き起こされており，電子のエネルギー損失はイオンと電子のクーロン衝突により生じている．さらに，ダイバータ板への熱流と中性粒子のガス圧の関係を調べたと
ころ，中性ガス圧が上がることにより熱流が逓減されることも示されている．このように，中性粒子導入によるダイバータ板への熱流逓減効果を，原子分子過程を含む第一原
理的手法により確かめることができた．
今後は，コードの最適化，高速化による計算の大規模化を進め，直線装置の実機サイズの計算が可能となるようにするとともに，非接触過程におけるプラズマ密度の逓減において重要である再結合過程の導入，さらに，非接触プラ
ズマへの高エネルギープラズマ流の影響などの運動論的コードでなければ解けない課題の解明にあたっていきたい．
１．
１１ JT-60SA，ITER，DEMO プラズマのシミュレーショ
ンのための統合モデリングコード開発

図２
２ ITER の L モードプラズマにおける NBI と !粒子による衝突・j × B トルク，NTV によるトルク（トロイダル磁場コイ
ル（TF）とフェライト鋼（FI）に加えて，真空容器の NBI
ポート（port）とテストブランケットモジュール（TBM）有無の３ケース）の分布．

JAEA integrated modeling code development for simulation of JT-60SA, ITER and DEMO plasmas（JIMS）
プロジェクト代表：林 伸彦［量子科学技術研究開発機構］
７０

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

!ダイバータ板熱負荷低減用入射不純物の炉心蓄積の評価

マイクロテアリングモードは本質的に磁場揺動を伴う電磁ダイバータ板の熱負荷低減のための外部からの不純物入的不安定性であり，微細モード構造や衝突オペレータへの射制御はITERやDEMOの定常運転シナリオでは必須であ敏感性などから取り扱いが難しかった．そこで，先行研究

り，制御時の炉心への不純物蓄積とプラズマへの影響評価の経験を持つGENEコードとのベンチマークを国際協力のが重要である．炉心不純物輸送コード IMPACT を開発し下で実施した．ASDEX-Upgrade の実験的な MHD 平衡配てTOPICSに結合し，ダイバータ統合コードSONICと連携位の下で，他粒子種衝突項なども加えた包括的なシミュ

したシミュレーシ ョ ン を JT-60SA の シ ナ リ オ に 適用しレーションを行い，構成された磁気座標や線形成長率・周た．入射不純物として Ar を用いて，その炉心での蓄積を最波数，非線形乱流輸送レベルなどがよく一致することが確

大 化 す る 計 算 条 件 で，ダ イ バ ー タ 板 熱 負 荷 を 許 容 値

かめられた．さらに，マイクロテアリングモードの線形固

2（10 MW/m ）以下にした場合の，蓄積量とその影響を評価

有関数における電流構造について詳しい解析を行い，これした（図２３）．蓄積による放射損失で炉心温度は下がるが，まで知られていた電子共鳴構造に起因する局在した電流装置の最大加熱パワー内の追加熱で回復できるレベルであ

シート構造だけでなく，トロイダルプラズマでは磁場ドリり，シナリオに対する不純物入射制御の成立性を確認した

フトによる非捕捉粒子の有限軌道幅による電流シート構造［４７］．

の変形が起こることを明らかにした（図２４）［４９］．

今後必要とされる計算として，物理トピックを限定した

次に，電磁揺動効果を含んだイオン・電子温度勾配駆動

統合モデリングによる物理解明や性能予測だけでなく，そ乱流の解析におけるマルチスケール相互作用の詳細解析を

れらが連携した計算が制御手法や運転シナリオの開発のた行った．そのために，巨大自由度の非線形相互作用の解析めに必要となる．

を実現するための流体近似解析や，異なるスケール間の相互作用を分離して評価するための部分空間伝達解析などの

１．
１２ 高ベータジャイロ運動論シミュレーション

新たな解析手法を考案した．これらを適用することで，イHigh-beta gyrokinetic simulations（HBGK）

オンスケール乱流における電子分布関数揺動が比較的短波

プロジェクト代表：前山伸也［名古屋大学］

長の帯状流を生成すること，このような短波長帯状流はイ

国際枠

オン温度勾配駆動乱流を抑制する働きを持つが，電子ス第４−第５サイクル

ケール乱流の影響を受けて減衰されやすいことなどを新た
本プロジェクトの目的は，電磁的ジャイロ運動論シミュに発見した［５０］．

レーションに基づいて高ベータプラズマにおける乱流輸送

本研究で GKV コードの適用性が確認されたマイクロテを調べることである．第４，５サイクルでは，新たにマイク

アリングモードと電子温度勾配モードの間のマルチスケー

ロテアリングモードに関するコード間ベンチマークや線形ル相互作用解析は，イオン・電子温度勾配駆動乱流以外のモード構造の解析を行うとともに，第２，３サイクルで展

系でのマルチスケール相互作用の例として興味深く，現在開したマルチスケール相互作用の詳細解析を行った．

別の計算機環境を利用して研究が進められている．これら

ジャイロ運動論シミュレーションコード GKV を用いたの知見を基に，マルチスケール相互作用のパラメータ依存［４８］．プラズマ乱流を記述する第一原理的方程式である電性や有効条件の明確化が今後の課題として挙げられる．

磁的ジャイロ運動論に基づいて，局所フラックスチューブ
モデルにおける線形成長率や固有関数構造，および，乱流
揺動の非線形時間発展を解析した．
高ベータ領域において電子熱輸送機構の候補の一つとして考えられているのがマイクロテアリングモードである．

図２
３ JT-60SA の高ベータ定常運転シナリオにおける Ar 入射時の，（a）
ダイバータ領域の Ar 密度分布，（b）
外側ダイバー
タ板熱負荷分布，（c）
炉心での価数毎の Ar 密度分布，（d）
Ar を入射しない炭素 C 不純物のみの場合と比較した電子温度と放射損失分布．

図２
４ マイクロテアリングモードにおける電子電流，静電ポテン
シャル揺動，ベクトルポテンシャル揺動の径方向・磁力線
方向分布．左列（a）‐（c）
は無矛盾なシミュレーション結果
だが，右列（d）‐（f）
では磁場ドリフト項を人為的に除いた
［４
９］
．

７１

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019１．
１３ ジャイロ運動論シミュレーションによる新古典トロ
イダル粘性研究
Gyrokinetic Simulation for Neoclassical Toroidal Viscosity（GSNTV2）
プロジェクト代表：松岡清吉［日本原子力研究開発機構（当時所属．現所属は核融合科学研究所）］
国際枠

第４−第５サイクル

核融合プラズマにおいて，プラズマ回転効果が閉じ込め性能の改善や安定性制御の観点から注目されている．実際の実験での回転分布を予測・評価する上で，コイルの誤差
磁場や共鳴摂動磁場印加等による非軸対称摂動磁場によっ
図２
５ GT5D（赤）及び FORTEC-3D（青）により得られた NTVの規格化衝突周波数依存性．橙（黒）の破線はそれぞれ磁気シア効果を含む（含まない）SBP 理論によるもの．

て駆動される新古典トロイダル粘性（Neoclassical toroidal
viscosity: NTV）が重要である．
NTV については，特に Superbanana-plateau（SBP）領!
）領域において，大域的シ域と呼ばれる低衝突周波数（!
!

!
依存性が標準的新古典理論の予ミュレーションによる !
!

測で定性的に一致しないことが示され［５１］，その不一致の原因が問題となってきた．本研究では，大域的 full- !ジャイロ運動論コードである GT5D
［５２］を用い，NTV の衝突周波数依存性の物理機構について明らかにすることを目的として研究を行った．
GT5D は大域的 full-f モデルに基づき，ジャイロ運動論方
程式を５次元空間における移流拡散方程式として解く．
GT5D において移流項は４次精度無散逸保存型差分法，衝突項には６次精度中心差分法を用いて離散化され，時間積分には半陰解法を採用している．また，本研究では静電ポ

図２６ 粒子が重水素（m = mi）及び電子（m = me）それぞれの場
合に，共鳴捕捉粒子軌道に沿って観測されたプラズマ密度
分布構造．ここで，同一の粒子軌道をポロイダル断面で見た場合，m = mi はバナナ幅が大きく，m = me ではバナナ幅
は小さくほぼ磁気面に沿った軌道となる．

テンシャルはゼロとしており，粒子種として重水素を用いた．
まず，GT5D コードによるNTV評価手法の妥当性検証とし て，［５１］で 用 い ら れ た 新 古 典 輸 送 コ ー ド で あ る
FORTEC-3D と の ベ ン チ マ ー ク を 行 っ た．GT5D と数ほどそれが顕著になることがわかる．これらから，大域!
FORTEC-3D により得られた NTV の !
依存性を図２５に示!

!
的シミュレーションで見られる NTV の !
依存性は有限の!

す．図には，標準的新古典輸送モデルである SBP 理論（磁バナナ幅から生じる位相混合による微細構造形成により生気シア効果有り／無し）［５３，５４］から得られた予測も併せてじることが明らかになった．

示されている．ここで，SBP 理論では，バナナ捕捉粒子の本研究で得られた成果をもとに，大域的シミュレーショ

トロイダル歳差ドリフト周波数がゼロとなる共鳴粒子からンに基づくより一般の NTV 評価モデルを確立するために!
には依存しないの寄与によって NTV が決定される為，!
!

は，今後は本研究で取り入れられていない径電場の効果等!
結果が得られている．図２５から，GT5D の NTV の !
依存!

についての系統的な NTV シミュレーションが必要である．

SBP 理論とは異なる依存性を示すことがわかる．次に，大１．
１４ プラズマの炉心−エッジ結合の研究!
域的シミュレーションに現れる !
依存性の原因を調べる!

Core-Edge Plasma Coupling（CEPC）

ため，粒子軌道解析を行った．共鳴捕捉粒子軌道に沿って

プロジェクト代表：宮戸直亮［量子科学技術研究開発機構］

観測されたプラズマ揺動分布構造をバナナ軌道幅の異なる

国際枠!
性が FORTEC-3D とよく一致しており，!
に依存しない
!

第４−第５サイクル

重水素及び電子質量に対して評価した結果を図２６に示す．
図から，バナナ軌道が小さな電子質量では分布は軌道に４場（プラズマ密度，磁力線方向イオン流速，静電ポテ

沿ってほぼ一定となるが，バナナ幅が大きい重水素では分ンシャル，磁場ポテンシャル）の簡約化電磁流体（RMHD）

布に有限のモード構造が生じている．このことはバナナ幅モデルに基づくグローバルシミュレーションで，過渡的なが大きい場合には，有限モード構造によって捕捉粒子領域

密度ソースを与えるとプラズマの非局所応答が観測される内に位相混合が生じることを示唆しており，実際に重水素［５５］．本プロジェクトではこれまで，この非局所応答に対での速度空間構造を観測すると，電子質量では見られない

してイオン温度勾配（ITG）駆動乱流が及ぼす影響につい

位相混合（微細構造形成）が観測され，さらに低衝突周波てのシミュレーション研究を行い，ITG 乱流が非局所応答
７２

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

を妨げる働きをすることを明らかにしてきた［５６］．非局所

Mode）振動する帯状流を観測した．帯状流の GAM 振動

応答を引き起こすには，局在化した密度ソースによって与は，主要には径電場と圧力揺動の)&'&成分の間の振動のたえられる，密度揺動の %()&成分（&はポロイダル角でトカめ，ソース／シンクを &!!で与えていた以前のシミュ

マクの外側中立面を０とする．ソースは &!!を中心に与レーションでは駆動されていなかった．もちろん乱流にえている）が必要だが，乱流がこれを散逸し非局所応答はよってもGAMは駆動され，乱流駆動のGAMでは局所的な観測されなかった．

GAM 周波数と磁力線方向のイオン音波の周波数の中間領

この結果から，乱流が弱まり，密度揺動の%()&成分に対域に周波数のピークを持つ．一方，シンク駆動の GAM は，する散逸が抑えられれば，非局所応答が起きる可能性があ

局所的な GAM 周波数にピークを持っていた．今後は乱流ると考えられる．ITG 不安定性の線形成長率は有限ベータ

輸送に及ぼす影響などを明らかにしていきたい．

効果により減少することがわかっている［５７］．そこで，上１．
１５ 磁気核融合プラズマにおけるマルチスケール乱流シ記の４場にイオン温度を加えた５場モデルによる電磁的グ

ミュレーション

ローバルランダウ流体コードを用いて，プラズマのベータ
値を変化させたシミュレーションを実施した．プラズマのMulti-Scale Turbulence Simulation in Magnetic Fusion

密度，温度，安全係数分布は前回のプロジェクトレビュー

Plasma（MS̲TURB）の報告のものと同じものを用い，ベータ値を 0.2％ から

プロジェクト代表：矢木雅敏［量子科学技術研究開発機構］

0.5％ まで変化させると，ITG 不安定性の成長率の減少に伴

国際枠

第１−第５サイクル

い，非線形飽和後の準定常状態における乱流の振幅は減少した．なお，ベータ値を 0.6％ にすると ITG 不安定性に代核融合プラズマにおいては異なる時間スケール，空間ス

わり，運動論的バルーニングモードがより不安定になる．

ケールの現象が数多く存在する．これまではスケール分離

!!!
!
#，&!!
準定常状態に達した後，$"（トカマクの弱磁を仮定し，それらの相互作用は摂動論的な取り扱いを行っ場側）を中心に局在化した，過渡的な密度ソース（と温度てきた．しかし，ある条件下では，スケール分離が仮定でシンク）を与えた．ソースを与えた場合と与えなかった場

きない場合がシミュレーションにより明らかにされつつあ合を比較すると内側コア領域での密度分布にほとんど違い

る．例えば，中立安定な長波長モードが短波長モードによはなく非局所応答は起きなかった（図２７）．しかし，ベータ

り励起される場合［５８，５９］やバースト現象［６０］のような事

値が高くなるにしたがい密度揺動の %()&成分の乱流によ例である．本プロジェクトにおいてはプラズマ乱流に注目る散逸が小さくなることを確認した．温度分布を変え乱流してマルチスケール相互作用の研究を進めてきた．特に第をさらに弱くする，あるいはソースを大きくすることで５サイクルにおいては full- "グローバルジャイロ運動論

RMHD で観測されているような非局所応答を引き起こす

コード GKNET に非断熱効果として捕捉電子の運動論的応ことができるかもしれない．

答を実装し，ITG/TEM 不安定性の線形解析を行った．

次にソース／シンクの位置（ポロイダル角 &）の影響に

full- "
ジャイロ運動論シミュレーションにおいて TEM 乱流

ついて調べた．簡単のため，電子の断熱応答を仮定し，密を取り扱うためには，運動論的電子を取り扱う必要があ度と磁場ポテンシャルを除いた静電的な３場モデルで，る．この場合，運動論的アルヴェン波の静電極限である&!%""に温度シンクだけを与えるシミュレーションを実'$ モードは通過電子の応答に起因するがこれが CFL 条件()成分ではなく)施した．これにより %&'成分の温度揺動をを律速し，計算コストが増加する問題が指摘されていた．

与 え る こ と に な る．こ の と き GAM（Geodesic Acoustic

この問題を回避するために GKNET にハイブリッド電子モ
デル［６１］を導入した．このモデルではトロイダルモード数#!
" !に対しては，通過粒子の断熱応答を仮定し，#!!
モードの寄与においてのみ通過粒子の非断熱応答を導入する．これにより '$ モードを回避する一方で新古典輸送と捕捉電子の非捕捉過程を考慮する．新しいモデルを用いて
ITG/TEM の $&スキャンを行った．ここで，イオンと電子の質量比は１００を仮定し，密度と電子温度勾配は固定した．

図２
７ 密度ソースを加えた場合（赤）と加えない場合（緑）の密度分布の比較（#= 0.5 ％）．内側コア領域（r/a<0.5）で密度にほとんど違いはない（密度ソースを加えた場合の方が中心密度はわずかに大きい）．

図２
８ ITG/TEM不安定性に対する（a）
周波数の$i 依存性および（b）
成長率の $i 依存性．R0/Ln=2.22，R0 /Ln=2.22 は固定．

７３

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

図２８に（a）実周波数，（b）線形成長率を示す．多くの文献が!
"付近で ITG から TEM へのモー示しているように，$'#"
ド変化が確認された．また，モード変化に伴いバルーニン
グ角の反転が見られた．さらに ITG/TEM 不安定性に対する磁気シア依存性に関しても調べた．磁気シアの減少に供い，成長率が減少することを確認した．これは磁気シアが小さい場合はバルーニング角が大きくなる事実と関連して
おり，非局所バルーニング理論と矛盾しない．TEM は磁気
シアにより強く安定化される．これは歳差ドリフトシアに起因すると考えられる［１３］．
１．
１６ トカマクプラズマにおける乱流拡散と帯状流抑制
Turbulent spreading and zonal flow suppression in Tokamak plasmas（TSTOK）
プロジェクト代表：矢木雅敏［量子科学技術研究開発機構］
国際枠

図２
９ 安全係数分布と各加熱源に対する定常イオン温度分布．

第５サイクルトカマクプラズマの乱流輸送の近年の理解は一義的には
フラックスと局所勾配に関連する局所パラダイムに基づいている．しかしながら，多くの実験的証拠から局所的描像には限界があることが判明している．非局所モデルによるアプローチはトカマクの乱流輸送の理解を深め，予測の確度を高めることが期待できる．トカマクの閉じ込め性能に対する乱流輸送のインパクトを考えると ITER の性能予測を行う上で最重要課題であり緊急性の高い研究課題と位置づけられる．乱流拡散は非局所輸送機構に対する第１候補のひとつと考えられる．乱流拡散により，揺動エネルギーが径方向の離れた領域に位置する異なるモード間結合をもたらし，局所モデルが与える輸送レベルを超える可能性がある．磁場の形状やトロイダル流により乱流拡散が制御できることが明らかになりつつあり，分布の硬直性を示す実験結果［６２］とよい一致を示している［６３］．本プロジェクト
で は ジ ャ イ ロ 流 体 モ デ ル［６４］や グ ロ ー バ ル #"コ ー ド
図３
０ 加熱分布モデル．シフトされたガウス分布で加熱源をモデ
ル化し，その印加場所を変化させ，イオン温度分布が定常に落ち着くまで時間発展させた．

gKPSP［６５］，グローバルfull- "コードGKNETを用いて乱流拡散やそれに関連するシミュレーション研究に取り組ん
だ．今年度は特に輸送障壁形成の問題［６６］や輸送障壁内へ

１．
１７ 複数のフラックスチューブを用いたジャイロ運動論の乱流拡散に関して研究を進めた．ジャイロ流体モデルを的シミュレーション

用い，反転磁気シアプラズマにおいて局所加熱源の位置を変化させて内部輸送障壁の形成に関して調べたのでその結

Gyrokinetic

simulation

withmultiple!"!
!
$において安全係数は最小値をとる果を報告する．#"（GMTUBES）（図２９）．加熱の大きさは固定してその場所を変化させた

プロジェクト代表：渡邉智彦［名古屋大学］!"!
!
#の熱源をおいた場合に弱い内（図３０）．その結果，#"

国際枠

flux

tubes

第４−第５サイクル

部輸送障壁が形成されることが明らかとなった．この結果
は反転磁気シアによる ITG に対する安定化効果のみなら

磁場閉じ込め核融合プラズマに生じる乱流輸送の評価おず，局所加熱源により増長された &!%シア流による安定

よび予測性能の向上には，そこで用いられるジャイロ運動化効果も輸送障壁の形成に重要であることを示している．

論的シミュレーションの精度を高めることが必要である．

一方，現モデルには加熱に伴うシャフラノフシフトの効果与えられた平衡分布における局所的なパラメータに対してやベータの上昇にともなうITG不安定性の抑制効果等が考

輸送係数の見積もりを得るために，いわゆる flux tube モデ慮されておらず，実験と比較するためにはモデルの高度化

ルを用いた解析が行われ，実験解析にも広く活用されていを進める必要がある．

る．一般に，ジャイロ運動論では，磁場に垂直方向にはジャイロ半径程度の空間スケールをもつゆらぎを想定しているが，磁力線平行方向にはトーラス大半径程度の波長を
７４

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

持つバルーニング型のモード構造を有した揺動を取り扱

力線平行方向の相関の影響評価を継続している．これまでう．磁気シアが弱い場合や線形不安定性の安定限界に近いの結果では，従来に比べ２０％程度以上の輸送フラックスの場合などに，ポロイダル角 &が $%を超えて磁力線方向の増大が認められており，現在，用いる flux tube の数に対すゆらぎが広がることが知られている．同様に，運動論的電る収束性を詳しく調べている．

子が存在する場合にもゆらぎの構造は磁力線方向に引き伸ば さ れ る．こ の た め，バ ル ー ニ ン グ 表 示 に お い て

２．高速粒子!& !&!"& に広がる揺動を扱う際には，乱流ゆらぎの２．
１ ITER１５MA 標準シナリオにおけるアルヴェ ン固有

磁力線平行方向への相関を注意深く扱う必要がある．本研

モードと LHD における高エネルギー粒子駆動測地的究では，そのために開発されたflux tube trainモデル［６７］を音響モード用いて，ITG および TEM が不安定となる場合における乱Alfvén eigenmodes in the ITER 15 MA baseline scenario and流輸送への磁力線平行方向の相関の影響を解析した．

EGAM in LHD（ITERLHD）
プロジェクト代表：藤堂flux tube train モデルは，従来の flux tube モデルの磁力

国際枠

線方向の取り扱いを拡張するもので，第１から第３サイク

泰［核融合科学研究所］

第４−第５サイクル

ルで進めた課題 GKNAXIS において開発され，GKV コードに実装された．従来の手法では，磁力線方向に長い相関長

高エネルギー粒子・MHD 連結シミュレーションコード

を扱うには，シミュレーション領域を !!%%&%!%に拡MEGA とマルチフェーズシミュレーション手法［６８］を用張していた．これに対して新しいモデルでは，flux tube

いて，トカマクと LHD における高エネルギー粒子駆動不を磁力線方向に複数接続し，バルーニング表現における安定性に関するシミュレーション研究を実施した．具体的

image mode に関する対称性を数値的にも保持したまま，な研究課題は，（１）ITER１
５MA 標準シナリオにおけるアル磁力線方向に相関長の長いゆらぎを取り扱うことができ

ヴェン固有モード（AE），（２）AE の時間発展に対する中性

る．また，&が増加することによる磁力線垂直方向波数の粒子ビーム入射（NBI）
パワーと衝突減速時間の効果［６９］，

secularな増大もなくなるため，従来の手法に比べて時間ス（３）LHD における高速イオン駆動 MHD 不安定性
［７０］，（４）テップ幅の設定においても優位性をもっている．

LHD における高エネルギー粒子駆動測地的音響モード

DIII-D 実験に対するジャイロ運動論的シミュレーション（EGAM）［７１］，（５）AE がもたらす高速イオン分布硬直性の検証でしばしば用いられるパラメータに対してこの新し［７２］である．課題（１）では，ITER１
５MA 標準運転シナリオ

いモデルを適用し，磁力線方向の相関が輸送係数に及ぼすにおいて高エネルギーアルファ粒子の有限ラーモア半径効

影響を調べた．ここではイオンとともに電子についても果を考慮すると，トロイダルモード数 "#$!程度のトロイ

ジャイロ運動論的方程式を解き，ITG/TEM 不安定性が乱ダルアルヴェン固有モード（TAE）は安定化されることを流輸送を駆動する状況を扱う．イオン熱輸送係数 #&の時間見出した．課題（２）では TFTR トカマクでの実験に基づい変 化 を，flux tube train で 用 い る flux tube の 数 !'をた物理条件を使用してマルチフェーズシミュレーションを

１，２，４と変化させた場合に得られたシミュレーション結

実行し，実験と同様の NBI パワー 10 MW，減速時間 100 ms

果を比較している（図３１参照）．ここで !'#"のケースは，の場合には複数 AE のバースト的な成長と減衰が 2 ms の時!%%&%%を扱う従来のシミュレーションと同じものであ間間隔で繰り返し発生し，実験とよく一致した結果を得る．!'##や !'#%の場合に，従来に比べて明らかに輸送た．高速イオン圧力の増大とともに高速イオン閉じ込めが係数が増大していることが見て取れる．

単調に劣化し，TFTR 実験に相当する 10 MW，100 ms の場

上述の結果を受けて，重水素‐電子プラズマに不純物（炭合のように AE バーストが発生する場合では，NBI パワー素）を加え，DIII-D の実平衡形状を用いたジャイロ運動論

や減速時間を増加しても高速イオン圧力分布が変化しない

的シミュレーションを実行し，ITG/TEM 乱流に対する磁「分布復元性（profile resiliency）」を見出した．課題（３）ではマルチフェーズシミュレーションを LHD プラズマに適用し，トロイダルモード数 "#"の複数のAEが現れることを見出した．主要な AE はプラズマ中心付近に存在するグ"$磁気面
ローバルアルヴェン固有モードと回転変換 $##付近に存在する TAE であり，両者の周波数は実験での観測値とよく対応している．
課題（４）では，HINT コードによる LHD の３次元実形状
平衡を用いて EGAM のシミュレーションを実行した．シミュレーションによって得られた EGAM の MHDポロイダ
ル流速とバルク圧力揺動の空間分布を図３２に示す．LHD
実験では第１モード（primary mode）と呼ばれる EGAMの周波数が上方に変調した後に，その周波数の半分の周波

図３
１ Flux tube train モデルを用いたジャイロ運動論的シミュ
レーションから得られたITG/TEM乱流によるイオン熱輸送
係数の時間発展．

数の第２モード（secondary mode）が突然励起されるという興味深い現象が観測されている．本研究ではMEGAを用
７５

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

す．臨界勾配は場所に依存する．
３．高速イオン圧力分布は，臨界勾配程度の勾配を保ちながらビームが入射されているプラズマ中心付近から外側へ向かって広がっていく．
４．運動論的ポアンカレ断面図解析により，輸送束の急激な増加の原因は複数の AE の 共 鳴 の 重 な り（Resonance Overlap）であることがわかった．
文献［７２］は Nuclear Fusion 誌の LabTalk 記事［７３］及び２０１６
年の年間ハイライトに選ばれた［７４］．
２．
２ トカマクにおけるディスラプション時の逃走電子生
成シミュレーション
Fokker-Planck Simulations of the Runaway Electron Generation in tokamak disruptions（FPSREG）
プロジェクト代表：奴賀秀男［核融合科学研究所］
国際枠

第１−第５サイクル

ITER や DEMO など大型のトカマク型装置では，ディスラプション時に大量の逃走電子が発生することが予想される．大電流の逃走電子ビームは第一壁に深刻な損害を与えることが懸念されており，逃走電子生成の抑制や緩和が大型トカマクにおける重要な研究課題となっている．本研究では Fokker-Planck コードを用いて，ホットテイル効果を考慮した逃走電子生成シミュレーションを行った．特に，逃走電子抑制目的で行われる不純物入射制御時におけるホットテイル効果の影響について注目した．本研究ではホットテイル効果を含めた逃走電子生成率を求めるため，
逃走電子の一次生成を電子の運動量分布関数から決定する必要がある．このため，Fokker-Planck 方程式（一次の逃走
電子生成），二次の逃走電子生成モデル，電場の拡散方程
式，熱クエンチによる温度変化モデル，不純物入射による電子密度変化モデルを自己無撞着に計算する．
不純物入射による密度変化は電子の衝突周波数を上昇さ
図３
２ （a）
LHD プラズマにおける EGAM のポロイダル流速分布と（b）
圧力揺動分布．増大した部分と減少した部分がそれぞれ青色と赤色で示されている．

せるため，逃走電子の生成を抑制することができる．しかし，一方で熱クエンチ時間を短くする傾向あるため，ホットテイル効果が発生しやすくなる．密度を変化させた場合

いたシミュレーションによってこの現象を再現することにの逃走電子生成量をホットテイル効果の有無で比較したも成功した．第５サイクル終了の時点では第２モードの振幅のが図３３である．高密度領域では，ホットテイル効果によが実験と比較して小さいことが課題であったが，その後の研究では実験と同様の結果を得ている．
課題（５）では米国の DIII-D トカマクの実験データを用いて，NBI パワーが 1.56 MW から 15.6 MW の範囲で異なる８つの場合のマルチフェーズシミュレーションを実行し，高速イオン分布と高速イオン輸送束を比較した．その結果，多数の AE との相互作用によって決定される高速イオン分布について以下の結論を得た．
１．ビームパワーの上昇とともに高速イオン閉じ込めが劣化する．
２．高速イオン圧力勾配が臨界勾配を超えると輸送束が急激に増加し，高速イオン圧力分布の硬直性をもたら

図３３ 電子密度に対する逃走電子電流量［７
５］
．

７６

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

る一次の逃走電子電流量の増加は数十 mA と非常に小さい

用い，波と粒子の共鳴相互作用を詳細に解析した［７６］．まが，二次生成によって増幅されるため，ホットテイル効果

た周波数挿引モードに対する新しい機構を提案した［７７］．

により数 MA の差が生じうることが明らかとなった．ま周波数挿引モードの長時間シミュレーションと JT-60U にた，ホットテイル効果は不純物入射による密度上昇の速さ

よる実験結果を首尾よくベンチマークすることができ，揺に依存し，逃走電子生成が始まる時点で十分に密度が上昇

動振幅の急激な脈動に対する説明を提案した［７８］．最後

していれば，逃走電子の生成は抑制できるという結果が得に，多相法を用いて ALE（図３５）を多サイクルで再現するられた（図３４）．

とともに，関連する高エネルギーイオンの雪崩的輸送を再より詳細な検証には不純物入射による密度上昇のより正現することができた［６０，７９］．

確なモデルが必要である．

ALE のシミュレーションを継続し，物理的過程，特にこれらの緩和現象のトリガー機構を同定する必要がある．シ２．
３ トカマクプラズマにおける高エネルギー粒子輸送のミュレーション結果の再現性と数値的敏感性を注意深く検大局的非線形シミュレーション

証しつつ，この方向の研究が進展している．

Global Nonlinear Simulation of Energetic Particle Trans２．
４ 垂直移動現象と逃走電子ビームのシミュレーション

portin Tokamak Plasmas（EPTRANS）
プロジェクト代表：ビアワーゲ

研究アンドレアス［量子科学技術研究開発機構］

Simulation of Vertical Displacement Event and Runaway

国際枠

Electron Beam（JVDRE）

第４−第５サイクル

プロジェクト代表：松山顕之［量子科学技術研究開発機構］
本プロジェクトはトカマクプラズマの高エネルギービー

日本枠

第４−第５サイクル

ムイオンの空間分布と速度分布関数をソース，衝突，MHD
波との相互作用を考慮することにより高精度で予測するこ本プロジェクトでは，ディスラプションがトカマク装置とを目的とする．負イオン源中性粒子ビーム入射によるに与える負荷を決める上で重要な役割を果たす垂直移動現

JT-60U 実験において常に観測されていた，突発的大事象

象（VDE: Vertical Displacement Event）のシミュレーショ

（ALE）と呼ばれる雪崩的な高エネルギーイオン輸送をシ

ン手法の研究を行った．VDEのダイナミクスを模擬する物

ミュレーションにより再現することに強い関心があった．

理モデルとしては 1.5 次元輸送モデルや２次元／３次元確立した MHD-PIC ハイブリッドコード MEGA を出発点

MHD シミュレーションが知られており，代表的なコードとして利用した．コードに現実的な N-NB イオンソースをに DINA［８０］および TSC［８１］がある．いずれのコードも組み込んだ．OFMC コードで採用された衝突モデルを実装

ITER の設計に用いられるなどトカマク装置のディスラプ

し，ベンチマークにより確認した．現実的な連続したハイション負荷や緩和手法の研究において重要な役割を果たしブリッドシミュレーション，および MHD が有／無（よりている．一方，近年の研究ではこれらの VDE シミュレー

効率的な計算）の組み合わせで構成されるマルチフェイズ

ションコードの物理モデルを拡張し，高エネルギー電子法と呼ばれるシミュレーションの両方を行った．典型的な（逃走電子）の効果や渦電流やハロー電流の非対称性の物

シミュレーションでは４０, ９６コアを使用し，数週間から数ヶ理を組み込む試みが進められており，これらの研究は月にわたって実行した．

ITER や将来の核融合炉におけるディスラプションの装置

ALE による雪崩的な輸送シミュレーションに取り組む

負荷をより精度よく予測する上で重要な課題となってい前に，JT-60U において ALE が発生する間隔で相対的に静る．本プロジェクトでは，量研で開発を進めている 1.5 次元

かな期間に周波数挿引モードのバーストとして観測されて輸送モデルと外部回路方程式を結合した VDE シミュレー

いる小さい振幅変動のダイナミクスをシミュレーション

シ ョ ン コ ー ド ETA（Evolutionary Toroidal-equilibrium

し，解析した．磁場ドリフト効果を考慮した新しい方法をAnalysis Code）の高速化・高精度化に Helios の本体システ
ムおよび増強システムを利用するとともに，逃走電子ビー
ム発生時の VDE を扱うためのモンテカルロシミュレー
ションと MHD 平衡の結合手法の開発や，ハロー電流の３次元性に重要な役割を果たすと考えられている外部キンク
モードの非線形 MHD シミュレーションなど，従来の VDE

図３
４ 一次の逃走電子生成率を密度変化の時定数で比較したもの
［７
５］
．

図３
５ MEGA コードにより再現したビーム駆動された JT-60U における３回の突発的大事象の時系列．

７７

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019の間に，主に線形 MHD モデルの拡張，および同モデルにシミュレーションを拡張するための研究開発を行った．
第４サイクルの成果として，運動論的 MHD モデルに基づいた数値コードの開発をすすめることで，以下のようよって定式化された MHD 平衡ソルバーを新たに開発し，な核融合炉の実現に重要となる MHD 安定条件の正確な予高エネルギー電子のモンテカルロシミュレーションコード

測の実現のみならず，核融合プラズマで発生しうる MHD

ETC-Rel によって評価した軸対称の逃走電子電流の２次元安定性に対する新たな物理的性質の発見を行ってきた．

分布 を MHD 平 衡 計 算 分 布 に 反 映 す る 手 法 を 開 発 し た

まず，ELM に関して，同安定性に対するプラズマ回転お［８２］．ETC-Rel は逃走電子の発生率を模擬した電場の関数

よびイオン反磁性ドリフト効果を自己無撞着に考慮した線で与えられるソース項によってテスト粒子を発生させ，そ形安定性解析を実現するために，反磁性 MHD モデルを考の後のドリフト軌道および衝突過程を追跡することで，逃案し，これを用いて線形 MHD 安定性解析コー ド MIN-

走電子のエネルギー・ピッチ角分布や空間分布を直接軌道

ERVA を拡張することで MINERVA-DI コードを開発した

損失や衝突拡散の影響を考慮して評価することが可能であ［８４］．この MINERVA-DI を用いて，JT-60U 装置および欧る［８３］．従来，逃走電子ビームプラズマの平衡磁場の計算

州 JET 装置における ELM 安定性に対する回転・イオン反には逃走電子電流分布を熱電子電流に置き換え，同時に磁磁性ドリフト効果が与える影響を解析した結果，JT-60U

気面関数として扱う仮定が用いられるのが普通であり，逃および ITER like wall を設置後の JET では，イオン反磁性

走電子電流が平衡に与える効果としては電流ピーキングに

ドリフトによる安定化効果よりも回転による不安定化効果

よるフープ力のみが考慮されてきたが，本研究で開発したが打ち勝つことで ELM が発生している可能性が高いこと

手法では図３６に例示したように逃走電子のエネルギー・を示した［８５，８６］．さらに，このような不安定化を引き起こ

ピッチ角分布を反映した２次元電流分布を考慮してビーム

すには，トロイダル・ポロイダル両方向の回転を考慮する

慣性によるシャフラノフシフトの効果を評価できる．このことが重要であることを明らかにした．

拡張は，平衡モデルの精度を向上させるだけではなく，逃また，RWM に関して，同安定性の回転分布に対する依

走電子ビームを含めたプラズマが VDE を生じるときに，存性を詳細に調べた結果，プラズマが特定の回転周波数を

逃走電子の軌道損失やビームの壁への接触に対する径方向持つ場合に新たな MHD モードが不安定化することを発見

ドリフトの効果を考慮できるため，ITER にとって重要とした．この不安定化は，RWM とプラズマに存在する安定なる逃走電子ビームの影響を考慮したポストディスラプ

な理想 MHD 固有モードの振動数が近づいた場合に観測さションフェーズにおける装置負荷の予測手法を高度化するれており，また，不安定化した MHD モードのモード構造ことが期待できる．

が，RWMではなく理想MHD固有モードとほぼ一致したことから，これら２つの MHD モードが相互作用を起こすこ

３．MHDとが原因であることを示した［８７］．さらに，RWM 安定性

３．
１ トカマクプラズマにおけるプラズマ 回 転・３次元に対する高エネルギー粒子効果について，プラズマ回転の性・運動論効果がMHD安定性に与える影響に関する影響を自己無撞着に考慮するためのモデルの考案，およびシミュレーション研究

MINERVA/RWMaC コードの拡張を行うことで，高エネ

MHD simulation for understanding rotation, three dimen-

ルギー粒子効果によって RWM がより効率的に安定化されsional effectsand kinetic effects in tokamaks（MHDRTK）

うることを発見した［８８］．

プロジェクト代表：相羽信行［量子科学技術研究開発機構］
国際枠

３．
２ ３次元平衡の線型及び非線型 MHD 解析

第２−第５サイクル

Linear and Nonlinear MHD Study of 3D Equilibria（LNMS3E）
本プロジェクトは，トカマクプラズマにおける MHD 安

プロジェクト代表：市口勝治［核融合科学研究所］

定性，特にエッジローカライズドモード（ELM）および抵

国際枠

第４−第５サイクル

抗性壁モード（RWM）の安定性に対するプラズマ回転，３次元性，および運動論効果に関する研究を進めることを目へリオトロンプラズマの MHD 安定性特性を解明するた的として進めてきた．第２サイクルから第５サイクルまで

めに，交換型モードに対するプラズマフローの交換型不安定性に対する効果を研究している．平衡及び非線型ダイナ
ミクスを３次元的に取り扱うこととしているが，現時点では，フローを矛盾なく含んだへリオトロンプラズマの３次元平衡を計算する手法はまだ開発されていない．そこで本
研究では，３次元静止平衡を利用し，そのダイナミクスの時間発展を追跡する際に，初期条件にプラズマフローを印
加することによって，その効果を調べることにした．３次元平衡及び非線形ダイナミクスの計算には，それぞれ，

図３
６ モンテカルロコードによって評価した 平 均 エ ネ ル ギ ー
20 MeV，電流値 0.5 MA，ビーム径 50 cm の逃走電子電流の２次元分布．トロイダル電流分布（右図）に関して外向
きの径方向ドリフトによる平衡磁気面（赤線）からのずれが観測される．

HINT コード［８９］及び MIPS コード［９０］を用いた．
今回，実験データを模擬したモデル的なフロー速度分布を仮定して計算を行った［９１］．ここでは，フローとしてポ
７８

Project Review（a）Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

ローのポロイダル成分とトロイダル成分が１次元データと（b）

して計測されている．これに対して，今回開発した手法を用いることによって，プラズマ全体での３次元のフロー成分が求められるようになった．
このフローが LHD プラズマでの交換型モードに対してどのような影響を与えるかを調べ，予備的な結果を得た．
（c）

この時，MHD 不安定性の効果を強調させるために，実験（d）でのプラズマよりも非常に不安定なモデル平衡を採用した．そして，初期条件にフローを加え，時間発展を追跡することによって非線型挙動を解析した．その結果，実験に対応するフローの大きさでは，プラズマの振る舞いへの影
響は小さく，フローがない場合と同様の崩壊が得られた．
これに対し，速度の大きさを実験データの２００倍にした場
合には，崩壊が少し抑制されている．このことから，プラ

図３７ ダイナミクス計算で得られた非線型領域での（a），（b）
全圧
力の鳥瞰図，及び（c），（d）
圧力の摂動成分のモードパター
ンと磁力線のポアンカレプロット．（a）
，（c）
はフローな
し，（b）
，（d）
はフローの最大値がアルヴェン速度の 0.1 倍の場合の計算結果を示す．

ズマフローは交換型モードに対して，安定化のポテンシャ
ルを持っていることがわかる．２
００倍もの大きさが必要だったのは，極めて不安定な平衡を採用したことと粘性によってフローが減衰したことによると考えられる．した

ロイダル成分のみを考慮し，その大きさが磁気面関数であがって，今後，定量的な評価のためには，実験に対応したると仮定した．また，交換型モードの存在を明確にするた平衡の採用とフロー減衰抑制モデルが必要である．

めに，回 転 変 換 が１の 有 理 面 に お い て メル シ エ 指 標 が!%!"!
$となる非常に不安定なサンプル的LHD平衡を採用

３．
４ 高ベータトロイダルプラズマの非線形 MHD 現象

した．この場合，図３７（a）及び（c）に示す通り，フローのな

Nonlinear MHD phenomena for high beta toroidal plasmas"#!#モードが成長し，非線型領域におい場合には，"!#（NMPHBTP）

いて，圧力分布の崩壊が生じている．一方，フローを加え

プロジェクト代表：佐藤雅彦［核融合科学研究所］

た場合，その絶対値が大きくなるにつれて，MHD 揺動の国際枠

第１−第４サイクル

成長が抑制されていく結果が得られた．そして，図３７（b）及び（d）に示すように，絶対値の最大値がアルヴェン速度

大型ヘリカル装置（LHD）の内寄せ配位での高ベータプ!
"倍の大きさとなった時には，モードの成長が完全にの!

ラズマにおいては，常に周辺部が磁気丘となるため，理論抑制される結果となった．

的には MHD 不安定であることが知られている．しかしな
がら，MHD 不安定であるにもかかわらず，体積平均ベー３．
３ ３次元 MHD 解析の拡張

タ値が約５％の安定なプラズマが LHD 実験では得られてExtension of 3D MHD Analysis（E3MA）いる．このことから，MHD 不安定性が線形不安定であっプロジェクト代表：市口勝治［核融合科学研究所］ても，その非線形飽和状態は閉じ込め性能に大きな影響を日本枠与えていないものと考えられる．本研究は，シミュレー

第４−第５サイクル

ション解析により，LHD 内寄せ配位における高ベータプラ
LHD 実験では，プラズマ回転の停止と交換型モードの成ズマの線形及び非線形現象を調べることで，安定な高ベー長が同期して生じる現象が観測されている［９２］．そこで，

タプラズマが LHD 実験で得られているメカニズムを明ら

ヘリオトロンプラズマにおける交換型モードに対するプラ
ズマフローの安定化効果に対する研究を行っている［９１］．
まず，へリオトロンにおいて実験データと矛盾のないフローの計算方法を検討した．磁気面量の径電場によってフローが生じていることを仮定すると，Hamada 座標系上では比 較 的 簡 単 な 表現 が 得 ら れ る．一 方，安 定 性 解 析 には，３次元 MHD 平衡計算コード HINT
［８９］と３次元 MHD
ダイナミクスコード MIPS
［９０］を用いる．これらのコード
では，円柱座標が採用されているので，円柱座標への座標変換が必要となる．さらに，実験データを導入する手法を組み合わせることによって，３次元フロープロファイルを求める計算手法を確立した［９３］．図３８に，LHD プラズマで
図３
８ LHD での１次元実験データに基づいて計算した３次元フローの流線図．

計測されたデータに基づいて計算したフローの流線図を示す．LHD 実験では，一つのポロイダル断面の水平線上でフ
７９

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

かにすることが目的である．

本研究で得られたシミュレーション結果では，コア崩壊

本シミュレーション解析では，一流体モデルを数値的にが引き起こされており，高ベータプラズマを安定に維持す

解く MIPS コード［９４］を用いて行った．MIPS コードでは，る結果が得られていない．その原因の一つとしては，磁気円柱座標系のもとで４次精度の中心差分法を用いて離散化

レイノルズ数が実験値よりも小さい値を用いてシミュレーを行っている．本研究では，非線形計算の数値安定性のたションを行ったことが挙げられる．本研究で見られた不安

めに，移流項に３次精度の風上差分法を適用した修正を定性は抵抗性バルーニングモードであることから，実験値行っている．時間積分は４次精度のルンゲ・クッタ法を用に対応した高磁気レイノルズ数を用いたシミュレーション

いて行われている．LHD プラズマの MHD 平衡では，最外では，本研究で得られた結果よりも MHD 不安定性の成長

殻磁気面の外側に磁力線がストキャスティックな領域が存が穏やかになるものと推測される．したがって，高磁気レ在し，ベータ値が高くなるほど磁力線がストキャスティッ

イノルズ数では，本シミュレーションで見られたコア崩壊

クな領域がコア側へ広がってくる性質がある．そのため，が抑制される可能性がある．今後は，実験値に対応した高磁力線がストキャスティックな領域も含めた計算を行う必

磁気レイノルズ数に対して計算を行うことにより，飽和レ

要がある．本研究では，磁気面の存在を仮定しない HINT

ベルが低下してコア崩壊が抑制されるかを確かめていく必

2 コードを用いて，中心ベータ値が 7.4％ と 9.2％ の２つの要がある．

MHD 平衡を構築し，磁力線がストキャスティックな領域も含めてシミュレーションを行った．モード解析を行うた

３．
５ 高ベータトカマクプラズマにおけるフィッシュボー

めに，最外殻磁気面の内側の領域に対して，VMEC コード

ンモードのようなシミュレーションのための拡張ハを用いて MHD平衡の再構築を行った後，NEWBOZコード

イブリッドモデル開発を用いて最外殻磁気面の内側の領域での Boozer 磁気スペ

Development of Extended Hybrid Model for Simulation of

クトルを取得し，揺動成分のモード分解を行っている．

Fishbone-Like Modes in High-Beta Tokamak Plasmas

中心ベータ値が 7.4％ の時の縦長断面における圧力の等（MEGAKTI）

値面の時間発展の様子を図３９に示す．線形段階では周辺部

プロジェクト代表：ビアワーゲ

にてMHD不安定性が発生する．このMHD不安定性のモー技術研究開発機構］

ド構造，および，線形成長率の磁気レイノルズ数依存性の日本枠

アンドレアス［量子科学第４−第５サイクル

解析から，不安定モードは抵抗性バルーニングモードであることがわかった．非線形段階になると，不安定性の影響

本プロジェクトは，イオンランダウ減衰のような熱イオ

がコア領域まで広がり，最終的にはコア崩壊が引き起こさンの運動論効果（KTI）を含む大域的非線形 MHD-PIC ハイ

れる結果が得られた．このコア崩壊メカニズムについて詳

ブリッドコードの開発をめざしている．これはベータ誘起細な解析を行ったところ，異なるトロイダルモード数を持

アルフベン固有モード（BAE）ギャップの周波数領域におつバルーニングモード間のモードカップリングにより，いて，高エネルギーイオンと低周波モード間の相互作用を

!!
"#
!
!
!#を含む低次のモードが発生したことが原因"!"

シミュレーションし，活性化された MHD モードによるバであることがわかった．特に，"
!!
"#
!
!
!#のモードは磁!"ルクイオンの加熱率を定量化するために必要である．

気軸にて有限の速度を持ち，コア領域のプラズマを外側へ確立した MHD-PIC ハイブリッドコード MEGA を出発点シフトさせるように働く．コア領域に見られる揺動成分はとして利用した．高エネルギーイオンは４点補間ジャイロ

非線形的に励起された低次モードによるものであり，コア平均を考慮したドリフト運動モデルを利用してシミュレー

領域にて低次モード自身のモードカップリングが発生する

ションを行った．バルクプラズマのダイナミクスは抵抗とことで，コア崩壊が引き起こされることがわかった［９５］．

粘性拡散を含む完全 MHD モデルを使いシミュレーションを行った．典型的なシミュレーションでは４，０９６コアを使用し，数日間にわたって実行した．
拡散係数と他のパラメータに関する MEGA シミュレー
ションの感度を調べた［９６］．研究の結論として，現在利用可能な計算資源では本来の目的に対し，熱イオンのジャイ
ロ運動的取り扱いは可能でないことが判明した（［９６］の付
録 B 参照）．主な理由として，熱イオンの有限ラーマ半径効果をとらえるためにはたいへん高い空間分解能が必要であること，抵抗スケールを有限ラーマ半径スケール以下にしなければならないことがあげられる．有限ラーマ半径効果を無視し，現在用いている空間分解により磁場に沿った熱イオンの流れの運動論的圧縮性のみ考慮したシミュレー
ションを行ったが，抵抗散逸は運動論的熱イオン効果と競

図３
９ 中心ベータ値が "0 ＝ 7.4％ の MHD 平衡に対する非線形シ
ミュレーションより得られた（a）
非線形初期段階と（b）
飽和状態での縦長ポロイダル断面における圧力の等値面．

合し，それらを阻害する細かい空間スケール及び低周波数
領域で主要な役割を演じ続ける．一方，発見はこのプロ
８０

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

ジェクトに新しい展望を与えた．MEGA による高ベータ

BA Project Planning and Coordination Group et al.

［１０１‐１０３］は以下の通りである．

JT-60U プラズマのスペクトル分析により，BAE ギャップ

拡張された MINOS コード（MUTSU/MINOS コード）に（図４０）［９７］の十分上の周波数でアルフベン波と音波結合がは，Miura et al.［１０１］で一様等方非圧縮性 Hall MHD 乱流お明らかになった．そのような結合は燃焼核融合プラズマによび磁化された一様非圧縮性 MHD 乱流のために開発されおいて，熱イオンの非衝突加熱に対する未知のチャンネル

た Sub-Grid-Scale（SGS）モデルを用いて，LES を実行するを提供する可能性がある．この発見は新しい２つの研究活機能が付加された．このため，圧縮性を伴う拡張 MHD 方動のきっかけとなる．線形ジャイロ運動シミュレーション

程式にFavre平均を施したGS方程式を導出し，速度による［９８］を利用したアルフベン波と音波結合の物理を研究する

圧力の移流を表す項に起因する SGS 項を新たに用意した．

国際共同研究を行うこと，及び MEGA を利用し，バルクイ

MUTSU/MINOS コードを用いて，LHD 中のバルーニン

オン の 異 常 加 熱 の 実 験 的 証 拠 を 探 査 す る た め JT-60SA

グ不安定性の非線形時間発展の LES を行った．この問題の

［９
９，
１
０
０］
における高ベータシナリオを研究することである．

LES は初めてであるため，過去のデータ［１０４］と計算結果

運動論的熱イオン効果のシミュレーションは，新しい高を比較しつつ，コードおよび LES 機能の確認を行った．そ性能計算機資源が利用可能になったとき，再考される．その典型的な結果を図４１（a）に示す．図４１（a）はバルーニングの間，線形ジャイロ運動シミュレーション［９８］の利用，お不安定性の線形段階終盤における圧力の等値面・等高線とよび MEGA［９９］を使った高ベータ JT-60SA の研究に集中

流線の，ポロイダル断面近傍での拡大図である．（b）
は不する．

安定性の飽和段階における圧力の等値面・等高線と磁力線（太管），流線（細管）である．どちらも文献
［１０３］
の図４か

３．
６ 磁化プラズマの圧力駆動型モードの拡張MHDシミュら引用した．このデータの詳細解析から，LES は線形段階

レーション

においては比較的低波数の線形不安定性を良く解像しつHigh Accuracy Extended MHD Simulation of Pressure

つ，少ない格子点数でも発散することなく非線形飽和を実Driven Modes in Magnetized Plasmas（HMHDPDM）

現するなど，LES の有用性を示した．
今後は，ジャイロ粘性項と LES の関係を確認した上プロジェクト代表：三浦英昭［核融合科学研究所］
国際枠で，LHD など短波長不安定性が成長し得る不安定性の非線

第４−第５サイクル

形過程について，SGS モデルの改良を行い，不安定性解析に LES を用いる予定である．

本研究の目的は，１）圧力駆動型モードなど磁化プラズ
マに現れる短波長不安定性の線形・非線形発展と非線形過程における短波長成分の役割を解明すること，２）短波長成分の役割を数値モデルで代替し，方程式をグリッドス
ケール（GS）のみ解く Large Eddy Simulation（LES）を実現することである．
こ れ ら の 目 的 を 達 成 す る た め，大 型 ヘ リ カ ル 型 装置
（LHD）中のプラズマについての MHD 方程式シミュレー
シ ョ ン コ ー ド MHD In Non-Orthogonal System（MINOS）
コードを，Hall 項およびジャイロ粘性項を含む拡張 MHD
方程式に拡張 し た．こ の コ ー ド を 用 い た 主 な 研 究成果

図４
１ LHD の磁気軸位置 Rax ＝ 3.6 m，!＝ 3.6％ 平衡についてのバルーニング不安定性の LES 結果．（a）
線形段階終盤における圧力の等値面・等高線と流線の，ポロイダル断面近傍での拡大図．（b）
不安定性の飽和段階における圧力の等値
面・等高線と磁力線（太管），流線（細管）．どちらも文献［１０３］
図４から引用．

図４
０ スペクトル分布はアルヴェン波と音波結合を示している．

８１

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019３．
７ 統合コードを用いた ITER ディスラプションのシミュ
レーション研究
Simulation Study on Disruption in ITER Plasmas using Integrated Disruption Code（DISRUPT）
プロジェクト代表：矢木雅敏［量子科学技術研究開発機構］
国際枠

第１−第５サイクル
図４
２ EXTREMコードによる逃走電子電流分布形成のシミュレー
ション：抵抗性キンクモードを考慮しない拡散型シミュ
レーション（a）
では逃走電子発生に伴い電流分布がピーキングするのに対し，非線形 MHD シミュレーション（b）
では抵抗性キンクモードの励起が繰り返されて電流分布はマージナルに保たれる．

トカマク型核融合炉の研究においてディスラプションは装置の健全性を脅かす難題であり，ITER においては熱負荷緩和と逃走電子発生回避の両立が信頼性の高いディスラ
プション緩和システム構築の最大の課題となっている．本
プロジェクトでは ITER におけるディスラプション時の逃
走電子発生メカニズムの研究を主課題とし，トカマクプラてはじめて明らかにした［１０５］．開発したシミュレーショ

ズマのディスラプションや過渡現象に関する理論・シミュ

ンコードEXTREMは，逃走電子発生時のMHD不安定性のレーション研究を高度化するための数値計算手法，物理モ

励起とそれに伴う逃走電子の輸送，逃走電子発生率への影

デル，並列計算技術に関する研究開発を行った．第４‐第５

響を自己無撞着に解析できるので，今後，不純物の効果なサイクルでは，新たに簡約化 MHD コードを用いた逃走電

どを考慮して発生メカニズムに関するモデルをさらに詳細

子発生現象の研究［１０５］，抵抗性 MHD 不安定性と ITG モー

化することで ITER における逃走電子発生機構の予測と緩ドを統一的に記述する６場簡約化 MHD シミュレーション

和手法の開発に貢献することが期待できる．

を用いた非局所輸送のシミュレーション研究［１０６］，逃走電子モンテカルロシミュレーションコードETC-Relの開発

４．周辺プラズマ［１０７］などが進展した．

４．
１ 共鳴摂動磁場を加えたトロイダルプラズマの粒子・エネルギー輸送解析

このうち，第４‐第５サイクルの主要な成果として，抵抗性 MHD 不安定性を扱う簡約化 MHD コード EXTREM にParticle and Energy Transport Analyses of toroidal plasma

ディスラプション時の逃走電子発生を模擬するように工夫

affected byresonant Magnetic Perturbations（PETAMP）

した高エネルギー電子（逃走電子）のモデルを実装するこ

プロジェクト代表：菅野龍太郎［核融合科学研究所］とで，ディスラプションによってプラズマ冷却が開始され国際枠

第４−第５サイクルてから定常的な逃走電子ビームが形成されるまでの数１０ミ
リ秒の長時間にわたる非線形 MHD シミュレーションを実

近年，トカマクプラズマにおける ELM の制御法の１つ現したことが挙げられる［１０５］．従来，逃走電子の分布形

として，共鳴摂動磁場（RMP）の利用が試みられている成のシミュレーションは１次元拡散方程式に基づく計算が［１０８］．十分強い RMP を静的に加えたプラズマエッジ領域

ITER の予測として議論されてきたが，本研究でははじめ

では磁力線はエルゴダイズするが，衝突度 #!〜0.04−0.12て MHD 不安定性による電流分布の大域的な変化を含んだの場合の径方向輸送係数が，従来の標準的な輸送モデルで形での分布形成のシミュレーションを可能とした．図４２ある Rechester-Rosenbluth（以下，R-R）の理論［１０９］に基は，（a）MHD 不安定性を考慮しない従来の拡散方程式にづいた予測と大きく異なっている［１０８］．本課題では，第よるシミュレーションと，（b）今回開発した非線形 MHD

１〜３サイクルに引き続き，このような領域におけるプラ

コードによるシミュレーションで得られた逃走電子電流分

ズマ輸送解析に必要となるコードの開発を行った．ここ布の時間発展を比較した結果を示したものである．このシで，計算手法としては，前サイクルまでの成果に基づき，

ミュレーションでは，逃走電子発生時に起こる不安定性とドリフト運動論的 "［１１０］を採用した．
!シミュレーション

して，逃走電子発生に伴う電流ピーキングによって中心の磁場の３次元性が強い場合，電場は，両極性条件によって安全係数が１を切った時に生じる ""!抵抗性キンクモー

決まる．グローバルにイオンと電子のドリフト運動論方程ドを考慮した．図４２（b）は逃走電子のビーム電流分布は鋸式を解きつつ，同時に電場を求めることは，極めて多くの歯状振動に類似した電流分布の飽和現象を伴うことで図４２計算資源と計算時間を必要とし，その実行は現実的でな（a）にみられるような電流ピーキングが阻害され，抵抗性い．そこで，３次元磁場配位における新古典輸送解析におキンクモードがマージナルに安定（##!）になるような形

いて，しばしば仮定されるように，時間固定したモデル電状を保ったままビーム電流分布が形成されることを示して場を与え，その下で案内中心分布関数 "!を 求 め，条 件

いる．このとき，単に逃走電子電流の分布が変化するだけ!#"!$から両極性電場を決定する手法を開発した．ここ

ではなく，（１）ビームのインダクタンスが低く保たれるこで，!#および !$は，それぞれ電子およびイオンの動径方向とで総ビーム電流値が間接的に増加する働きがあること，

粒子フラックスである．手法の詳細については，［１１０］を（２）逃走電子ビーム電流発生の初期において抵抗性キンク

参照．本課題における "!シミュレーションでは，０次オー

モードが与える強い誘導電場が逃走電子量を直接的に増加

ダーの分布としてマクスウェル分布 !
" を仮定している．実させる働きがあることを本研究のシミュレーションによっ験では，径方向熱拡散係数が R-R 理論よりもオーダー小さ８２

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

く，温度分布は完全に平坦化することなく勾配が残ってい４．
２ 原型炉設計に向けたシミュレーション研究

ることが知られているが，磁場が乱れた領域において勾配

Simulation Study for Fusion DEMO Reactor Designのある温度分布を仮定した "#シミュレーションが，その仮（SSDEMO）

定に反さずに実行できるのか，確認が必要である．そこで，プロジェクト代表：星野一生［量子科学技術研究開発機構］

まずは電場無しの条件下で，数密度の勾配はゼロ（つまり，国際枠第１−第５サイクル!('%$一定）だが，有限の勾配をもつ温度分布（これを "とする）を設定し，電子の案内中心分布関数 "#から，"'!(!('

からの温度のずれ "" を求め，"!""('

本プロジェクトの目的は，原型炉設計研究における重要

!(
'

と"
"
$" ""
を比較した．ここで，"'!(は，RMP を無視した場合の磁気課題及びそのオプション検討のための数値シミュレーショ

面 を リ フ ァ レ ン ス 面 と し て，そ の 関 数 と 仮 定 し，最重要課題のひとつであるダイバータの検討や，炉心プラ，+#,は 隣 り 合 う リ
!($')
$)'#!
"
#
"
"$'
!
%(
$(
"
'
,
!"'!(*+
ファレンス面の囲む微小体積による平均である．温度分布

ズマへの粒子供給の検討，保守概念検討のための炉内線量ン研究である．全サイクルを通して，原型炉設計においてマップの解析，ブランケット第一壁設計検討等の数値シ!"に対して，磁気島同"'!(のスケール長 !' &$,$&-"'!(!
士 の 重 な り 幅 !/ が 同 程 度 の オ ー ダ ー に 近 づ く と
!""(
'

（!/!
!' . "），"

ミュレーション研究を進めた．以下では，最も多くのリソースを使用したダイバータ検討の結果について述べる．

!('

は " から大きくずれてしまい，"
#
成
立
し
な
く
な
る（図
４
３
）
．一
法の前提条件"
#
!
#
%"が
&

シミュレーションコードSONICを用いた．ポロイダル断面方，磁場はそのままで，"'!(を !/!
!'"1/10 となるようににおける２次元の SOL・ダイバータプラズマ輸送は流体

勾配を緩めた温度分布に変更して計算すると "" が小さくモデルにより扱い，中性粒子及び不純物の輸送はモンテカ

なることがわかった．ここで，エルゴダイズした領域におルロ粒子モデルにより解いている．典型的な原型炉ダイ

けるチリコフ（Chirikov）パラメータは #%*+.!1 である．以

バータの解析では，１ケースあたり３２ノード（５１２コア）を上より，エルゴダイズした領域において勾配がゼロでない用いて１５０時間程度の計算リソースを使用する．

原型炉ダイバータプラズマの解析には，統合ダイバータ

温度・密度分布を仮定しても，"#シミュレーションが実行

原型炉ダイバータにおける最重要課題は，主プラズマかできることがわかり，今後，電場の効果を含め，輸送特性

ら排出される膨大な熱の制御である．主プラズマから排出の理解を進める予定である．

される熱は ITER の数倍に及ぶ一方，ダイバータでの除熱
は，中性子照射環境における材料選択の制約から ITER ダ
イバータと同程度かそれ以下となる．第１〜第３サイクル
では，主に原型炉概念 SlimCS（核融合出力 3 GW，装置サ
イズ 5.5 m）を対象とし不純物ガス注入による放射損失と非接触ダイバータプラズマを併用した熱制御シナリオの検討［１１２，
１１３］を進めると共に，核融合出力と不純物放射量をパラメータとした運転領域の解析［１１４］を行った．
これらのダイバータ物理検討に加え，原型炉設計全般の検討結果を踏まえ，原型炉の設計方針が核融合出力 1.5 GW
程度，装置サイズ 8 m 級に見直された［１１５］．装置の大型化は受熱面積増加による熱負荷低減が期待できる．しかし，現在の設計案では主プラズマ密度が 6.6×1019 m−3 と低
く，この結果，熱負荷低減に必要な非接触ダイバータプラ
ズマの形成が不利になることが懸念された．そこで，新しい原型炉概念案を対象に，まずは周辺領域への燃料ガスパ
フを行わず，可能な限り SOL 密度を低くし，ダイバータプ
ラズマ性能予測と熱制御シナリオの再検討を行った［１１４］．SONIC 解析の結果，1.4×1019 m−3 と，現案の主プ
ラズマ密度と整合しうる低い赤道面 SOL 密度が得られた．この時のダイバータ熱負荷分布を図４４に示す．核融合
出力が以前の 3 GW から半分の 1.5 GW になったことに加え，この解析では周辺領域へ排出される熱の８０％を Ar 不

図４３ （a）
Lw/LT→１になると，電子温度のずれ "T により，温度分
布は平坦化する（図中の破線）．ここで，r/a は規格化小半
径．一方，（b）
Lw/LT ! 1/10 となるよう，（a）の破線の r/a〜0.8 における勾配とほぼ同じである温度分布に変更して計算すると，"" は小さくなる．ここで，両ケースとも，エルゴダイズしている領域は，0.6 ＜ r/a ＜ 1 の範囲である．また，RMP 無しの磁場として，円形トカマク磁場を仮定した
［１
１
０，
１
１
１］
．

純物により放射させているため，低 SOL 密度であっても部分非接触ダイバータが形成され，熱負荷は当面の目標である 10 MW/m2 を下回る結果が得られた．したがって，低SOL 密度でも十分に不純物放射させることで，熱制御シナリオが成立できることを示した．しかし，この時のダイバータにおけるピークイオン温度は 300 eV を超えており，
８３

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

燃料プラズマや Ar 不純物による損耗が問題となる可能性

３要素により L-H 遷移が記述される．本プロジェクトでが高い．この解決策の一つとして，100 Pa m3/s の D2 ガス

は，電磁プラズマ乱流系における L-H 遷移シミュレーショ

19−3

パフを行い，赤道面 SOL 密度を 3×10 m

程度まで増加ン手法の構築の前段階として，前述の２場の L-H 遷移モデ

させることで，ダイバータ板のピークイオン温度を 100 eV

ルを三次元周辺 MHD コード BOUT++［１１９］へ実装して程度まで減少させることができる結果を得た．得られたダ

L-H 遷移シミュレーションの予備計算を行うとともに，３

イバータプラズマ分布の元で，正味の損耗量の評価を行う

場の簡約化モデルによる無衝突バルーニングモード

ために，IMPGYRO コードを用いてダイバータの損耗・再（CBM）／抵抗性バルーニングモード（RBM）不安定性の堆積の解析も進めている［１１６］．今後，SONIC コードの予

線形成長率に対する抵抗率と電子表皮長のパラメースキャ

測精度向上へ向けた改良を進めると共に，原型炉ダイバーンを実施した［１２０］．

タプラズマの物理特性の理解を深めながら，炉心プラズマ

図４５は L-H 遷移シミュレーションの予備計算で得られた設計およびダイバータ工学設計と整合したダイバータ物理

平均流シアと乱流強度の径方向分布であり，周辺部で平均設計の最適化を進めていく．

流 シ ア に よ る 乱 流 強 度 の 抑 制 が 確 認 で き た．一 方 で，
BOUT++コードの数値スキームの検証から，フルート近４．
３ トカマクプラズマにおける L-H 遷移のシミュレー

似を用いた二次元ポアソンソルバーは揺動の平均流成分のション研究

ポロイダル構造を正確に表現できないため数値振動を引きSimulation Study on L/H Transition on Tokamaks

起こすことが明らかになった．周辺輸送障壁形成時に平均（LHTOK）

流圧力と平均流磁場から平衡磁場を更新する必要がある電プロジェクト代表：矢木雅敏［量子科学技術研究開発機構］

磁プラズマ乱流系における L-H 遷移シミュレーション手法

国際枠の構築には，平均流成分の整合良い取り扱いが不可欠であ第４−第５サイクル

るため，本プロジェクトで得られた知見を基に平均流成分
プラズマ境界近傍の周縁領域に輸送障壁が形成され，炉を整合良く取り扱うことが可能な三次元ポアソンソルバー

心の閉じ込め性能が改善する H モード放電は ITER・原型の開発を行っている．

炉の標準運転モードであり，L-H 遷移の物理機構の解明は

図４６は 円 形 断 面 の 解 析 的 平 衡 に お け る #"#!-CBM/

ITER・原型炉の炉心性能予測の高度化における重要課題

RBM 不安定性の線形成長率の電子表皮長・抵抗率依存性の一つである．これまでの L-H 遷移の理論研究としては径を表している．実電子表皮長を用いた場合は，電子温度ペ

電場とプラズマ回転シアが周辺乱流を抑制して L-H 遷移に

2
3）
デスタルに特有の電子温度領域（10［eV］!!"!10［eV］

至るモデルが複数提案されているが，実験を定量的に説明に相当する抵抗率 "!!%!#!"!!$で RBM/CBM 不安定性

出来るモデルの確立には至っていない．そのため，従来のの遷移が起こるが，現実よりも大きな抵抗率や電子表皮長

L-H 遷移モデルと実験との直接比較による検証に加えて，を用いた場合は局所理論では十分に小さいと近似されていより第一原理的な物理モデルに基づく L-H 遷移の数値シる磁力線方向の構造の影響が無視できなくなるため，電子ミュレーション手法の構築が課題となっている．

表皮長・抵抗率依存性は弱くなる結果が得られた．これは近年，熱源駆動型の２場の簡約化 MHD 方程式に基づく

非理想効果を導入した実形状平衡の周辺安定性解析・乱流

L-H 遷移モデルによる数値シミュレーション
［１１７，
１１８］が

計算を行う場合は流れの相似性を仮定せずに現実の抵抗率報告されており，これらのモデルでは，静電抵抗性乱流，と電子表皮長を用いる必要があることを示唆する結果であ

渦度方程式に新古典粘性に由来するダンプ項として導入さる．

れる径方向の力の釣り合い，圧力の式のソース・シンク項と乱流輸送の釣り合いとして導入される動的圧力分布，の

図４５ ２場の L-H 遷移モデルの予備計算で得られた平均流径電場―&2〉
Er と乱流強度〈"の径方向分布．

図４４ 核融合出力 1.5 GW の 8 m 級原型炉におけるダイバータ熱
負荷分布［１１
４］．

８４

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

に加え，背景プラズマに温度勾配及び密度勾配が存在する
場合に生じる背景プラズマの速度分布関数のゆがみを考慮し，モンテカルロ二体衝突モデルを適用した．これにより，
磁力線平行方向はもちろん，従来，SOL／ダイバータプラ
ズマの運動論的シミュレーションでは，全く考慮されてこなかった磁力線垂直方向の古典・新古典輸送による不純物輸送の体系的な理解を行うためのモデリング基盤が，ほぼ整備されたと考えられる．
（２）JT-60U を対象とした計算メッシュ体系におけるW輸送
シミュレーション［１２６］成果（１）に基づき，IMPGYRO コードに磁力線平行方向の運動論的熱力の効果を導入し，JT-60U 計算メッシュ体
図４６ n ＝ 40 − CBM/RBM 不安定性の線形成長率の電子表皮長・抵抗率依存性，ここでde*/de は現実の電子表皮長で規格化された電子表皮長の大きさを表す．

系での不純物輸送シミュレーションを行った．計算の結５．統合モデリング

ル面上の不純物分布の上下非対称性が生じると，磁気ドリ

５．
１ SOLPS-IMPGYRO 統合コードによる不純物輸送シ

フト由来の内向きの粒子束が発生し，SOL からコアへと不果，不純物は熱力・摩擦力による磁力線方向の輸送によって SOL 上部に滞留する場合がある．さらに，このポロイダ

ミュレーション

純物が浸入することが示された．本結果は，成果（１）で述べImpurity transport simulation of the ITER-like plasma by

た新古典輸送効果の重要性を一部示唆している．

the coupled SOLPS-IMPGYRO codes（IMPGYRO）（３）新古典輸送の効果の検証［１２７］

プロジェクト代表：畑山明聖［慶應義塾大学］
国際枠

成果（２）
で得られたコア浸入過程の詳細な検証を行うた

第５サイクル

め，磁力線垂直方向の W 輸送過程の解析を行った．その結果，新古典的な輸送が，従来 SOL／ダイバータ解析で用い

ITER およびそれに続く原型炉では，タングステンがダられてきた異常輸送係数による異常輸送と比較し無視できイバータ板および第一壁候補材になっている．炉心プラズない領域が存在することがわかった．これは，SOL におけ

マの性能は，炉壁から発生する不純物によって大きく影響る新古典輸送の効果が無視できないことを示唆している．

を受ける．したがって，タングステンのような高 Z 材料の（４）IMPGYRO コードによる ITER 運転シナリオでの W

周辺プラズマでの輸送過程およびコアプラズマへの混入過不純物輸送予測計算［１２８，１２９］程の理解とその制御手法の確立は，ITER および原型炉設

IMPGYRO コードを初めて ITER に適用し，現在考えら計，ひいては，核融合実用化にとって必須の課題と考えられている代表的な運転モードにおける W 不純物輸送の予れる．本研究は，ITER，JT-60SA 等の予測シミュレーショ

測計算を行った．日本のコードとしては，ITER の SOL／

ンが可能な信頼性の高い不純物輸送シミュレーションコー

ダイバータの不純物輸送解析ははじめての試みである．

ドを開発することを目的とする．

ITER プラズマ計算に用いたメッシュを図４７に，計算結果

慶應義塾大学にて開発されたモンテカルロ高 Z 不純物輸

例として得られた W 不純物密度分布を図４８に示す．IMP-

送コード IMPGYRO［１２１］は，１）不純物粒子のジャイロ

GYRO による ITER プラズマにおける W 不純物輸送の計算

運動を直接計算，２）背景イオンとのクーロン衝突を二体基盤が整い，ITER への貢献が可能となった．今後，H モー

衝突法により模擬，３）モンテカルロ法による多価電離・ドやドリフト等を考慮した詳細な計算を行う上での，リ再結合過程のモデル化，４）EDDY コード［１２２］との結合ファレンスケースが確立できた．

によるスパッタリング・反射機構を考慮している等の特徴（５）SOLPS-ITER IMPGYRO 統合コードの開発を持っている．以上の特徴から，IMPGYRO はダイバータ

正確な W 不純物のコアへの混入量評価のためには，自己領域，スクレイプ・オフ層（SOL）における大域的な輸送

無撞着な統合コード開発が必要である．コード開発にあたのみでなく，壁表面でのプロンプト・リディポジションやり，ITER 機構サイトにて６ヶ月間，本プロジェクトの主

セルフスパッタリングの効果をも考慮可能である．IMP-

要メンバーである ITER 機構 Xavier BONNIN 博士と慶應

GYRO のような運動論的輸送コードは並列化計算との相性

義塾大学博士課程学生

がよく，Helios を利用した場合，128 CPU の本研究室のク

SOLPS-ITER IMPGYRO 統合計算コードの開発に成功しラスターシステムを使う場合と比較し，計算時間を 1/10

た．これにより，背景燃料への W のイオン化や放射冷却の以下に抑えることが可能である．

影響を考慮した，ELM を含む詳細な実験解析・予測計算

矢本昌平氏が共同で作業を進め，の基盤が確立された．

主として以下の成果を得た．ページの制約から，詳細にH モード時のペデスタルを考慮した背景プラズマ分布，

ついては参考文献に譲る．

および背景プラズマのドリフトを考慮した計算を行い，よ（１）タングステン輸送シミュレーションモデルの高度化とIMPGYRO への実装
［１２３‐１２５］

り詳細な W 輸送過程の解析を進めていく．また，特に，モデルの高度化では上で述べた従来のIMPGYROモデル

ITER での ELM 発生時の W 輸送についても，開発された
８５

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019（a）部分デタッチ状態（a）全体メッシュ（b）ダイバータ領域拡大図
図４
７ ITER プラズマ計算に用いたメッシュ形状［１２９］．
（b）高リサイクリング状態

SOLPS-ITER IMPGYRO 統合コードを利用し，計算を行っ

図４
８ 計算結果例．ITER 形状での W 不純物密度分布［１２９］．

ていく予定である．
５．
２ 統合ディスラプションコードのためのモジュール開発

ギーを計算するサブルーチン tmenr（FFTW ライブラリを

Module Development for Integrated Disruption Code使用），セットしたマトリックスをLU分解により解くサブ（MODEV）

ルーチン tmpus（LAPACK ライブラリを使用），非線形項

プロジェクト代表：矢木雅敏［量子科学技術研究開発機構］

を擬スペクトル法を用いて評価し右辺ベクトルにセットす国内枠

るサブルーチン tmrhs（FFTW ライブラリ（スレッド並列

第２−第５サイクル版）を使用）．tmrhs における経過時間のスレッド数依存性
トカマク型核融合炉の研究開発において，ディスラプを調べたところスレッド数１に比べ，スレッド数８の時のションは装置の健全性を左右する問題であり，その回避法

ほうが，FFTW のスレッド並列効果により，性能向上が見や緩和法の開発は急務とされている．本プロジェクトは計られた．参考のため，Helios 本体システム（インテル社製

算科学の観点から，ディスラプション統合コードに関連す

CPU の Sandy Bridge を搭載する）でも同様の計測を行っる物理モジュールの開発やコードの高速化を進めた［１３０‐た．傾向は増強システムと同様であったが，計測のばらつ

１３２］．特に第５サイクルにおいては，Helios 本体システム

きは小さかった．Elapse 時間は Helios 本体システム 26.75

用にチューニングした簡約化 MHD モデルを Helios 増強シ秒，Helios 増強システム 234.34 秒であり，約 8.8 倍程度の性

ステム（インテル社製 MIC コプロセッサ Knight Corner

能差が確認された．また，欧州の High Level Support Team

（KNC）を搭載する）に移植し，Nativeモードで実行するこ（HLST）の好意により，伊にある MARCONI の MIC パーとで性能評価を行った．VTuneを用いてホットスポット解テション（インテル社製 Knight Landing（KNL）搭載）に

析 を 行 っ た．こ こ で MPI 並 列 数（-n）＝１２８，ノ ー ド 数

よ る 性 能 比 較 を 行 う こ と が で き た．KNL に 対 し て は（-N）＝６４とした．解析の結果，以下の３つのサブルーチン-xMIC-AVX512のオプションを追加した．KNLではElapse

においてコストが高いことが判明した．各モードのエネル

時間は 43.62 秒であり，KNC と比較すると約 5.37 倍程度性
８６

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

の形成をある程度再現することができた［１３３，
１３４］．

能が向上している．今回，以下の２点についてプログラム
をチューニングし，スレッド数をパラメータとした経過時

その本質は，材料中のヘリウム原子の拡散を KMC に担間計測を行った．（１）サブルーチンtmpusでLAPACKライ

当させたことで，現実の 0.5 秒に相当する長時間スケール

ブラリを call するループを，OpenMP により並列化した．

のシミュレーションを可能としたことにある．これによ（２）サブルーチン tmenr の FFTW をスレッド並列版に変更

り，プ ラ ズ マ 粒 子 の 入 射 フ ラ ッ ク ス を 現 実 同 様 のした．この結果，増強システムで FFTW のスレッド数を８1022 m−2s−1 まで下げることができたことにある．MD シにした場合，121.5 秒，Helios 本体システムで 20.15 秒まで

ミュレーションだけを用いた場合は，その再現時間の短さ高速化することができた．性能差は約６倍まで短縮するこから入射フラックスを数桁上げざるを得ず，粒子の入射とができた．今後の課題として，キャッシュミスヒットにと，材料内部拡散を経た再放出の競争が正しく再現できなより，性能飽和が判明したため，キャッシュチューニング

くなってしまう．ハイブリッド化によりこの点を改善できを行う必要がある．具体的には配列の引数の入れ替え等をたわけである．

行うことで，メモリのアクセスパターンを変更しメモリコ

この知見をもって，第４サイクルから第５サイクルでピーの速度改善を図る．また，NAG LAPACK もしくはは，より汎用性の高いシミュレーション手法へと拡張すべ

ScaLAPACK はスレッド並列に対応しているので，それにく，水素照射によるリテンションの算出を課題とした．こ置き換えることで，さらなる速度改善が見込めるかもしれ

こでは，BCA-KMC ハイブリッド手法を開発した［１３５］．

ない．個人的な印象としてMARCONIのMICパーティショ

材料中の原子拡散をKMCで扱う点は以前の手法を継承し，ンは Helios 本体システムと同程度（約８０％程度）の性能を今度は照射部分を BCA で解き（図４９），より高エネルギー示しているが，性能を引き出すためにはデータ配置等を工かつより長時間スケールへと対応させた．入射フラックス

夫しメモリーコピーを最適化する必要があり，一般ユーを 1020 m−2s−1 に設定した場合，最長で３０秒に相当するシザーが利用するには少し閾値が高い計算機と言えるだろミュレーションが２４時間程度で完了する．
また，KMC 部分もMD-MCに比べるとより詳細なモデル

う．

へと進化し，現実の格子間の原子移動経路に沿った拡散を

６．炉材料扱えるようになった．特に，粒界やアモルファスの様な対６．
１ プラズマ−物質相互作用の分子動力学研究

称性の低い構造中での原子移動経路と，そこでの移動障壁

Molecular Dynamics Simulation for Plasma-Material Inter-

エネルギーを自動的に探索する手法を開発した［１３６］．本action（MDPMI）

手法では図５０のように材料を構築する任意の原子配置に対プロジェクト代表：伊藤篤史［核融合科学研究所・総合研

して，局所的に切り出した領域内で構造緩和の MD を行う

究大学院大学］

ことで拡散原子の安定点を見つける．空間的に満遍なく切

日本枠

り出して計算することで，全ての安定点を洗い出す．続い第１−第５サイクルて，安定点の間が移動可能な経路となりえるかどうかを調
核融合炉の内壁，特にダイバータ部分におけるプラズ

べるため，MD による Nudged Elastic Band（NEB）法を行

マ−物質相互作用では，材料表面のミクロな損耗や破壊がう．これにより移動障壁エネルギーが算出される．各局所重要である．損耗によってプラズマ中へ放出される不純物

領 域 の 計 算 は 独 立 に 実 行 可 能 で あ り，Multi-Program

も分子レベルのミクロな物質であるし，材料表面のミクロ

Multi-Data 形式でロスの少ない並列化が行われ，高速な移な凹凸による表面積の増加は水素同位体のリテンションに

動経路探索が可能となる．

大きな変化を与える可能性がある．よって実験研究においても照射試験をした試料表面の電子顕微鏡画像などによる解析が必要不可欠である．それに呼応して，理論的な研究においても，原子レベルの挙動を解くシミュレーション研
究が必要である．
本プロジェクトでは，主にタングステン材料表面のプラ
ズマ−物質相互作用現象を解明すべく，分子動力学（MD）
を中心に，二体衝突近似（BCA），密度汎関数理論（DFT），
動的モンテカルロ法（KMC）を使い分けたマルチスケール解析を実施した．この中で，DFT および KMC に必要な知見は後に述べる FPSPFM プロジェクトと密に連携することで大いに得られた．
第１サイクルから第３サイクルの間は，ヘリウム照射で誘起されるタングステン表面の繊維状ナノ構造を対象として，マルチスケール解析に必要な各種のコード開発から適

図４９ BCA-KMC ハイブリッド法，MD を利用した KMC 経路の自動探索，ならびに，MD 用ポテンシャルモデル開発の関連
図．モデルにおいても階層構造が成り立っている．

用までを幅広く行った．特に，MD と KMC を連結させた
MD-MC ハイブリッド手法の開発により，繊維状ナノ構造
８７

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

空孔に戻るような過程を経ることがわかった．従って，水素が多いほど移動エネルギーも大きくなる．空孔と水素同位体の結合エネルギーを計算した（図５２）．タングステン
空孔には水素は最大で１２個まで捕獲される．さらに水素同位体の零点振動エネルギーの違いにより軽い水素同位体ほど空孔との結合エネルギーは大きい．そのため，残留量が問題となっている放射性の三重水素が軽い水素同位体の影
響で量が少なくなることが期待される．今後の課題としてはこのような複数の水素同位体が共存する系の研究が必要になると考えられる．また，図５２が示すように水素の影響
で空孔（空孔水素複合体）の生成エネルギーは著しく小さくなる．空孔内に３個以上の水素が捕獲されれば空孔の生成エネルギーはマイナス，すなわち自発的な空孔の生成が起こる［１３８］．このような状態での水素の捕獲量の研究も行った．

図５０ MD を利用した KMC 経路の自動探索の概念図．局所安定サ
イトおよび移動経路は区切られた局所領域内での MD 計算
によって求められる．各 MD 計算は MPMD 並列化によって，個々の CPU プロセス上で非同期に実行される．

これらに加え，本プロジェクトでは MD に必要な原子間相互作用ポテンシャルモデルの開発も行ってきた［１３７］．
ま た 現 在 で は，MD-MC と BCA-KMC を 合 わ せ た，
BCA-MD-KMC 三 連 ハ イ ブ リ ッ ド コ ー ド を 開 発 中 で ある．このような一連の取り組みにより，プラズマ照射過程
から材料内部の現象までのマルチスケール現象を扱える
コード群が構築されてきている（図４９）．今後は，より汎用性を高め，ダイバータを始めとした将来の炉材料開発に役立つ実用性のあるコードへと成長させたい．
６．
２ 第一原理計算によるタングステン中の空孔拡散の研究
Study of vacancy diffusion in tungsten by first principle cal図５
１ タングステン空孔中の水素（H）の個数と移動エネルギーの関係．

culation（SVDT）
プロジェクト代表：大澤一人［九州大学応用力学研究所］
日本枠

第４−第５サイクル

タングステンは融点が高く水素溶解度が低いために，核融合炉材料の中でも最も激しい照射を受けるダイバータの表面を保護するために使われる予定である．タングステン
結晶に水素が溶解することはないが，照射で生成する格子欠陥に捕獲される可能性はある．そこでタングステン材料の照射損傷や水素吸蔵に関する研究が行われている．本研
究では，第一原理計算を使ってこれらの現象の研究を行った．タイトルには空孔拡散とあるが空孔と水素同位体の多岐にわたる相互作用のシミュレーションを行った．
タングステン空孔に水素が !個捕獲された空孔水素複合体（VH!）の移動エネルギーを計算した（図５１）．空の空孔の移動エネルギーは 1.8 eV であるが，捕獲される水素が増えるほど移動が難しくなる．移動の特徴は空孔と水素が常に一緒に移動するのではなく，空孔に捕獲されていた水素

図５
２ タングステン空孔中の水素同位体の個数とその空孔との結合エネルギー．

が一旦全部外に出て空の空孔が移動し，最後に水素が元の
８８Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

６．
３ 第一原理計算によるウラン合金中水素の挙動解析
Analysis of Hydrogen Behavior in Uranium alloys by firstprinciples calculation（AHBU）
プロジェクト代表：鬼塚貴志［福井大学附属国際原子力工
学研究所］
日本枠

第２−第５サイクル

DT 核融合炉開発において，余剰のトリチウムを安全に

図５
４ （a）
UNiAlHx 合金および（b）
UPtAlHx 合金の格子定数の水
素 濃 度（H/U）依 存 性．計 算 値［１
４
１，
１
４
２］
を白抜き記号で，実験値［１３９，１４０］を黒塗り記号で示す．H/U ＞ 0.5の範囲の UPtAlHx 合金は実験的に得られていない．

蓄積可能で良好な吸蔵特性をもつトリチウム貯蔵材料の開発は，安全性と経済性向上の観点において重要課題である．現在，有力な候補材として UNiAl 合金に代表される
UTX 系（U：ウラン，T：遷移金属，X：半金属）ウラン化した網羅的な計算評価による，水素吸蔵特性の優れた合金合物が挙げられている［１３９］．しかしながら，UTX 系化合の探索が必要である．現在，Selene の計算機システムを利物であってもUPtAl合金のように貯蔵能の低い合金も存在

用して引き続きこれらの課題に取り組んでいる．

し，経験的手法によって合金の水素吸蔵特性を予測するこ６．
４ プラズマ対向材料および構造材料の照射・熱時効下とは困難である．それゆえ，吸蔵特性のより優れた合金をでの微細構造発達

探索するには，計算科学的手法が不可欠である．本プロ

Microstructure evolution of plasma-facing and structural

ジェクトでは，第一原理計算を用いて UTX 系ウラン合金の水素同位体吸蔵特性の予測法を確立することを目的とし

materials under irradiation and thermal aging（MEPFSM）

た．

プロジェクト代表：鈴土知明［日本原子力研究開発機構
システム計算センター］

計算パッケージとして密度汎関数法に基づく第一原理計

日本枠算コードである VASP を用いた．まず実験値が比較的充足第４−第５サイクル

している UNiAl 合金と UPtAl 合金を対象とし，系の全エネ
ルギーの計算結果から水素化物の構造安定性を検討した．

本プロジェクトはプラズマ対向材料で起こる素過程の中

また，水素吸収にともなう，格子定数および最安定な原子で材料劣化に影響を与えるものを切りだして，精度の高い

位置等のミクロ構造の変化に関する計算を行い，それぞれモデリングを行うことを目的とした．以下では特に，タン実験結果と比較・検討した［１４０，
１４１］．

グステン（W）材中で生成される核変換元素の影響につい

図５３に UNiAlHx 合金の安定構造の水素濃度（H/U）によての研究を報告する．W に中性子照射されると核変換にる変化を示す［１４３］．UTX 系合金は一般に六方晶に属する
Fe2P 型構造をとるが，UNiAlHx 合金は H/U=0.67 を越える

よってレニウム（Re）やオスミウム（Os）が生成されるがそれらの元素が BCC 相の W 中で "相や !相として析出し，と c 軸上の Ni の安定位置が顕著に変化するという実験結果

材料の延性低下を引き起こすことが知られている．またこがある［１３９，１４０］．水素化物の構造安定性および，最安定なれとは別に，W に Re や Os を加えると照射スウェリングが原子位置等のミクロ構造に関する計算結果は，このような抑制されることが知られている．このような現象の発生メ

実験知見を良く再現する結果を得た．一方 UPtAl 合金の水カニズムはまだ知られておらず，今後 W 材の適用を考える

素化物はエネルギー的に不安定で，水素吸蔵能力が小さい

上で，解決すべき問題である．このような背景から第一原ことが計算からも示された［１４２］．図５４に UNiAlHx 合金

理計算をベースに，W 結晶中で Re や Os がどのように移動および UPtAlHx 合金の格子定数の水素濃度による変化をするかを正確にモデル化することを試みた．

示す．UNiAlHx 合金と UPtAlHx 合金の変化はそれぞれ異Re や Os は従来から空孔と位置交換することによって比なる傾向を示すが，いずれの場合においても実験と概ね一

較的ゆっくり移動することが知られていたが，格子間原子

致する傾向が得られた［１４０，
１４１］．以上の成果から，UTXとしての移動メカニズムは明らかにされていなかった．こ系合金の水素吸蔵特性を計算によって十分に予測可能であれらの格子間原子は W 原子と安定な混合ダンベルを形成ることが確認できた．

することが知られていたため，多くの研究者はその移動障

今後は実験的知見の少ない多様な UTX 合金系を対象と壁は高い，すなわち移動が困難と考えていた．しかしながら，Re に関して様々な移動様式を探索した結果，図５５に示したように第一近傍へのジャンプとその場での回転によって非常に容易に拡散することが明らかになり
［１４３］
，Os に
ついても類似の結果を得た．これにより，Re や Os は照射下では容易に移動し，析出もしやすいことがわかった．その上，純 W 中の格子間原子が１次元運動であるのに対して，これらの拡散次元は３次元になるため，周辺の空孔との再結合確率が増えることから，ボイドの形成すなわちスH/U ＝ 0.667，（c）
H/U ＝ 2.0
図５
３ UNiAlHx 合金の（a）
H/U ＝ 0，（b）の場合の安定構造の計算結果［１４３］．

ウェリングを抑制する効果も明らかになった［１４４］．
８９

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

る，という手順により計算を実行した．初期状態を設定する際，望まない準安定状態に収束することを避けるため，想定される平衡原子位置から適切な変位でランダムにずらす工夫をした．
タングステン空孔に捕捉されたヘリウム原子系に対する
図５
５ 〈１
１
１〉混合ダンベルの移動パターン；第一近傍間でのジャンプと回転によって，Re の格子間原子は容易に移動できる．

研究［１４５］を踏まえ，タングステン以外のいくつかの遷移

今後は，スウェリング抑制効果を別の析出を起こさない

は，核融合炉材料研究のための基礎データのみならず工学

金属中の単一空孔にヘリウムが捕捉された場合に対する計算を行い，結合エネルギーとして図５６を得た．このデータ

元素で代替できないか，そして W 材のもう一つの大きな問

応用を考える上でも有用となりうる．

題点である，脆性を合金元素添加によって改善できないか

空孔のない（BCC 構造をもつ）タングステン中にヘリウを計算科学的手法で探っていくことになる．

ムおよび水素が格子間原子として存在する場合についてシ
ミュレーションを行った．また，得られた安定状態に対し

６．
５ 第一原理計算によるプラズマ対向材料の研究て電子密度分布解析を行った．その結果，ヘリウムは（空First Principle Simulation for Plasma Facing Materials

孔がなくても）格子間原子としてクラスターを形成しうる（FPSPFM）

一方，水素は格子間原子としてはクラスターを形成し難い

プロジェクト代表：!山有道［核融合科学研究所・総合研

がヘリウムが介在すれば複数の水素が格子間原子として集究大学院大学］

合しうることが明らかとなった．これらの結果はタングス

日本枠

テン結晶格子間に捕捉されたヘリウムの周辺に電子密度の第２−第５サイクル

低い領域が形成されそれが新たな捕捉サイトとなることで核融合炉においてプラズマと対向するタイバータの材料説明しうることを示した．ここで得られた結果は，ヘリウ

としてタングステンは有力な候補であり，精力的に研究がムリッチなクラスターが水素の捕捉サイトとして機能する

進められている．ダイバータは，核融合の燃料である水素ことにより水素の拡散挙動に影響を及ぼしうること，水素および核融合反応生成物であるヘリウムに曝される．タンによるタングステン材のブリスタリングがヘリウム照射にグステンは融点が高く，スパッタリング耐性も高いが，ダより抑制されるという実験事実と関係しうること，を指摘

イバータのような過酷な環境下では強く損傷を受けることした［１４６］．

が懸念され，プラズマ−物質相互作用を理解することは重タングステン中の空孔に捕捉されたヘリウムの動的挙動要である．本プロジェクトは，タングステンとヘリウムをを調べるため，単空孔および二重空孔の移動障壁エネル

はじめとする希ガスおよび水素の共存系に対して密度汎関

ギーおよび単空孔に捕捉されたヘリウム１原子および２原

数法に基づく第一原理シミュレーションを実行し原子ス

子の移動障壁エネルギーを第一原理計算により評価した．

ケールでの特性およびポテンシャルエネルギー面を得るこ

その結果，単空孔は［１１１］方向へ移動しやすいこと，二重空と，得られたポテンシャルエネルギー面に基づいて各種状孔の移動は単空孔の移動が連続して起こるという見方で説

態間の遷移確率モデルを構成し動的モンテカルロ法による明しうるとの理解を得た［１４７］．これらの第一原理シミュ

シミュレーションが実行可能なコードを開発すること，をレーションによって得られた結合エネルギーおよび移動障目的として遂行した．最終的にめざすところはタングステ

壁エネルギーのデータを用いて，タングステンにヘリウム

ンダイバータでのプラズマ−物質相互作用を理解することを照射した際にナノ構造（ファズ構造）形成の前段階で観にあるが，これは核融合研の伊藤が代表としてIFERC国際測されるヘリウムバブルの生成過程に対する動的モンテカ

枠で実施されたプラズマ−物質相互作用のための分子動力学シミュレーション（MDPMI）プロジェクトと密に連携することにより取り組んだ．
シミュレーションは，東京大学の尾崎らが主として開発しているオープンソースコードOpenMXをMDPMIおよび本プロジェクト遂行のため伊藤がチューニングしたものを用いた．OpenMX は密度汎関数法に基づく第一原理シミュ
レーションコードであり，交換相関ポテンシャルとしてPerdew-Burke-Ernzerhof 汎 関 数 に よ る 密 度 勾 配 近 似（GGA-PBE）を，コアクーロンポテンシャルに対して擬原
子局在軌道およびノルム保存擬ポテンシャルを用いた．初期状態として計算対象とする系の原子配置を設定し，自己
無 撞 着（SCF）収 束 さ せ，さ ら に Broyden-FletcherGoldfarb-Shanno（BFGS）法により体積保存下で全原子位

図５６ 各種遷移金属の単空孔中に捕捉されたタングステン原子の結合エネルギー．

置を緩和させて構造最適化を行うことにより最終状態を得
９０

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

ルロシミュレーションを行い定性的な再現をみた（図５７）．

BA Project Planning and Coordination Group et al.

"
!"!
κ2" !
!
!
!
"
"
"
!!!
#"

本プロジェクト終盤には結晶粒界近傍に対する第一原理
計算に取り組んだが，結果をまとめるまでには至っておら（３）

ここで #は温度，!"" と !"!はそれぞれ熱フラックスのず，今後の課題として残っている．

自己相関関数と熱フラックス‐質量フラックスの相関関数
６．
６ マルチスケールモデリングによるトリチウム増殖材の現象論的係数であり，MD 計算により評価できる．２式料の性能評価Multi-scale modelling for long-term performance of tritiumにより計算される熱伝導率の比較を，２元系（Li2O，TiO2）で進めたところ，κ2 は実験結果と類似の温度依存性を示し

breeding materials（MMLP̲TBM）

た一方で，κ1 は高温において上昇するという不自然な挙動

プロジェクト代表：谷川 尚［量子科学技術研究開発機構］を示した（図５８）．２式の相違を解析的に分析した結果，ポ国際枠

テンシャルエネルギーの絶対値の物理的意味は乏しいにも関わらず，κ1 ではそのエネルギーの輸送が熱伝導率に寄与

第２−第５サイクル

Li2TiO3 などの固体増殖材料中では，高速中性子によるする形になっていることがわかった．そのため，κ1 は原子

弾き出しやトリチウム生成に伴う Li 燃焼により，核融合炉輸送に伴うポテンシャルエネルギー輸送が発生する状況，

運転時に大量の照射欠陥が形成される．そして照射欠陥

つまり拡散係数が大きくなる高温や欠陥を含む系においは，熱伝導率の低下やトリチウム放出速度の低下などの性て，不適切な結果を与える可能性がある．

能劣化を引き起こす．本プロジェクトでは，固体増殖材料

本プロジェクトで実施した他の計算において，Li2TiO3の照射損傷過程を理解し，照射欠陥の蓄積がどのように材の重要な照射損傷は，弾き出しにより形成されるアンチサ

料の熱機械的特性や化学的安定性に影響を与えるかを解明

イト型欠陥（TiLi）と，Li 燃焼により形成される Li2O 欠損することを目的として，量子力学計算，古典分子動力学であることが示唆されている．そこで，これらの欠陥が（MD）計算，速度論計算等を利用したマルチスケールモデ

Li2TiO3 の熱伝導率に与える影響を（３）を利用して評価しリングを展開した．計算は，量子力学計算コード VASP

た（図５９）．まず，いずれの欠陥も熱伝導率を低下させるこや，MD 計算コード LAMMPS 等を利用した並列計算を主とが確認できた．そして，アンチサイト型欠陥はその形成に実施した．以下では，熱伝導率に照射欠陥が与える影響の容易さにも関わらず，Li2O 欠損と同程度の大きさの影響について得られた知見を概説する．

を与えることがわかった．

まず，MD 計算による熱伝導率評価において利用されて既往の増殖ブランケット性能評価においては，Li2TiO3いる，グリーン・久保公式に基づく手法の比較検討を行っの物性値として未照射材の値が利用されてきた．本研究でた．その結果，異なる以下の２つの式が利用されているこ得られたような照射損傷の影響を含んだ物性値を利用するとがわかった［１４８，１４９］．

ことで，長期的な性能評価の信頼性向上が期待できる．しκ1" !""#"

かし，計算結果を実際の設計に利用するためには，さらなる計算精度の向上が求められる．特に，MD 計算で利用す（２）

るポテンシャルモデルの精度向上は重要である．本プロ
ジェクトにおいても，ポテンシャルモデルの精度向上のために，フーリエ級数展開を用いて２体力モデルを構築する手法を開発した［１５０］．現在，類似の方法を用いた多体力モデルの構築手法の開発を進めている．また，近年その利

図５
７ 動的モンテカルロシミュレーションによるヘリウムバブル
構造形成．

図５
８ 異なる計算式により評価した Li2Oの熱伝導率の温度依存性
の比較．

９１

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019

図６
０ Y2Ti2O7 および Fe 中の四面体（tetra.），八面体（octa.）
および空孔（sub.）
位置に占有した（a）
He および（b）
H の欠陥形成エネルギー．エネルギーが低いほど安定であることを示す．

と，He および H はそれぞれ Y2Ti2O7 および $-Fe 中に捕獲されることが明らかになった．この結果から，Y2Ti2O7 が分散された$-Fe では，He によるスウェリングが生じ難く，He
図５９ アンチサイト欠陥と Li2O 欠損が Li2TiO3 の熱伝導率に与える影響.

と H の同時蓄積による相乗効果も低いことが予測される．
６．
８ リチウムイオン伝導体を用いたリチウム同位体分離に関する研究開発

用が急拡大している機械学習を利用したモデルの開発も，より高精度の計算結果を得るには重要になると考えられ，

Research and Development of Lithium Isotope Separation

今後の研究課題である．

Technology by using a Lithium Ionic Superconductor（LISTLIS）

６．
７ ナノ微細粒子分散型フェライト鋼中のヘリウムおよプロジェクト代表：星野 毅［量子科学技術研究開発機構］

び水素の相対的安定性

日本枠

第４−第５サイクル

Relative Stability of Helium and Hydrogen in Nanostructured Ferritic Alloys（HEHFA）
プロジェクト代表：土屋
日本枠

核融合炉の燃料として必要なトリチウムは，核融合炉ブ

文［名城大学理工学部］

ランケット内に装荷されるトリチウム増殖材料中のリチウと中性子との核反応により生産する．しかしながム 6（6Li）

第４サイクルら，天然のリチウムには6Li が最大で約 7.8％ しか存在せず，核融合炉の定常運転には，6Li の存在比を約９０％に濃縮し１．緒言
フェライト鋼（$-Fe）にナノサイズの微細な酸化物粒子

たリチウムが必要となる．研究代表者は，これまでに，セ（Y2Ti2O7）を分散させた補強材料が核融合炉の構造材の候

ラミックス製のイオン伝導体をリチウム分離膜とした，海補の一つとして挙げられている．転位の運動の抑制が材料

水からのリチウム回収技術を発案した［１５１］．イオン伝導中に 添 加 さ れ た ナ ノ 微 細 Y2Ti2O7 粒 子 に よ っ て 促 進 さ

体中をリチウムのみが透過することに着目した技術であれ，材料強度が高まることが期待される．さらに，分散粒

る．6Li は7Li より移動度が高いことから，原理的には本技

子がプラズマ粒子として照射される He，H 等を捕捉し，ボ術を発展させることで，6Li と7Li の分離が可能となる．本イドスウェリング（膨張）や粒界脆化などを抑えることが研究では，イオン伝導体中の6Li と7Li の移動メカニズムを期待される．本研究では，第一原理計算を用いて，Y2Ti2O7

計算科学により探索するとともに，6Li と7Li の分離試験を

および $-Fe 中の四面体および八面体位置に捕獲された He

行った．

および H の欠陥形成エネルギーを求め，Y2Ti2O7 が分散さ

Li 拡散経路を nudged elastic band 法により探索し，フォれたフェライト鋼中の He および H の最も安定な占有位置

ノン計算を行うことにより，単位時間当たりに Li イオンがを決定することを目的とした．

サイト間をジャンプする頻度 #を計算した．また，キネ

２．実験

ティクス・モンテカルロ法を実行することにより6Li，7Li第一原理計算プログラム（VASP）を用いて，He あるいの自己拡散定数を評価した．計算には第一原理計算ソフト

は H の１原子をパイロクロア型 Y2Ti2O7 中の四面体およびvasp を使用した．Li 同位体拡散に関して得られた知見は以八面体型の格子間位置に占有させた時の形成エネルギーを下の二つである．

計算した．さらに，体心立方晶の $-Fe 中の四面体，八面体①：Li イオンが狭い領域から広い領域に移動する結晶構

格子間および空孔位置に捕捉された He および H の欠陥形造の場合に量子効果は同位体拡散比を上昇させる．

成エネルギーについても計算した．

②：空孔を介した6Li，7Li 間の相関効果により同位体拡３．結果・考察

散比が減少する．

Y2Ti2O7 中に捕捉されたHeおよびHの形成エネルギーを

イオン伝導体LiTi2O4 及びLiC6 における6Li，7Liジャンプ

図６０に示す．He および H はそれぞれ Y2Ti2O7 中の八面体および四面体型の格子間位置を占有することが判明した．さ

頻度比 #""##を評価した（図６１）．①を満たす構造を持つ

らに，$-Fe 中の He および H の形成エネルギーと比較する

LiTi2O4 では，室温においてジャンプ頻度比は #""
!
##!!
!
９２

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

６．
９ 材料照射相関の解明に向けた計算機シミュレーショ

と な り，古 典 統 計 力 学［１５２］から予測される 値!
$$$#"!&&を上回る結果となった．①と逆の構造を%$' %"持つ LiC6 では，#$"!
#%#"!$となり，同位体効果が小さくComputational Study for Materials Irradiation Correlation

なった．これらの結果は，結晶構造の違いにより，拡散に（COSMIC）

おけるゼロ点振動の寄与が異なる事に起因する．また，キプロジェクト代表：渡辺淑之［量子科学技術研究開発機構］

ネティクス・モンテカルロ法により求めた LiTi2O4 におけ

国際枠

ン研究

第４−第５サイクル

る同位体拡散定数比を求めた結果，!$"!%"1.05〜1.06 程
度の値となり，単一イオンの動きやすさを表すジャンプ頻

本研究では，高エネルギー粒子線照射に晒された材料の度 #の計算から予見される比（#$"!
#%#"""）を下回る値とミクロ構造変化の数理モデルの構築に必要となる，はじきなった．我々は相関因子等の計算から，空孔を介して生じ出し欠陥生成量や欠陥エネルギーなどの原子レベルの情報

た6Li，7Li 間の相関効果が同位体拡散比を小さくしているの獲得を目的としている．主な対象材料は低放射化フェラ

ことを明らかにした．以上の結果から，①を満たす構造で

イト鋼（RAFM），炭化ケイ素（SiC），タングステンでああり，且つ空孔機構が弱い物質において Li 同位体効果が大り，同材料中の点欠陥（原子空孔，自己格子間原子，ヘリきく現れることが期待される．これらの条件を満たす物質

ウム・水素原子など）やそれらの複合欠陥の挙動についてと し て，イ オ ン 伝 導 体 Li3xLa2/3−xTiO（LLTO）におい3の解析を実施した．ここで，ヘリウムや水素同位体は，材て，リチウム濃度（x）が小さくした組成が挙げられる．

料内の格子欠陥（転位，欠陥集合体，析出物，粒界，酸化67次に，実験により Li， Li の同位体分離効果を検証した．

膜など）と強く相互作用することで材料の特性・形状変化

0.1 M の LiOH 水溶液 10 L を原液として，海水リチウム回収を促進させることが懸念されているが，とりわけ水素同位技術と同様の手法にてイオン伝導体 LLTO に電圧を 5 V

体と格子欠陥の相互作用については未だ十分な理解が得ら

印加し，回収液中の6Li の割合を分析する透析試験を行っれていない．本稿では，RAFM 中の主要酸化膜である酸化た（図６２（a））．７日間の電気透析の結果，試験期間中，6Li

鉄（FeO）を対象とし，第一原理計算により同結晶中の水の割合は天然比の 7.8％ を大きく上回る約 8.0％ を維持して素同位体の存在状態を電子論的に評価した内容について報おり（図６２（b）），同位体分離係数は最大で 1.04 と高い値を告する．

得ることに成功した．

密度汎関数法に基づいた第一原理計算コードである SI-

今後は同位体分離効果に優れたイオン伝導体の材料探索

ESTA［１５３］を使用した．SIESTA では各原子に帰属する実を計算科学にて行うとともに，実験結果の計算へのフィー

在波動関数を基底関数として用いており，格子欠陥周辺なドバックを試みる予定である．

どの局所的な電子状態を知るための効果的手法である．本
計算では原子数６４個の FeO スーパーセルに水素原子を導入し，それぞれ０K・体積一定の下で系を構造緩和させ，水素原子の安定サイト（固溶エネルギーが低い配置）を調
べた．第一原理電子状態計算における構造最適化計算は原子の初期配置の影響を強く受け，初期配置近傍の極小値を解として得る傾向にある．そこで水素原子を導入する際には，FeO における結晶構造の対称性を考慮し，１２パターンの初期配置を設定し，それぞれについての計算を行った．
ま た，尚，水 素（H）原 子 の 固 溶 エ ネ ル ギ ー は（FeO-H）（FeO）"# ""('""('(H2 に て 算 出 し た．
#
!"('
!"
(
(
（FeO-H）
は H 原子を含む FeO 結晶のトータルエネル"('(
ギー，"('（FeO）はFeO 完全結晶のトータルエネルギー，("('(H２は真空中の水素分子の結合エネルギーである．

図６
１ （a）
LiTi2O4, LiC6 における6Li，7Li のジャンプ頻度比 #6/#7の計算結果．(b)Li イオンの拡散経路.

図６３は，計算より得られた FeO 中の H 原子の安定サイトとその時の固溶エネルギーを示している．最も安定なサイト（Site-1）は６個の Fe を頂点とした８面体の中心付近に酸素と近接して存在していた．このときの H-O 間の距離は約 1.1 Åである．また，２番目に安定なサイト（Site-2）は２個の Fe と２個の O を頂点とした菱形の中心付近に存在する．このときの H-O 間距離も約 1.1 Åであった．また，得ら れ た H 原 子 の 固 溶 エ ネ ル ギ ー は そ れ ぞ れ 0.074 eV，0.132 eV で あ り，純 鉄 中 で の H 原 子 の 形 成 エ ネ ル ギ ー
0.25 eV［１５４］よりも低い．これは，Fe ベースの母相を持つ
RAFM においても酸化膜が水素の拡散障壁として働き，水

図６
２ （a）
実験概念図（6Li が早く移動）（b）
回収液中の6Li の割合．

素同位体の強いトラップサイトになっている可能性を示唆
９３

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019ては，米ロスアラモス研で開発されている MCNP5/6 を主に用いた．並列化には MPI および MPI/Open MP ハイブ
リッドを用いており，２５６コアを標準として，最大４０９６コアまでを用いた．また，重陽子入射による核反応（中性子，ガンマ線生成）のシミュレーションのため，PHITS コードに実装された物理モデルの検証も実施した．様々な分散低減法を組み合わせて，統計的に有意な結果を得ることに注意を払う必要があるが，必要な計算精度が達成できることを確認した．
第４サイクルでは，5 MeV までの重陽子パルスビーム（デューティー比 0.1％）を用いた，いわゆる「フェーズ B」

図６
３ FeO における水素原子の安定サイト．

試験のための解析を行った．得られた空間線量分布の例を
図６４に示す．この結果を基に，運転時の線量が法令基準をするものである．続いて，FeO 中における水素の拡散挙動満たすことを確認し，原子力規制庁に許認可申請を行っを調べるため，水素が格子間位置をジャンプするための拡た．加速器の運転に必要な許可は２０１６年３月３１日で取得す

散障壁を調べた．ここでは Site‐１に位置する水素原子が隣ることができた．ま た，第５サイクルでは，「フ ェーズ

接する等価な Site‐１に移動する際の活性化障壁を評価した

B」試験時の機器（ビームダンプ）の放射化による運転停止

結果，活性化エネルギーは約 0.26 eV であった．これは，同後の線量率の解析を実施し，安全に保守作業を行える見通

酸化膜中において水素は室温のような比較的低温においてしを得た．停止後線量の解析結果を図６５に示す．このほか，も容易に拡散できることを示唆している．これらの知見

9 MeV ま で の 重 陽 子 直 流 ビ ー ム を 用 い た，い わ ゆ る

は，RAFM における水素同位体の拡散・放出メカニズムの「フェーズ C」試験のための解析を行い，必要な追加遮蔽等

究明やトリチウムインベントリー評価モデルの構築に有用の条件を明らかにした．

である．今後は，酸化膜だけでなく粒界やラス境界などの今後は，上記「フェーズ C」の許認可申請に向けて，追加複数の相が共存する場合における水素の存在状態を解析

遮蔽等を含む申請条件に基づく詳細な遮蔽計算，放射化計

し，RAFM 中の水素同位体欠陥エネルギーの場所依存性を算を行う必要がある．本計算においては，およそ１４桁の中

系統的に明らかにする必要がある．

性子束の減衰を精度よく計算する必要があり，チャレンジ
ングな取り組みとなる．また，加速器室内の空気の放射化，

７．炉工学

冷却水の放射化による即発ガンマ線等に注意を払う必要が

７．
１ 並 列 モ ン テ カ ル ロ 計 算 に よ る IFMIF 原 型 加 速 器

あり，それぞれ計算方法の工夫が必要である．

（LIPAc）の中性子工学解析
LIPAc Neutronic Analyses by means of Monte Carlo Parallel Compution（MCLIPAC）
プロジェクト代表：近藤恵太郎［量子科学技術研究開発機
構］
国際枠

第４−第５サイクル

現在，量子科学技術研究開発機構

六ヶ所核融合研究所

では，日欧共同で IFMIF 原型加速器（LIPAc）の研究開発が進められている．MCLIPAC プロジェクトは，原型加速器の中性子工学解析のため，大規模並列モンテカルロ計算
による放射線輸送計算を行うものである．このプロジェク
トの目的は，原型加速器の核解析のための計算ツールと方法論を確立し，加速器施設の運転許認可を取得するための，総合的な放射線安全解析を実施することである．
IFMIF 原型加速器では，ビーム加速に用いる大電力の高周波入射のための導波管と，加速空洞を冷却するための冷却水配管のため，建屋の貫通部が多く，遮蔽計算では放射線ストリーミング効果を厳密に考慮する必要がある．このため，解析計算ではなく，モンテカルロ法による遮蔽計算が必須となるが，およそ１１桁の中性子束の減衰を計算する

図６
４ 得られた IFMIF/EVEDA 原型加速器運転中の空間線量分布の例（地下ピット，ダクト貫通部を含む建屋の垂直断面の線量分布）．

必要があり，かつ建屋が比較的大きい（約 40 m×20 m）ため，大規模並列計算が不可欠である．輸送計算コードとし
９４

Project Review

Simulation Research Projects Using Helios at IFERC-CSC

BA Project Planning and Coordination Group et al.

つの配列の軸反転を行うアルゴリズムを用いている．チャ
ネル乱流の場合は，周期境界条件の方向にスペクトル法を用いているが同じコードで２つの座標系を扱うことができる．本研究の乱流直接数値シミュレーションのプログラミ
ングコードは MPI + OpenMP による並列化をしている．
円管内乱流において，格子数は，流れ方向，径方向，周方向にそれぞれ２
０４８×５１２×７６８を用いた．壁面摩擦速度と円管半径で無次元化されたレイノルズ数は，１０５０である［１５５］．十分発達した乱流場が得られ，速度場のみならずスカラー場の統計量を得ることができた［１５５］．図６６にハ
イブリッド並列に成功した計算例を示す．１
２８にランクの数を固定した時にホスト CPU，MIC，CPU+MIC を様々な組み合わせで稼働させた時の速度比較を示す．円管内乱流において流れ方向，径方向，円周方向にそれぞれ２０４８×
５１２×７６８の格子 数 で 最 大６４ノ ー ド ま で の 計 算 を 行 っ た

図６５ 得られた IFMIF/EVEDA 原型加速器運転後の停止後線量率の解析結果の例（低出力ビームダンプ周辺での停止後線量
分布）．

［１５６］．MIC は CPU+ MIC の組み合わせであるヘテロジニ
アス計算において MPI でランクを作ることができるので，演算コアとして使用することができる．つまり，コア

７．
２ 剪断乱流の大規模スケール MHD シミュレーション

内でランクを振ることができるため，全体の並列化効率を

Large-scale MHD simulation of turbulent shear flow向上させることができる．ホスト CPU，ネイティブモー（LSMHDTSF）

ド，シンメトリックモードの速度を比較すると，ホスト

プロジェクト代表：佐竹信一［東京理科大学基礎工学部］

CPU が最速となった．そして次は，シンメトリックモー

国際枠

ド，ネイティブモードであった．ノードの数が１６，３２の時

第４−第５サイクル

シンメトリックモードとホスト CPU 及び MIC が同じ数の
このプロジェクトの目的は，核融合炉において重要な第ランクを使用する場合はシンメトリックモードのほうがホ

一壁冷却特性を検討するため，冷却材として Flibe‐から Li

スト CPU より速くなる．ヘテロジニアススーパーコンを想定した MHD 乱流における熱輸送特性に関する大規模ピュータでシンメトリックモードが使える環境であるな

シミュレーションを実施することである．特にターゲット

ら，少ないノード数で高性能な計算を行うことができるととするレイノルズ数は，実験的研究で行われているレイノいうことが明らかになった．

ルズ数１０万以上の規模に相当するシミュレーションを行う

本研究では磁場下の乱流現象解明ということが一番の目ことを予定している．実験データでいわれている主流れの的であるが，次世代の計算機でも搭載されると考えられる

乱れのピーク値のレイノルズ数依存性や円管中心部の大規

アクセレータの搭載もしくはスレッド計算可能な並列計算

模構造の役割などをシミュレーション結果から考察する．

機のための並列化手法の開発という目的もある．将来ヘテ

さらに磁場を印加し，得られるDNSデータから詳細な乱流

ロジニアス並列を想定しているため本計算コードのチュー

エネルギー収支式の各項の詳細を得ることで MHD 乱流モニング目標は，MPI + OpenMP の組み合わせの並列化におデル構築するための基礎知見を得ることである．次期計算

いて最適なメモリ分散の方法を開発することである．また機で標準搭載されるアクセレータ計算機の最適並列化手法

この方法を用いて，より大型で高レイノルズ数の円管内乱の開発という意味で大きい．MIC で大規模並列を想定して流で磁場下における乱流統計量を得ることを今後予定してい る た め 本 計 算 コ ー ド の チ ュ ー ニ ン グ 目 標 は，MPI+
OpenMP のいわゆるシンメトリック並列の最適なメモリ分散の方法を開発することである．またこの方法を用い
２０４８×５１２×７６８の系での印加磁場下で乱流統計量を得ることを主たる目的としている．
計算コードを MPI から MPI+OpenMP のハイブリッド
コードに変更し高速化を図った．基礎方程式は円筒座標において３次元非圧縮性 Navier-Stokes 方程式，連続の式を用いており，圧力ポアソン方程式の周期方向にはフーリエ変換を径方向の離散化には２次精度中心差分を適用し３重対角行列解法を用いている．時間離散化には非線形項に３次精度 Runge-Kutta 法，粘性項は Crank-Nicolson 法により時間積分を行った．さらに並列化についてのアルゴリズム
は前述のフーリエ変換を行うために別々の並列軸をもつ２図６６ Helios における６
４ノードまでのヘテロジーニアス計算．

９５

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019いる．
７．
３ 核融合用負イオン NBI システムの開発に向けた負イ
オン源プラズマの数値シミュレータ
Numerical simulator project for negative ion source plasmas
for developing fusion N-NBI heating systems（N-IOS）
プロジェクト代表：宮本賢治［鳴門教育大学］
日本枠

第５サイクル

核融合プラズマの追加熱・電流駆動装置の１つに，負イオンを用いた中性粒子入射装置（N-NBI）がある．負イオ
ンビームには発散角の小さいビームコア以外に，ビームハ
図６７ 数値計算モデル．

ロと呼ばれる発散角の大きい成分が観測されている．NNBI で大電流負イオンビームを高エネルギーに加速する
表２

際，ビームハロは加速部での電極熱負荷や絶縁破壊の原因の１つと考えられるため，N-NBI を実現する上で，ビーム加速部の電極熱負荷（荷電粒子の加速に要したパワーに対する割合）．

ハロの低減は必須の課題である．本研究の目的は，３次元

加速電流

実験結果
（0 Paに外挿した値）

数値計算結果の Particle In Cell（PIC）コードを用いて負イオンビームの

A1G
A2G
GRG

2.4％
3.0％
5.0％

0.9％
2.4％
3.7％

引出・加速のシミュレーションを行い，ビームハロ発生の物理機構を明らかにすることである．
負イオン源の負イオン生成部から加速部までを，総合的にモデル化した（図６７）．負イオンビームの引出・加速方

対して，負イオン引出し面の端部から引き出された負イオ

向を !軸とし，またフィルター磁場，電子抑制磁場に平行

ンはビームハロになる．負イオン引出し面の端部の電界レな方向をそれぞれ "軸，#軸とした．計算モデル領域は実サ

ンズ効果は強いため，引出部で負イオンビームが過収束しイズで，309 mm×19 mm×19 mm であり，これを２４７２×てビームハロになる．
引出部からの２次電子が考慮されていないため，A1G

１５２×１５２個のメッシュに分割した．
電子，H＋，H−の粒子種の運動方程式を解き，各メッの電極熱負荷の計算結果は実験結果に比べて大幅に低いとシュ領域での電位や粒子密度を求めた．各粒子数は 2.96×

考えられる．今後，引出部からの２次電子をモデルに考慮

107 個（電子），3.11×107 個（H＋），0.15×107 個（体積生成する予定である．

されたH−）である．負イオン源における物理パラメータの測定値を基に，電子密度は 1.0×1018 m-3 とし，また電子温

参考文献

度と正イオン温度はそれぞれ 1 eV，0.25 eV とした．プラズ［１］中島徳嘉 他：プラズマ・核融合学会誌 91, 711 (2015).
［２］中島徳嘉 他：プラズマ・核融合学会誌 92, 612 (2016).
［３］日本原子力研究開発機構 IFERC-CSC 利用 委 員 会 他
編：プラズマ・核融合学会誌 92, 157 (2016).
［４］文 部 科 学 省：平 成３
０年 版 科 学 技 術 白 書 http://www.
mext. go. jp / b ̲ menu / hakusho / html / hpaa 201801 /1398098.htm［５］D. Borba et al., Proc. 26th IAEA Fusion Energy Conf. (Kyoto,
Japan, 2016) OV/P-9.
［６］N. Nakajima et al., IFERC-CSC Final Report (2017).
［７］A. Ishizawa et al., Nucl. Fusion 53, 053007 (2013).
［８］A. Ishizawa et al., J. Plasma Phys. 81, 435810203 (2015).
［９］A. Ishizawa et al., Phys. Plasmas 21, 055905 (2014).
［１０］A. Ishizawa et al., Nucl. Fusion 55, 043024 (2015).
［１１］A. Ishizawa et al., Nucl. Fusion 57, 066010 (2017).
［１２］Y. Asahi et al., Phys. Plasmas 24, 102515 (2017).
［１３］Y. Idomura, J. Comput, Phys. 313, 511 (2016).
［１４］Y. Idomura et al., Proc. 26th IAEA Fusion Energy Conf.
(Kyoto, Japan 2016) TH/P3-1.
［１５］Y. Idomura, Plasma Fusion Res. 11, 2403006 (2016).
［１６］Y. Idomura, Phys. Plasmas 24, 080701 (2017).
［１７］K. Imadera et al., Proc. 26th IAEA Fusion Energy Conf.
(Kyoto, Japan 2016) TH/P3-3.

マ電極（PG）表面から放出される表面生成負イオンの初速度分布は，1 eV の初期温度の半マクスウェル分布を仮定した．
PG 表面から放出された表面生成負イオンは，クーロン
衝突や荷電交換によって，エネルギー緩和される．これをモデル化するために，クーロン衝突については Binary Collision Model［１５７］を，また荷電交換については Null CollisionModel［１５８］をそれぞれ用いた．
まず，計算モデルの妥当性を評価するために，荷電粒子の加速に要したパワーに対する電極熱負荷の割合について，実験結果と計算結果を比較した（表２）．加速部の電極
熱負荷の要因として，１）残留水素ガス分子とH−との衝突によって生じた剥離電子，２）ビームハロ，３）引出部からの２次電子が挙げられる．0 Pa まで外挿した際の電極熱負荷は，ビームハロか，引出部からの２次電子による．表２では 0 Pa まで外挿した際の値を示した．両者がほぼ一致することから，計算モデルは妥当であると言える．
次に，ビームハロ発生の物理機構は，以下の通りである
ことを明らかにした．すなわち，負イオン引出し面の中央部付近から引き出された負イオンはビームコアになるのに
９６

Project Review

Simulation Research Projects Using Helios at IFERC-CSC［１８］Y. Kishimoto et al., Proc. 26th IAEA Fusion Energy Conf.
(Kyoto, Japan 2016) TH/P3-2.
［１９］K. Obrejan et al., Comput. Phys. Commun. 216, 8 (2017).
［２０］N. Kasuya et al., Plasma Sci. Tech. 13, 326 (2011).
［２１］N. Kasuya et al., Plasma Fusion Res. 8, 2403070 (2013).
［２２］N. Kasuya et al., Phys. Plasmas 21, 110701 (2014).
［２３］N. Kasuya et al., Plasma Fusion Res. 12, 1303005 (2017).
［２４］S. Satake et al., Plasma Phys. Control. Fusion 53, 054018(2011).
［２５］M. Honda et al., Nucl. Fusion 54, 114005 (2014).
［２６］M. Honda et al., Nucl. Fusion 55, 073033 (2015).
［２７］S. Matsuoka et al., Phys. Plasmas 22, 072511 (2015).
［２８］B. Huang et al., Phys. Plasmas 24, 022503 (2017).
［２９］W.I. van Rij and S. P. Hirshman, Phys. Fluids B 1, 563

BA Project Planning and Coordination Group et al.

sity, (Japan, 2016).
［６７］T.-H. Watanabe et al., Phys. Plasmas 22, 022507 (2015).
［６８］Y. Todo et al., Nucl. Fusion 54, 104012 (2014).
［６９］Y. Todo, New J. Phys. 18, 115005 (2016).
［７０］Y. Todo et al., Phys. Plasmas 24, 081203 (2017).
［７１］H. Wang et al., Phys. Rev. Lett. 120, 175001 (2018).
［７２］Y. Todo et al., Nucl. Fusion 56, 112008 (2016).
［７３］http://iopscience.iop.org/journal/0029-5515/labtalk/article/65845［７４］ http : / / iopscience. iop. org / journal / 0029-5515 / page /Highlights-of-2016［７５］H. Nuga et al., Nucl. Fusion 57, 086011 (2017).
［７６］A. Bierwage and K. Shinohara, Phys. Plasmas 23, 042511(2016).
［７７］A. Bierwage and K. Shinohara, Phys. Plasmas 23, 042512(2016).
［７８］A. Bierwage et al., Nucl. Fusion 57, 016036 (2017).
［７９］A. Bierwage et al., Proc. 26th IAEA Fusion Energy Conf.
(Kyoto, Japan, 2016) TH/4-3.
［８０］R.R. Khayrutdinov and V.E. Lukash, J. Comput. Phys. 109,
193 (1993).
［８１］S.C. Jardin et al., J. Comput. Phys. 66, 481 (1986).
［８２］A. Matsuyama et al., Bull. Am. Phys. Soc. 60, 295 (2015).
［８３］A. Matsuyama et al., Bull. Am. Phys. Soc. 58, 261 (2013).
［８４］N. Aiba et al., Plasma Phys. Control. Fusion 58, 045020(2016).
［８５］N. Aiba et al., Nucl. Fusion 57, 022011 (2017).
［８６］N. Aiba et al., Nucl. Fusion 57, 126001 (2017).
［８７］N. Aiba et al., Phys. Rev. Lett. 114, 065001 (2015).
［８８］J. Shiraishi et al., Sci. Rep. 6, 25644 (2016).
［８９］Y. Suzuki et al., Nucl. Fusion 46, L19 (2006).
［９０］Y. Todo et al,. Plasma Fusion Res. 5 S2062 (2010).
［９１］K. Ichiguchi et al., Plasma Fusion Res. 11, 2403035 (2016).
［９２］S. Sakakibara et al., Nucl. Fusion 53, 043010 (2013).
［９３］K. Ichiguchi et al., Proc. 26th IAEA Fusion Energy Conf.
(Kyoto, Japan, 2016) TH/P1-4.
［９４］Y. Todo et al., Plasma Fusion Res. 5, S2062 (2010).
［９５］M. Sato et al., Nucl. Fusion 57, 126023 (2017).
［９６］A. Bierwage et al., Nucl. Fusion 56, 106009 (2016).
［９７］A. Bierwage et al., Phys. Rev. Lett. 114, 015002 (2015).
［９８］A. Bierwage and Ph. Lauber, Nucl. Fusion 57, 116063(2017).
［９９］A. Bierwage et al., Plasma Phys. Control. Fusion 59,
125008 (2017).
［１００］A. Bierwage et al., Plasma Phys Control. Fusion 61,(1989).
［３０］M. Honda et al., Nucl. Fusion 57, 116050 (2017).
［３１］G.G. Plunk et al., New J. Phys. 14, 103030 (2012).
［３２］T. Tatsuno et al., Phys. Plasmas 19, 122305 (2012).
［３３］龍野智哉：日本物理学会誌 67, 640 (2012).
［３４］宮川大和 他：プラズマ・核融合学会誌 93, 271 (2017).
［３５］M. Nunami et al., 20th International Stellarator-Heliotron
Workshop (Greifswald, Germany, 2015) P1S2-25.
［３６］T. Görler et al., Phys. Plasmas 21, 122307 (2014).
［３７］R. Numata et al., J. Comput. Phys. 229, 9347 (2010).
［３８］R. Numata and N.F. Loureiro, J. Plasma Phys. 81,
305810201 (2015).
［３９］K. Alvelius, Phys. Fluids 11, 1880 (1999).
［４０］T. Pianpanit et al., Plasma Fusion Res. 11, 2403040 (2016).
［４１］S. Ishiguro and H. Hasegawa, J. Plasma Phys. 72, 1233
(2006).
［４２］H. Hasegawa and S. Ishiguro, Plasma Fusion Res. 12,1401044 (2017).
［４３］V. Vahedi and M. Surendra, Comput. Phys. Commun. 87,
179 (1995).
［４４］K. Nanbu, Phys. Rev. E 55, 4642 (1997).
［４５］G.A. Emmert et al., Phys. Fluids 23, 803 (1980).
［４６］M. Honda et al., Nucl. Fusion 57, 116050 (2017).
［４７］N. Hayashi et al., Nucl. Fusion 58, 066001 (2018).
［４８］T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
［４９］S. Maeyama et al., Plasma Fusion Res. 11, 2403011 (2016).
［５０］S. Maeyama et al., Nucl. Fusion 57, 066036 (2017).
［５１］S. Satake et al., Phys. Rev. Lett. 107, 055001 (2011).
［５２］Y. Idomura et al., Nucl. Fusion 49, 065029 (2009).
［５３］K.C. Shaing et al., Plasma Phys. Control. Fusion 51, 035009(2009).
［５４］K.C. Shaing et al., J. Plasma Phys. 81, 905810203 (2015).
［５５］M. Yagi et al., Contrib. Plasma Phys. 54, 363 (2014).
［５６］N. Miyato, Plasma Fusion Res. 9, 1203148 (2014).
［５７］N. Miyato et al., Phys. Plasmas 11, 5557 (2004).
［５８］M. Muraglia et al., Phys. Rev. Lett. 107, 095003 (2011).
［５９］S. Maeyama et al., Phys. Rev. Lett. 114, 255002 (2015).
［６０］A. Bierwage et al., Nature Communications 9, 3282 (2018).
［６１］K. Imadera, 20th ITPA Transport and Confinement Topical Group Meeting, NFRI (Daejeon, Korea 2018).
［６２］P. Mantica et al., Phys. Rev. Lett. 107, 135004 (2011).
［６３］S. Yi, et al., Nucl. Fusion 55, 092002 (2015).
［６４］Y.W. Cho et al., 2017 KSTAR conference (Korea, 2017).
［６５］J.M. Kwon et al., Nucl. Fusion 52, 013004 (2012).
［６６］M. Yagi et al., 14th BPSI meeting, RIAM, Kyushu Univer-

014025 (2019).
［１０１］H. Miura et al., J. Comput. Phys. 316, 385 (2016).
［１０２］H. Miura et al., Proc. 26th IAEA Fusion Energy Conf. (Kyoto,
Japan, 2016) TH/P1-5.
［１０３］H. Miura et al., Nucl. Fusion 57, 076034 (2017).
［１０４］H. Miura and N. Nakajima, Nucl. Fusion 50, 054006 (2010).
［１０５］A. Matsuyama et al., Nucl. Fusion 57, 066038 (2017).
［１０６］M. Yagi et al., Proc. 26th IAEA Fusion Energy Conf. (Kyoto,
Japan, 2016) TH/P3-21.
［１０７］A. Matsuyama et al., Nucl. Fusion 54, 123007 (2014).
［１０８］T.E. Evans et al., Nature Phys. 2, 419 (2006).
［１０９］A.B. Rechester and M.N. Rosenbluth, Phys. Rev. Lett.
40, 38 (1978).
［１
１０］R. Kanno et al., Contrib. Plasma Phys. 56, 592 (2016).

９７

Journal of Plasma and Fusion Research Vol.95, No.2 February 2019［１
１１］R. Kanno et al., Nucl. Fusion 58, 016033 (2018).
［１
１２］N. Asakura et al., Nucl. Fusion 53, 123005 (2013).
［１
１３］K. Hoshino et al., Plasma Fusion Res. 9, 3403070 (2014).
［１
１４］K. Hoshino et al., Plasma Fusion Res. 12, 1405023 (2017).
［１
１５］Y. Sakamoto et al., Proc. 25th IAEA Fusion Energy Conf.
(St. Petersburg, Russia, 2014) FIP/3-4Rb.
［１
１６］Y. Homma et al., Nucl. Mater. Energy 12, 323 (2017).
［１
１７］L. Chone et al., Phys. Plasmas 21, 070702 (2014).
［１
１８］G.Y. Park et al., Phys. Plasmas 22, 032505 (2015).
［１
１９］B.D. Dudson et al., J. Plasma Phys. 22, 365810104 (2015).
［１２０］H. Seto et al., Plasma Fusion Res. 11, 1203122 (2016).
［１
２１］A. Fukano et al., J. Nucl. Mater. 211, 363 (2007).
［１
２２］K. Ohya, Physica Scr. T124, 2006 (2006).
［１
２３］Y. Homma and A. Hatayama, J. Comput. Phys. 231, 3211
(2012).
［１
２４］Y. Homma and A. Hatayama, J. Comput. Phys. 250, 206(2013).
［１
２５］Y. Homma et al., Nucl. Fusion 56, 036009 (2016).
［１２６］S. Yamoto et al., J. Nucl. Mater. 463, 615 (2015).
［１
２７］S. Yamoto et al., Contrib. Plasma Phys. 56, 646 (2016).
［１２８］S. Yamoto et al., Proc. 26th IAEA Fusion Energy Conf.
(Kyoto, Japan, 2016) TH/P6-23.
［１２９］S. Yamoto et al., Nucl. Fusion 57, 116051 (2017).
［１３０］M. Yagi, J. Plasma Phys. 81, 435810201 (2015).

Science & Engineering (Jeju, Korea 2017).
［１
３３］A.M. Ito et al., J. Nucl. Mater. 463, 109 (2015).
４］A.M. Ito et al., Nucl. Fusion 55, 073013 (2015).
［１
３
［１
３
５］A.M. Ito et al., Proc. 26th IAEA Fusion Energy Conf. (Kyoto,
Japan, 2016) MPT/P5-23.
［１
３６］A.M. Ito et al., Nucl. Mater. Energy 12, 353 (2017).
［１
３７］A.M. Ito et al., Physica. Scr. T159, 014062 (2014).
［１
３８］K. Ohsawa et al., J. Nucl. Mater. 458, 187 (2015).
［１
３９］T. Yamamoto et al., J. Alloys Compd. 213/214, 533 (1994).
［１４０］T. Yamamoto et al., J. Alloys Compd. 269, 162 (1998).
［１４１］M. Yamawaki et al., Adv. Sci. Tech. 94, 32 (2014).
［１４２］M. Yamawaki et al., J. Nucl. Sci. Tech. 26, S10312 (2015).
［１
４３］T. Suzudo et al., J. Nucl. Mater. 467, 418 (2015).
［１
４４］T. Suzudo and A. Hasgawa, Sci. Rep. 6, 36738 (2016).
［１
４５］A. Takayama et al., Jpn. J. Appl. Phys. 52, 01AL03 (2013).
［１
４６］A. Takayama et al., J. Nucl. Mater. 463, 355 (2015).
［１
４７］Y. Oda et al., Plasma Fusion Res. 9, 340117 (2014).
［１
４８］G.D. Samolyuk et al., J. Nucl. Mater. 418, 174 (2011).
［１
４９］N. Galamba et al., J. Chem. Phys. 126, 204511 (2007).
［１
５０］T. Oda et al., Comput. Mater. Sci. 111, 54 (2016).
［１５１］T. Hoshino, Desalination. 359, 59-63 (2015).
［１５２］G.H. Vineyard, J. Phys. Chem. Solids. 3, 121 (1957).
［１５３］E. Artacho, Phys. Status. Solidi B 215, 809 (1999).
［１５４］E. Fromm and G. Hörz, Intern. Metals Rev. 25, 269 (1980).
［１５５］S. Satake et al., Fusion Sci. Tech. 68, 640 (2015).
［１
５６］K. Hosaka et al., Proc. 28th International Conference on Parallel Computational Fluid Dynamics (Kobe, Japan, 2016).
［１
５７］T. Takizuka et al., J. Comput. Phys. 25, 205 (1977).
［１
５８］K. Nanbu, IEEE Trans. Plasma Sci. 28, 971 (2000).

［１
３
１］矢木雅敏：原子力・量子・核融合事典，第１分冊原子核物理とプラズマ物理・核融合 ISBN978-4-621-088746, I-170 (丸善出版，2014).
［１３２］M. Yagi et al., Proc. M&C International Conference onMathematics & Computational Methods Applied to Nuclear

９８