኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃᱌ෆ㸺኱つᶍ⛉Ꮫィ⟬ࢩࢫࢸ࣒㛵㐃ᴗົࡣࠊࢧ࢖ࣂ࣮ࢧ࢖࢚ࣥࢫࢭࣥࢱ࣮ᮏ㤋ෆࡢ᝟ሗ㒊᝟ሗᇶ┙ㄢࡀᢸᙜࡋ࡚࠸ࡲࡍࠋ㸼KWWSZZZVVFFWRKRNXDFMS

㝵

ಀ䞉ᐊྡ

฼⏝┦ㄯᐊ

㟁ヰ␒ྕ䠄ෆ⥺䠅 
HPDLO  
VRGDQ#FFWRKRNXDFMS┦ㄯဨ୙ᅾ᫬

  䝃䞊䝡䝇᫬ 㛫ᖹ᪥ィ⟬ᶵ฼⏝඲⯡䛻㛵䛩䜛┦ㄯ 䡚኱ุ䝥䝸䞁䝍䚸฼⏝⪅➃ᮎ➼䛾 䡚฼⏝୺䛺䝃䞊䝡䝇ෆᐜ୍㝵ྛ䝉䞁䝍䞊ᗈሗ䛾㜀ぴ฼⏝⪅ㄯヰᐊ

 

ᒎ♧ᐊ ศᩓ䡶䢙䢇䢛䡩䡬䡼༤≀㤋 ぢᏛ䜢䛤ᕼᮃ䛾᪉䛿ඹྠ฼⏝ᨭ᥼ಀ䜎䛷䛤㐃⤡䛟䛰䛥䛔䚹Ṕ௦䛾኱ᆺィ⟬ᶵ➼䛾ᒎ♧䠖䡚䠖ྍど໬ᶵჾᐊ

 ୕ḟඖྍど໬䝅䝇䝔䝮䛾฼⏝䡚

⥲ົಀ  
VRP#FFWRKRNXDFMS⥲ົ䛻㛵䛩䜛䛣䛸䡚

఍ィಀ  
NDLNHL#FFWRKRNXDFMS఍ィ䛻㛵䛩䜛䛣䛸䚸㈇ᢸ㔠䛾ㄳồ䛻㛵䛩䜛䛣䛸䡚

䠄䠅ඹྠ฼⏝ᨭ᥼ಀ
  ཷ㻌 ௜ XNHWXNH#FFWRKRNXDFMS฼⏝ᡭ⥆䛝䚸฼⏝┦ㄯ䚸ㅮ⩦఍䚸䝷䜲䝤䝷䝸䚸ぢᏛ䚸䜰䝥䝸䜿䞊䝅䝵䞁䛻㛵䛩䜛䛣䛸䡚

ඹྠ◊✲ᨭ᥼ಀ

  
UVVHF#FFWRKRNXDFMSඹྠ◊✲䚸ィ⟬ᶵ䝅䝇䝔䝮䛻㛵䛩䜛䛣䛸

䡚

䝛䝑䝖䝽䞊䜽ಀ  
QHWVHF#FFWRKRNXDFMS䝛䝑䝖䝽䞊䜽䛻㛵䛩䜛䛣䛸䡚

᝟ሗ䝉䜻䝳䝸䝔䜱ಀ  
LVHFXULW\#JUSWRKRNXDFMS᝟ሗ䝉䜻䝳䝸䝔䜱䛻㛵䛩䜛䛣䛸䡚

ᅄ◊✲㛤Ⓨ㒊
㝵  ஬㝵 ➃ᮎᶵᐊ

 

3& ➃ᮎᶵ䠄; ➃ᮎ䠅⮬㈍ᶵ䡚


୕
㝵


䡚㸦 㸧ෆࡣᮾ໭኱Ꮫෆࡢࡳࡢෆ⥺㟁ヰ␒ྕ࡛ࡍࠋ㟷ⴥᒣ࣭ᕝෆᆅ༊௨እ࠿ࡽࡣ㢌࡟  ࢆຍ࠼ࡲࡍࠋ

ᮏㄅ䛾ྡ๓䛂㻿㻱㻺㻭㻯䛃䛾⏤᮶㻌᫛࿴  ᖺ䛻ᮾ໭ᆅ༊䛾᭱ึ䛾㟁Ꮚィ⟬ᶵ䛸䛧䛶䚸ᮾ໭኱Ꮫ㟁Ẽ㏻ಙ◊✲ᡤ䛻䛚䛔䛶᏶ᡂ䛥䜜䛯䝟䝷䝯䝖䝻䞁ᘧィ⟬ᶵ䛾ྡ๓䛷 6(1$& 6(1GDL$XWRPDWLF&RPSXWHU 䛛䜙䛸䛳䛶࿨ྡ䛥䜜䛯䚹㻌

— 1 —

SENAC Vol. 51, No. 1（2018. 1）

[巻頭言]
［巻頭言］

情報科学の研究推進における超高速コンピューティングの役割

情報科学研究科

研究科長
徳山

豪

情報科学と最先端のコンピュータは切っても切り離せない関係にあり、私が所属する情報科学
研究科は、コンピュータインフラストラクチャの基盤技術においての研究開発の担い手であると同時に、ヘビーな利用者でもあります。組織的な側面では、サイバーサイエンスセンターから協力講座としての参画をいただいて、ハイパフォーマンスコンピューティングや超高速ネットワー
ク関係の研究の核となっていただいております。一方で利用者としては、機関（部局）単位で利用負担金を支払う制度を活用して、本研究科から 103 名の研究者が現在サイバーサイエンスセン
ターの利用登録をしており、スーパーコンピュータ及び並列コンピュータでのデータ処理や計算機実験を実施しています。科学シミュレーションや、アルゴリズムやシステムの計算機実験などの教育研究はもちろん、毎日数億と発信される Twitter データの解析や、東北メディカルメガバ
ンクのゲノムデータ解析、交通情報解析、更には深層学習に代表される機械学習など、研究科が取り組んでいる学術の推進や、実社会の重要な問題解決における研究活動には、巨大計算が日常的に要求され、スーパーコンピュータの力は必要不可欠であります。
実際、計算機の高速化と情報科学の最先端技術は相互に刺激しあって発展し、世の中のいろい
ろな難問を解決していきます。私自身の専門はアルゴリズム理論ですが、1990 年代には、データ
数 n に対して 10000 log2 n の計算時間の理論保証を持つ新しい探索アルゴリズムを開発しても、それは理論だけで実用システムでは活用されませんでした。なぜならデータ数が十万程度だと、10000 log2 n は十万より大きいので、自明な逐次探索より効率性は上がらず、ただでさえ複雑な理論的アルゴリズムを用いるメリットがなかったからです。しかし、計算機が進歩し、テラサイ
ズ、つまり 10 の 12 乗のデータを扱えるようになると、この方法は自明な探索より 200 万倍以上速くなり、ペタサイズだと、20 億倍の効率化になります。つまり、計算機の能力が向上すると、より理論的に高度なアルゴリズムが活躍する、逆にそのような高度なアルゴリズムが計算機の更なる高速化を支えるという相互作用で、情報社会の大変革を推進しているわけです。
スーパーコンピュータというと、事業種分けの時の蓮舫氏の『世界一になる必要の理由は何ですか？二位じゃダメなんでしょうか？』という発言で、世界最高速への挑戦がともかく注目されました。オリンピックの百メートル走や自動車の F1 レースと同じで、もちろん技術者にとっては一位か二位かは大きな違いがあり、当時この発言に奮起した研究者たちが京コンピュータを世界最高速にして留飲を下げたのでした。
しかし、ビッグデータや IoT を使って第四次産業革命を推進する時代には、汎用的な高速計算が必須の要請となっています。このような時代のスーパーコンピュータは、広く世の中一般のユーザが、自由自在に活用してそれぞれの課題を解いていくのに使われないといけません。利用者の立場からしますと、レーシングサーキットでしか走れず、トップレーサ―しか運転できない車

— 2 —

SENAC Vol. 51, No. 1（2018. 1）のようなスーパーコンピュータは困るわけで、「より使いやすく、実用上重要な問題が効率よく解ける」ことが、世界最高速を競うコンテストでのベンチマークテストに特化する計算機を構築す
るよりも重要だと思いますし、世の中の期待でもあります。その点で、サイバーサイエンスセン
ターのスーパーコンピュータは汎用性が高く使いやすく、また、サポート体制も整備されていて、世の中の動向にあったシステム構築をしており、大いに感謝しております。
一方で技術革新は日進月歩であり、脳神経回路網の計算機構としての実現や、量子計算、量子アニーリング、スピントロニクス、分子コンピューティングなど、新たな科学を取り入れた計算機がどんどん出現し実用化されるのがこれからの 10 年間であると予想しております。それに追従するだけでなく、世界をリードするような画期的な計算システムが東北大学から出現する、それが待望されるところです。
情報科学研究科を核として 2017 年度に開始したデータサイエンス国際共同大学院では、このような技術革新を目指して情報科学を開拓し、ビッグデータを利活用して社会を推進する人材の育成を目標にしております。一方でその教育において、世界最高水準の高速計算の利用が大きな
基盤インフラストラクチャとなっており、今後とも Win-Win の協力関係を強化しつつ、サイバ
ーサイエンスセンターにおける超高速コンピューティングの維持及び躍進に大きく期待しております。

— 3 —

SENAC Vol. 51, No. 1（2018. 1）［共同研究成果］

X 線自由電子レーザーによる分子の超高速崩壊のシミュレーションモデル
中村公亮 1，落合宏平 1，花崎浩太 1，菅野 学 1，山崎 馨 2，
高梨 司 3，福澤宏宣 3，上田 潔 3，永谷清信 4，河野裕彦 1
1
東北大学大学院理学研究科化学専攻
2
東北大学金属材料研究所
3
東北大学多元物質科学研究所
4
京都大学大学院理学研究科物理学・宇宙物理学専攻
近年，X 線自由電子レーザー（XFEL）が生み出す高強度のフェムト秒 X 線パルスが利用できるようになり，物質科学に新たなイノベーションをもたらす光として期待されている。たとえば，
時間分解 X 線回折や化学反応の実時間追跡など様々な実験が行われるようになってきた。高強度の X 線に試料が曝されると，内殻イオン化に続いて起こるオージェ効果によって数フェムト秒のオーダーで多価イオンが生成し，クーロン爆発によって試料が損傷を受ける。一瞬にして電子を剥ぎ取られ多価イオン化した分子のクーロン爆発実験で測定される解離イオンの運動エネルギー分布や解離イオン間の角度相関は，分子の瞬間的な幾何構造の情報を有しており，分子イメージ
ングへの応用などの研究も始まっている。クーロン爆発過程の理解やそれを利用した分子イメー
ジング法の確立，さらには放射線損傷の機構解明を目指して，我々は内殻イオン化からオージェ崩壊による正電荷の発生と生成した電子励起状態からの電子緩和をモデル化し，原子核の動きを
半経験的電子状態計算法である密度汎関数強束縛法で評価する高速の動力学シミュレーション法を開発した。本稿では，この動力学モデルを解説し，5-ヨードウラシルの実験結果との比較を紹介する。

1. はじめに
X 線の位相を揃えてレーザー発振する X 線自由電子レーザー（X-ray Free Electron Laser）[1]の実用化により，X 線波長領域の強力なコヒーレント光源が実現している。電子ビームを相対論的な速度まで加速し，アンジュレータを通過させることにより，フェムト秒のパルス幅と高いピー
ク輝度をもつ可干渉性を有する X 線を発生させることができる。米国 SLAC 国立加速器研究所のLinac Coherent Light Source (LCLS) [2] ， 日 本 の SPring-8 Angstrom Compact Free Electron Laser(SACLA) [3]，ヨーロッパの複数施設[4]の光源が稼働中であり，そのほかの国でも建設が進められている。放射光より数桁高い輝度の XFEL パルス[1]を使うと，非結晶試料からでも X 線回折像が得られ，単一分子レベルでの構造が決定できると期待されている。時間分解の X 線回折も進んでおり[5]，高エネルギーのポンプ‐XFEL プローブ分光のスキームを使って[6-9]，化学反応等フェ
ムト秒オーダーの動力学を原子スケールの分解能で実時間測定する道が拓けつつある[10]。
しかしながら，高強度の X 線に試料が曝されると，内殻イオン化に続いて起こる一連のオージ
ェ崩壊によって数フェムト秒のオーダーで多価イオンが生成し，クーロン爆発によって試料が損傷を受ける。たとえば，ヨウ素などの重原子は XFEL に敏感に反応し，X 線光子を吸って内殻軌
道に 1 つの正孔ができると，原子内で連続的なオージェ崩壊（オージェカスケード）が引き起こされる。その後，重原子に局在した正電荷が速やかに分子全体に行き渡り[11, 12]，分子オージェ
崩壊によって多重イオン化が進行する。その結果，ヨードメタン CH3I などでは，+10 価に至る多価イオンが容易に生成し[11]，分子中の原子が正電荷を帯びてばらばらに飛散し（クーロン爆発），
解離イオンが生成する。数個の X 線光子を吸収すると，このようなイオン化過程が繰り返され，数十価のカチオンまでも生成する[12]。したがって，構造だけでなく，X 線回折に重要な寄与をす— 4 —

SENAC Vol. 51, No. 1（2018. 1）

る内殻の電子状態も大きく変化しうる。XFEL を使った時間分解 X 線回折が直面するこのような問題を克服する 1 つの方法として，10 fs 程度の超短 XFEL パルスを利用することが提案された。
2000 年，Neutze らはモデル計算により，フェムト秒オーダーの超短パルスを使うことで放射線損傷が測定に及ぼす影響を回避し得ることを予見し，X 線超短パルスを使った測定の端緒となった[13]。現在では，構造や電子状態の大きな変化が起こる前に X 線の回折信号を取得する diffractbefore-destroy scheme が提案，応用されている[14]。個別の測定結果の解釈においても試料の放射線損傷過程の理論的解明は極めて重要である。

図 1 XFEL によるヨードメタン CH3I のクーロン爆発過程。X 線領域の光子を吸収したヨウ素原子の内殻イオン化後，原子内での連続的な原子オージェ過程が起こる（オージェカスケード）。発
生した正電荷が分子全体に広がる間に，さらに分子オージェ過程で電子が放出される。この間，
原子核は電子緩和やクーロン反発により大きな運動エネルギーを獲得し，クーロン爆発に至る。
X 線強度が高くなるにつれて複数の光子を吸収しながらオージェカスケードを繰り返し，+20 価を越えるような超多価の親イオンもできる。一般に，親イオンの価数が高くなるにつれて，大きな運動エネルギーをもった解離イオンが生成する。
一方，XFEL によるクーロン爆発自体の研究も，解離イオンの運動量を利用した分子イメージング法の確立や放射線損傷の機構解明に寄与するものとして注目されている。同一分子から放出されたイオンを捉えるコインシデンス計測によって，分子の幾何構造を反映すると考えられる放出イオン間の角度などが'
測定できる。以上のような状況から，X 線光子の吸収に始まる内殻イオン化からオージェカスケードによる多重イオン化の全過程における電子や原子核の動きを理論的に予想することが求められている。しかしながら，そのためには無数の電子配置の関与やそれに伴う原子核の様々な動きを考慮する必要があり，その第一原理的な追跡は容易ではない（第4 節参照）。数原子以上からなる分子のクーロン爆発の追跡
図 2 5-ヨードウラシルの構造。
と，解離イオンの運動量分布の第一原理的な再現あるいは予
クーロン爆発に際して放出され言は，達成しなければならない大きな課題になっている。
る原子イオンの運動量の方向は，そこで，我々は，内殻イオン化後の多重イオン化過程を時それらの間の安定構造における間と共に正電荷が 1 つ 1 つ上昇していく逐次イオン化モデル
角度を反映する。たとえば，図のの導入によって記述し，イオン化に際して生成した電子励起
ヨウ素と酸素の放出イオン間の状態からの電子緩和の効果を原子核への余剰運動量注入に角度は安定構造の中心からの角
よってモデル化する高速動力学シミュレーション法を開発
度'を反映する。1 つの分子からしてきた[15, 16]。このようなモデルのもとで電荷が上昇して放出される複数のイオンの運動
いく分子の電子状態と原子核の動きを，密度汎関数（Density
量はコインシデンス計測によっfunctional theory: DFT）法に基づく半経験的電子状態計算法でて決定できる。文献[15]より改変
ある密度汎関数強束縛（Density-functional based tight-binding:
して転載。
DFTB）法[17, 18]によって評価し，正電荷間の反発だけでな

X 線自由電子レーザーによる分子の超高速崩壊のシミュレーションモデル

— 5 —く化学結合の効果も計算負荷をかけないで適切に組み込めるようにした。これらによって，様々な正電荷をもたらす多重イオン化とそれらのクーロン爆発に対して実験と比較しうる統計サンプ
リングが可能となり，解離イオンの運動量分布（運動エネルギー分布）や角度分布を実験と比較することができるようになった[15, 16]。本稿では，開発したシミュレーション法を概説し，放射線増感剤として使われる 5-ヨードウラシルなどを対象に，解離イオンの運動量分布や角度相関を求め，その有用性を実験結果[15, 19]と比較しながら検証した[15]。開発した動力学モデルの骨子と実験結果との比較を以下に報告する。

2.

クーロン爆発の動力学モデル

重原子含有分子が XFEL によってクーロン爆発に至る過程は，次の 3 つから構成される。
(i) ヨウ素のような X 線に敏感に反応する重原子の X 線内殻イオン化とその後の原子オージェ
カスケードによるイオン化の進行。
(ii) 重原子に局在した正電荷の分子全体への移動とさらなるオージェ崩壊による多重イオン化。
(iii) イオン化の進行に際して生成する電子励起状態からの電子緩和（振動エネルギーへの転換）。
過程(i)が引き金となり，(ii)や(iii)の過程を伴いながら，最終的に原子イオン（あるいは分子イオン）が解離種（フラグメント）として放出される。
ヨードメタン CH3I を標的とした SACLA の実験[11]では，光子エネルギー5.5 keV，パルス長約 10 fs，ピークフルエン
図 3 光子エネルギー5.5 keV，パスが平均 26 J m−2（ピーク光強度 I ~ 3×1017 Wcm−2）のルス長 10 fs，平均ピークフルエン
XFEL パルスによって，+10 価程度の電荷を持つ親イオンがス 26 J m−2 の XFEL パルスをヨ
最も多く生成した（このような多価の親イオンは最終的にードメタン CH3I に照射した際に生はクーロン爆発する）。ヨウ素原子が 3 光子を吸収すると成する親イオンの電荷 Qfin の分布。
+20 価を越える親イオンまで生成する（最大で+22 程度）[16]。同様のパルスに対して，ジヨードメタン CH2I2 も+10価を持つ親イオンが最も多く生成し，やはり+20 価を越える親イオンまで生成する[16]。2 つのヨウ素原子が独立に X
線光子を吸収するため，+10 価を越える親イオンの相対量は CH3I の場合より多くなる。このように，親イオンの正電
荷 Qfin（≤ 22）は広い分布をもち，Qfin が大きいほど，大きな運動エネルギーを持つ解離イオンが放出される。
ヨウ素含有分子の正電荷 Q(t)は，時間 t の関数と見なせ
図 4 ヨウ素含有分子の XFEL パる。ヨウ素原子が X 線光子を吸収して内殻イオン化が起こ
ルスによる内殻イオン化後の電荷った時刻を t = 0 とすると，オージェカスケードによってQ(t)が上昇し，最終電荷 Qfin に達すると考える。これまで， 上昇モデル。はオージェカスケー
ドの時間スケールを特徴づけるパ
電荷上昇の時間スケールを導入した次式が提案されてきラメーターである。逐次垂直イオ
た[11, 15]。
ン化モデルでは，分子の電荷 Q(t)(1)
Q (t )  Qfin 1  e  t /が整数値になったところでイオン
はパルス長とオージェ崩壊の寿命の両者の効果を含んで 化したと仮定して，分子の電子数いる。オージェカスケードは超高速で起こり，実験的・理を 1 つ減らした電子状態に瞬時に論的にはおおよそ 10 fs のオーダーと見積もられている
スイッチさせる。



— 6 —

SENAC Vol. 51, No. 1（2018. 1）[12, 20]。ヨウ素原子に局在して生成した正電荷は分子全体に移動していくが，この時間スケール
 M はより短いことが示唆されている。たとえば，5-ヨードウラシルでは，点電荷古典クーロンモ
デル（次節参照）を使った実験結果の解析から， Mが 2 fs 程度と見積もられている[15, 19]。
我々が開発した動力学モデルは，以上の物理的考察に基づき，以下の 3 点を骨子としている[15,16]。
(i) 逐次垂直イオン化モデル
まず，生成した Q(t)の電荷が一瞬にして分子全体に広がると仮定する。つまり，(1)式に沿ってQ(t)が増加し，たとえば，Q(t) = 3 のように整数電荷になった時刻 t で分子をその構造を保ったま
ま 2 価から 3 価にイオン化させる。Q(t) = Qfin になるまでこの手続きを繰り返す。
(ii) 密度汎関数強束縛（DFTB）法による動力学
逐次垂直イオン化モデルに従って多重イオン化する分子の構造は，電子の脱離や電荷を帯びた
原子間のクーロン反発で，初期構造から大きく膨張していく。この動力学過程における原子核の動きには，半経験的電子状態計算法である DFTB 法を使う。この方法は，高速ながらも DFT 法に近い精度で計算でき，化学結合やクーロン反発の効果を取り入れることができる。
DFTB 法では，系の電子密度を中性原子で構成される電荷の偏りがない参照系の電子密度 0 と
電子密度揺らぎの和で表す。つぎに，の汎関数である DFT 法のエネルギーを 0 の近傍での多項式として展開する。1 次の項は 0 となり，系のエネルギーは参照系の Kohn-Sham ハミルトニ
アンと原子間反発から得られる 0 次の主要項，およびに依存する 2 次以上の補正項に分けられる。電荷の偏りを無視して 0 次解を評価する手法は Non-self-consistent-charge (Non-SCC) DFTB 法とよばれている。強束縛近似（化学の用語では LCAO）に基づいて，参照系の Kohn-Sham 軌道を
原子軌道の線形結合で表すと，Kohn-Sham 方程式は原子軌道間の行列要素からなる永年方程式に帰着される。
DFTB 法は，これらの行列要素を事前に DFT 法で見積もられたパラメーターにして，計算の高速化を図っている。一方，系のエネルギーに対して 2 次の補正項も考慮する手法は Selfconsistent-charge DFTB (SCC-DFTB) 法[17]とよばれている。SCC-DFTB 法では，電荷とエネルギーを自己無撞着的に求める。それでもパラメーターを使用しない DFT 法と比べればはるかに低負荷であり，異核分子や極性分子の電荷分布を適切に記述できる。したがって，電荷の分布が重要となるクーロン爆発を扱う本研究では，Slater-Koster パラメーターセット halorg-0-1 [18]と組み合わせた SCC-DFTB 法を使う。以上の DFTB 法については，SENAC の既報でも解説している[21]。
クーロン爆発に至る動力学過程においては，内殻空孔状態に加えて，その後のオージェ崩壊に
よって価電子の励起状態も生成する。価電子の励起状態は有効電子温度 Te で特徴づけられるフェ
ルミ-ディラック統計に従って分布すると考えられる。本モデルでは，複数の励起状態を経由する
非断熱的なダイナミクスを追跡するのではなく，フェルミ-ディラック分布で平均化された SCCDFTB 法の 1 つのポテンシャル上の動きで評価する。その際，励起状態の分布を適切に与える有効電子温度 Te をパラメーターとして導入する。一方，多重イオン化が始まる 2p 空孔状態などの
内殻空孔状態は無視している（DFTB 法には内殻空孔状態は含まれていない）。内殻空孔状態の寿命は 2 ~ 3 fs 程度と短く[20]，かつ数 fs で価電子励起状態への電子緩和が進むので，最終的に生成する解離イオンの運動エネルギーに与える影響は相対的には小さいと考えられる。
(iii) 逐次イオン化に際する余剰振動エネルギーの発生
分子は垂直イオン化に際して様々な電子励起状態に分布することになり，そのポテンシャル形
状の違いを有効電子温度の導入によって平均化することを説明した。いろいろな電子励起状態にある分子は，非断熱遷移によって複数のポテンシャル面を乗り移り，電子緩和を経て余剰の振動エネルギーを獲得する。ここでは，+n の電荷をもつ段階で獲得する余剰エネルギー En をパラメ

X 線自由電子レーザーによる分子の超高速崩壊のシミュレーションモデル

— 7 —ーターとして，この過程を次のような現象論的モデルで扱う。まず，Q(t)が整数値に達して垂直イオン化した瞬間に，原子核の運動エネルギーにその電荷の状態で最終的に獲得すると思われるエ
ネルギーを追加してしまう。具体的には，各イオン化段階で原子核 j の運動量 p j を新たな運動量p j に変える。

j p j  p j
p
(2)



ここで， Δp j は運動量ベクトルの変化量である。 Δp j の方向の選び方としてはランダムにとることも考えられるが，ここでは，結合軸方向に沿って原子同士が及ぼす力を考慮した結合軸反跳モ
デルを採用する。たとえば，CO2 分子の C 原子の 1s 軌道に空孔ができると全対称伸縮振動が誘起され，O の 1s 軌道に空孔ができると逆対称伸縮振動が誘起されることが知られている[22]。これに倣って，このモデルでは，原子 A と B それぞれに対して，化学結合 A-B に沿った方向に正規分布に従った大きさをもつ変化量 ΔpA (A-B) と ΔpB (A-B) を考える。さらに両者の間に，孤立した二原子分子 AB に対して成立する次の運動量保存則が近似的に成り立つとする[15, 16]。



ΔpA (A-B) + ΔpB (A-B) = 0(3) B Δp j ( j -B) を求める。最終的な Δp j(3)式をすべての結合に適用して，1 つの原子 j に対する Δp j =の絶対値は，原子核の運動エネルギーの総和が En だけ増えるようにスケールし直す。

3.

5-ヨードウラシルへの適用

本節では，開発した動力学モデルを放射線増感剤として知られる 5-ヨードウラシルに適用した結果を紹介する。CH3I の実験と同じ XFEL パルス（光子エネルギー5.5 keV，パルス長 10 fs，ピー
クフルエンス 26 J m−2）を 5-ヨードウラシルに照射した際に得られた原子様の解離イオンの運動エネルギー分布を図 5 に示す[15, 19]。これらは様々な電荷を持つ親イオンが解離して生成したものである。質量から予想されるように，最も重いヨウ素イオンの運動エネルギーが最も小さい。

図 5 XFEL パルスによってクーロン爆発した 5-ヨードウラシルの解離原子イオン運動量分布。青線は SCC-DFTB 法に基づいた本モデルの計算結果，赤破線は点電荷古典クーロンモデルによる結果を示している。実験では，最も重いヨウ素イオンの運動エネルギーは小さいが，最も軽い H+の運動エネルギーのピーク値は N+や O+と大きな違いはない。N+や O+の運動エネルギーはそれらよ
り軽い C+の運動エネルギーよりピーク値で比べると 2 倍ほど大きく，運動エネルギーが質量だけでは決まっていないことを示している。また，(e)の運動量分布は，主要な電荷 q が+1 から+4 までのヨウ素イオンの和である。文献[15]より改変して転載。

— 8 —

SENAC Vol. 51, No. 1（2018. 1）

一方，最も軽いプロトン H+の運動エネルギーのピークは窒素イオン N+や酸素イオン O+のピークと大差なく，また，N+や O+のほうがそれより軽い炭素イオン C+より高い運動エネルギーを有しており，解離イオンの運動エネルギーが質量だけで決まっていないことがわかる。
図 5 の実験の運動エネルギー分布を再現するため，第 2 節で説明した動力学モデルを使ったシ
ミュレーションを行った。初期状態は中性分子の電子基底状態とし，温度 T = 300 K で SCC-DFTB法に基づいた古典トラジェクトリーを走らせた。トラジェクトリーからランダムに初期の原子核の位置と速度を選び出し，その時刻を最初に X 線光子を吸収した時刻 t = 0 として，最終電荷 Qfin
でクーロン爆発に至る動力学シミュレーションを始めた。統計平均を取るため，各最終電荷に対して，1000 本の動力学計算を行った。その際，5-ヨードウラシルに対する Qfin の分布は，ヨウ素原子を同じように 1 つ有する CH3I の分布 P(Qfin)と同じとした（図 3）。物理的な考察[11, 20]から，
電荷上昇の時間スケール = 10 fs，電子温度 Te = 6 eV，n 番目のイオンの段階で獲得する余剰エネ
ルギーを n に関わらず En = 6 eV とした。図 5 の青線が SCC-DFTB 動力学モデルから得られた結果である。各イオンの広いエネルギー領域で実験値と一致している。
注目すべきは，シミュレーションの結果も N+や O+の運動エネルギーがそれより軽い C+より高いことである。ウラシル環の外側に位置する O+は，環内の C+より早く放出される。実験やシミュ
レーションの結果は，クーロン爆発の早い段階で放出されるほど，大きなクーロン反発によって高い運動エネルギーを得ることを示している（同じ電荷を持ったイオンであれば，早く放出されたものほど高い運動エネルギーをもつ）。後から放出されるイオンは，先に放出されたイオンの正電荷分だけ弱くなったクーロン反発の環境にあり，運動エネルギーは小さくなる。これは超多価分子イオンのクーロン爆発の一般則と考えられる。ただ 1 種類の原子からなる C60 分子のクーロ
ン爆発でも，早く放出される高エネルギーの炭素イオンと後から放出されるものと大きく分けると 2 成分あることがわかっている（C60 の場合，前者はクーロン反発による解離，後者は統計的解離と見なすこともできる）[23, 24]。
N+はウラシル環内に位置するが，その外側はすぐに放出される軽い H+である。したがって，N+は C+よりも平均的には早く放出され，そのため C+より大きな運動エネルギーをもつことになる。
逆に，C+はその外側の重いヨウ素や酸素のために遅れて放出されることになり，クーロン反発エ
ネルギーが低下した環境の中で放出されるので，N+や O+よりも小さな運動エネルギーをもつともいえる。H+は環の外にあり，しかも軽いので，他の相対的に重い原子からなる分子の残りの部分
から 10 fs もかからずに大きく離れていく。一方，そのような短い時間 t = では，(1)式から，Q(t)は Qfin の 63%の電荷にしか達していない。そのため，H+の運動エネルギーはそれよりずっと重い
N+や O+と同程度になっている。
図 5 には，原子を点電荷で扱い，その電荷間のクーロン反発だけを取り入れた古典動力学モデ
ルの結果も示している（赤破線）。この点電荷古典クーロンモデルの場合も，(1)式を使って各原子の電荷を上昇させている。点電荷古典クーロンモデルには化学結合の効果は全く入っていないので，点電荷間の反発だけでイオンの放出のタイミングが決まり，イオンの質量が軽くなるにつれてより大きな運動エネルギーを持つ傾向が顕著になる。ほぼ同じ質量をもつ C+，N+，O+の運動エ
ネルギーのピークは 20 eV あたりにある。これらのシミュレーションの結果から，クーロン爆発
過程にも分子構造の影響や化学結合の効果が残っており，それらを適切に取り込むことが重要であることがわかる。
図 6 は，各イオンのヨウ素イオンに対する放出角の分布である[15, 19]。2 つの放出イオン間の角度は次式で定義されている。





 
p p
cos    A  B
pA pB

(4)

(4)式の pA と pB の一方をヨウ素イオン Iq+の運動量ベクトルにし，もう一方をコインシデンス実

X 線自由電子レーザーによる分子の超高速崩壊のシミュレーションモデル

— 9 —

験で得られた相手の運動量ベクトルとして分布をとっている。SCC-DFTB 法に基づいた動力学シ
ミュレーションの H+，N+，O+の角度のピークは実験結果とほぼ一致している。これらのピーク
は，図 2 に示された安定構造におけるヨウ素原子と各原子との角度'の値に対応しており（安定構造における角度のおおよその領域を，図 6 では縦帯で示してある），2 つの放出イオン間の角度分布が分子の幾何構造の直接的な情報をもっていることがわかる。つまり，XFEL を使った超高速のクーロン爆発を利用すれば，時々刻々変化する反応中の分子の瞬間的な構造の情報が得られる。

図 6 XFEL パルスによってクーロン爆発した 5-ヨードウラシルの解離原子イオンのヨウ素イオンに対する角度分布。青線は SCC-DFTB 法に基づいた本モデルの計算結果である。(4)式で定義されている角度の分布のピークは，図 2 の安定構造におけるヨウ素原子と他の原子との角度'（縦帯の領域）と相関している。H+，N+，O+のピークは'と一致しているが，C+の場合は明瞭なピー
クはない。
このように XFEL 誘起のクーロン爆発は時間分解分子イメージングの新たな方法としての可能性を有している。しかし，構造を反映する放出イオン間の角度分布は 2 体の相関なので，1 体の物
理量である運動エネルギーの場合に比べて，実験結果を高い精度で理論的に再現することは容易ではない。実際，図 6 の実験と理論の結果では，ピークの幅など微細な構造が異なっている。また，C+は明確なピークを示しておらず，遅い段階で放出されるイオンの角度分布から分子構造の情報を得ることが容易でないことがわかる。どのようにすれば実験で得られた解離イオンの運動
エネルギー分布や角度分布からその瞬間の分子構造を抽出できるかに答えるには，より第一原理的な方法からの理解がまず必要になる。

4.

XFEL によるオージェ崩壊とクーロン爆発に対する第一原理手法開発の現状

これまで説明してきたように，Ｘ線が物質に照射されると，主に重元素の内殻準位からのイオン化が起こり，続いてオージェ崩壊を繰り返し，さらにカスケード状に分岐しながら崩壊する。
また，オージェ崩壊と競合する蛍光 X 線を放出してより安定な状態に移る過程も起こる[25, 26]。
これら一連の過程は，多数の電子の放出を伴い，試料は最終的に正電荷間の斥力によって速やかにクーロン爆発する。特に強力な FEL 光源の下ではＸ線光子の吸収が複数回起こり，オージェ崩
壊が何段階も生じて，数 10 価もの多価イオンを生じる場合もある[12, 20]。このような連続状態を

— 10 —

SENAC Vol. 51, No. 1（2018. 1）

含む高エネルギー過程の計算は，基底状態付近の低エネルギー状態の高精度計算に特化してきた
従来の量子化学的手法では極めて困難である。このような過程を半経験的・現象論的に記述するのが第 2 節で説明した動力学モデルであった。
これに対して，近年，第一原理的に計算する手法が開発され[12, 27-29]，実験結果と比較し得る
結果を挙げつつある。Son らは，波動関数を数値グリッドで表して[30, 31]，ハートリー-フォック-スレイター方程式[32]を解くことにより，任意の内殻占有状態に対する数値波動関数を高速に計算する手法を開発した[27]。モンテカルロ法[33]と結びつけて，上述のカスケード状の動力学を再現する方法を提唱した。XATOM と名付けられたこのプログラムは，最初，原子に適用され，イオ
ン収率や運動エネルギー分布などの実験結果の解釈に使われた[20, 28]。特筆すべきことは，個々の過程を on-the-fly で第一原理計算し，全ての崩壊カスケード分岐をバイアス無く統計的に探査した結果，個々の電子状態計算は平均場近似の範疇でありながら，いくつかの観測量において定量的にも満足しうる結果が得られていることである。彼らは次いで，この手法で得られた数値的（原子）波動関数を基底関数とした分子の計算法（XMOLECULE）を開発し[29]，CH3I 分子などの内殻イオン化からクーロン爆発に至る一連の放射線損傷過程の計算を実現している[12]。
このような第一原理計算の結果を実験結果と照合することによって，直接観測が難しい反応の電子・原子核動力学の情報などを得ることができ，多くの系への適用が期待される。しかしながら，XATOM/XMOLECULE においても，さらに取り込むべき効果が残っている。分子のクーロン
爆発は電子励起状態からの原子核動力学であるため，核運動に伴って電子状態の非断熱遷移[34]を経由していくはずであるが，電子状態に関して単参照の XATOM/XMOLECULE にはこの効果は顕に入らない。また，高強度Ｘ線のもとでは，イオン化の速さがオージェ崩壊の速さを上回る現象も観測されており，オージェ崩壊より短い時間スケールを扱わなければならないことが示唆されている[35]。これに対して，XATOM/XMOLECULE は静的な平均場計算の時系列上の組み合わせであるため，過渡状態を含めた真の実時間ダイナミクスを扱えない[36]。これらは複数の電子状態が同時に関与する現象であり，対応できる計算手法をうまく開発できたとしても，相当の計算
コストの増大を伴ってしまう。オージェカスケードからクーロン爆発に至る全過程の包括的な理解には，計算負荷が低く大きな分子にも適用できる半経験的手法から，計算負荷が高くともより
多くの物理過程を正確に記述できる第一原理手法まで，複数の理論手法が必要である。クーロン
爆発の実験データから分子の幾何構造を決定する分子イメージング法の確立も，それらの相補的な組み合わせによって初めて可能になるはずである。

5.

おわりに

本稿では，XFEL 誘起のクーロン爆発をシミュレーションする動力学モデルについて解説し，その 5-ヨードウラシル分子への適用結果を紹介した。開発した動力学モデルでは，内殻イオン化からオージェ崩壊による正電荷の発生と生成した電子励起状態からの電子緩和をモデル化し，原子
核の動きを半経験的電子状態計算法である SCC-DFTB 法で評価した。SCC-DFTB 法では，大きな計算コストをかけることなく，化学結合の効果と分子中の電荷分布を適切に記述できる。そのため，得られた結果は測定された解離イオンの運動エネルギー分布や解離イオン間の角度分布をうまく再現した。シミュレーションと実験結果の比較から，クーロン爆発の早い段階で分子から放出される位置にある原子ほど，高い運動エネルギーをもって原子イオンとして放出されることがわかった。つまり，運動エネルギーがクーロン爆発の時系列を反映している。また，解離イオン間の角度分布も分子の幾何構造を反映しており，これらの測定量から逆に分子の構造を推定する
方法，つまり，分子イメージング法の実現はもはや夢ではない。理論面からの今後の寄与としては，大きな分子にも適用できる半経験的手法からより正確な第一原理手法まで様々な理論手法の開発と相補的な利用をさらに進めることが重要である。

X 線自由電子レーザーによる分子の超高速崩壊のシミュレーションモデル

— 11 —

本研究では DFTB 法を使って，分子の動きを評価した。その計算の一部は，東北大学サイバーサイエンスセンターの並列コンピューターにおいて，DFTB+ 1.2.2 プログラムを用いて実行した。
並列化効率などについては，SENAC に報告している[21]。すでに DFTB+は Ver. 1.2 から 1.3 にアップデートされており，下記の機能などが新たに追加されている。
(1) 電子エネルギーを電荷揺らぎに関して 3 次まで展開した DFTB3 が使える[37, 38]。
(2) 分子間力の分散力補正の方法が増えた。従来の Lennard-Jones と Slater-Kirkwood に加え，
DFTB3 用に最適化した DFT-D3[39]という補正が選べる。
(3) 線形応答時間依存 DFTB 法による励起状態計算が可能になった[40]。
(4) 構造最適化アルゴリズムの種類や力場計算用の数値微分法の種類が増えた。
2017 年 6 月にリリースされた DFTB+の最新版 Ver. 17.1 は，内容的には 1.3 とほとんど変わりは無
いが，誰でも登録なしで入手できるフリーソフトになっており，ワークステーションなどにも簡単に移植できるようになっている。

謝辞
本研究の計算の一部は，東北大学サイバーサイエンスセンターの並列コンピューターで実行した。なお，本研究の計算は科研費（河野：JP16H04091）の助成の下に行われた。

参考文献[1] C. Peliegrini, A. Marinelli, and S. Reiche, “The physics of x-ray free-electron lasers”, Rev. Mod. Phys.
88, 015006 (2016).
[2] Y. Ding et al., “Measurements and Simulations of Ultralow Emittance and Ultrashort Electron Beams
in the Linac Coherent Light Source”, Phys. Rev. Lett. 102, 254801 (2009); C. Bostedt et al., “Linac
Coherent Light Source: The first five years”, Rev. Mod. Phys. 88, 015007 (2016).
[3] T. Ishikawa et al., “A compact X-ray free-electron laser emitting in the sub-ångström region”, Nat.
Photonics, 6, 540 (2012); K. Tono et al., “Beamline, experimental stations and photon beam diagnostics
for the hard x-ray free electron laser of SACLA”, New J. Phys. 15, 083035 (2013).
[4] E. Allaria et al., “Highly coherent and stable pulses from the FERMI seeded free-electron laser in theextreme ultraviolet”, Nat. Photonics 6, 699 (2012).
[5] K. J. Gaffney and H. N. Chapman, “Imaging atomic structure and dynamics with ultrafast x-ray
scattering”, Science 316, 1444 (2007); J. Xu et al., “Time-resolved molecular imaging”, J. Phys. B 49,
112001 (2016).
[6] B. K. McFarland et al., “Ultrafast X-ray Auger probing of photoexcited molecular dynamics”, Nat.
Commun. 5, 4235 (2014).
[7] M. P. Minitti et al., “Imaging Molecular Motion: Femtosecond X-Ray Scattering of an Electrocyclic
Chemical Reaction”, Phys. Rev. Lett. 114, 255501 (2015).
[8] K. H. Kim et al., “Direct observation of bond formation in solution with femtosecond X-ray
scattering”, Nature 518, 385 (2015).
[9] Y. Obara et al., “Femtosecond time-resolved X-ray absorption spectroscopy of anatase TiO2
nanoparticles using XFEL”, Struct. Dyn. 4, 044033 (2017).
[10] J. Ullrich, A. Rudenko, and R. Moshammer, “Free-electron lasers: new avenues in molecular physics
and photochemistry”, Annu. Rev. Phys. Chem. 63, 635 (2012).
[11] K. Motomura et al., “Charge and Nuclear Dynamics Induced by Deep Inner-Shell Multiphoton
Ionization of CH3I Molecules by Intense X-ray Free-Electron Laser Pulses”, J. Phys. Chem. Lett. 6,2944 (2015).

— 12 —

SENAC Vol. 51, No. 1（2018. 1）[12] A. Rudenko et al., “Femtosecond response of polyatomic molecules to ultra-intense hard X-rays”,
Nature 546, 129 (2017).
[13] R. Neutze et al., “Potential for biomolecular imaging with femtosecond X-ray pulses”, Nature 406, 752(2000).
[14] C. H. Yoon et al., “A comprehensive simulation framework for imaging single particles andbiomolecules at the European X-ray Free-Electron Laser”, Sci. Rep. 6, 24791 (2016).
[15] K. Nagaya et al., “Femtosecond Charge and Molecular Dynamics of I-containing organic molecules
Induced by Intense X-Ray Free-Electron Laser Pulses”, Faraday Discuss. 194, 537 (2016); P. Decleva
et al., “Structural dynamics: general discussion”, Faraday Discuss. 194, 583 (2016).
[16] T. Takanashi et al., “Ultrafast Coulomb explosion of a diiodomethane molecule induced by an X-ray
free-electron laser pulse”, Phys. Chem. Chem. Phys. 19, 19707 (2017).
[17] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai, and G. Seifert,“Self-consistent-charge density-functional tight-binding method for simulations of complex materials
properties”, Phys. Rev. B 58, 7260 (1998).
[18] B. Aradi, B. Hourahine, and Th. Frauenheim, “DFTB+, a sparse matrix-based implementation of the
DFTB method”, J. Phys. Chem. A 111, 5678 (2007). 詳細は，http://www.dftb-plus.info/ 参照。
[19] K. Nagaya et al., “Ultrafast Dynamics of a Nucleobase Analogue Illuminated by a Short Intense X-ray
Free Electron Laser Pulse”, Phys. Rev X 6, 021035 (2016).
[20] H. Fukuzawa et al., “Deep Inner-Shell Multiphoton Ionization by Intense X-Ray Free-Electron Laser
Pulses”, Phys. Rev. Lett. 110, 173005 (2013).
[21] 菅野 学，小林倫仁，W. C. Chung, 山崎 馨，瀬高 渉，河野裕彦，結晶性分子ジャイロスコー
プの構造と回転動力学の密度汎関数強束縛法によるシミュレーション，SENAC Vol. 48，No.
4，1 (2015).
[22] M. Hoshino, K. Nakagawa, T. Tanaka, M. Kitajima, H. Tanaka, A. De Fanis, K. Wang, B. Zimmermann,
V. McKoy, and K. Ueda, “Vibrationally resolved partial cross sections and asymmetry parameters foroxygen K-shell photoionization of the CO2 molecule”, J. Phys. B 39, 3047 (2006).
[23] 山崎 馨，上田 潔，河野裕彦，X 線自由電子レーザーパルスによるフラーレン超多価カチオ
ン C60q+の爆発解離の動力学シミュレーション ―2 段階爆発機構の提案―，SENAC Vol. 48，
No. 3，1 (2015).
[24] K. Yamazaki, T. Nakamura, N. Niitsu, M. Kanno, K. Ueda, and H. Kono, “Communication: Two-Step
Explosion Processes of Highly Charged Fullerene Cations C60q+ (q=20‒60)”, J. Chem. Phys. 141,
121105 (2014).
[25] W. Bambynek, B. Crasemann, R. W. Fink, H. U. Freund, H. Mark, C. D. Swift, R. E. Price, and P. V.
Rao, “X-Ray Fluorescence Yields, Auger, and Coster-Kronig Transition Probabilities”, Rev. Mod. Phys.
44, 716 (1972).
[26] R. Santra, “Concepts in x-ray physics”, J. Phys. B 42, 023001 (2009).
[27] S.-K. Son, L. Young, and R. Santra, “Impact of hollow-atom formation on coherent x-ray scattering athigh intensity”, Phys. Rev. A 83, 033402 (2011).
[28] B. Rudek et al., “Ultra-efficient ionization of heavy atoms by intense X-ray free-electron laser pulses”,
Nat. Photonics 6, 858 (2012).
[29] Y. Hao, L. Inhester, K. Hanasaki, S.-K. Son, and R. Santra, “Efficient electronic structure calculation
for molecular ionization dynamics at high x-ray intensity”, Struct. Dyn. 2, 041707 (2015).
[30] A. D. Becke, “A multicenter numerical integration scheme for polyatomic molecules”, J. Chem. Phys.
88, 2547 (1988).
[31] X. M. Tong and S. I. Chu, “Theoretical study of multiple high-order harmonic generation by intense
ultrashort pulsed laser fields: A new generalized pseudospectral time-dependent method”, Chem. Phys.

X 線自由電子レーザーによる分子の超高速崩壊のシミュレーションモデル

— 13 —

217, 119 (1997).
[32] J. C. Slater, “A Simplification of the Hartree-Fock Method”, Phys. Rev. 81, 385 (1951).
[33] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics 3rd Ed.
(Cambridge University Press, New York, 2009).
[34] 非断熱動力学に関しては数多くの総説があるが，J. C. Tully, “Mixed quantum–classical dynamics”,
Faraday Discuss. 110, 407 (1998).
[35] L. Young et al., “Femtosecond electronic response of atoms to ultra-intense X-rays”, Nature 466, 56(2010).
[36] 電子の実時間動力学を含めて扱った試みの 1 つとして，T. Matsuoka and K. Takatsuka, “Dynamics
of photoionization from molecular electronic wavepacket states in intense pulse laser fields: A
nonadiabatic electron wavepacket study”, J. Chem. Phys. 146, 134114 (2017).
[37] M. Gaus, Q. Cui, and M. Elstner, “DFTB3: Extension of the Self-Consistent-Charge Density-Functional
Tight-Binding Method (SCC-DFTB)”, J. Chem. Theory Comput. 7, 931 (2011).
[38] 菱沼直樹，及川啓太，岡田朝彦，菅野 学，山崎 馨，W. C. Chung，齋藤敦子，河野裕彦,
DNA 鎖切断の反応動力学シミュレーションと AVS/Express を用いた三次元可視化, SENAC
Vol. 50, No.1, 3 (2017).
[39] S. Grimme, S. Ehrlich, and L. Goerigk, “Effect of the damping function in dispersion corrected density
functional theory”, J. Comput. Chem. 32, 1456 (2011).
[40] T. Niehaus, S. Suhai, F. Della Sala, P. Lugli, M. Elstner, G. Seifert, and T. Frauenheim, “Tight-binding
approach to time-dependent density-functional response theory”, Phys. Rev. B 63, 085108 (2001).

— 14 —

SENAC Vol. 51, No. 1（2018. 1）［共同研究成果］[共同研究成果]

ベクトル型スーパーコンピュータを用いた
八木・宇田アンテナアレーの高利得化
今野 佳祐 † , 陳 強 †東北大学大学院工学研究科 通信工学専攻 †1

まえがき
電磁波を用いた無線電力伝送に関する研究は 20 世紀初頭から行われており，長い歴史がある[1]-[3]．無線電力伝送は，近傍界を用いたものと遠方界を用いたものに大別される．近傍界を用いた無線電力伝送は高効率だが伝送距離が短く，遠方界を用いた無線電力伝送は低効率だが伝送距離が長いという特徴がある．また，近傍界を用いた無線電力伝送では送受電アンテナ間の整合が重要な一方，遠方界を用いた無線電力伝送ではアンテナの高利得化が重要になる．近傍界を用いた無線電力伝送は，マサチューセッツ工科大学のグループの研究により，非常に高い効率で電力伝送可能なことが実証されて以来，盛んに研究されている [4]．その一方で，遠方界を用いた無線
電力伝送は，宇宙太陽光利用システムなどの送受電システムとしての応用が古くから期待されていた [5]-[7]．
遠方界を用いた無線電力伝送では，アンテナの高利得化を図るために，アレーアンテナが用いられる．例えば，宇宙太陽光利用システムでは，太陽電池で発電した電力を地上に送電するために，数億素子規模のアレーアンテナが利用されることが想定されており，その給電分布の研究 [8]や数値解析法の研究 [9], [10] がこれまでに行われてきた．一般に，アンテナの利得はその大きさとトレードオフの関係にあるので，高利得のアンテナは大型のアンテナになる．その一方で，アンテナの大きさが無制限に許容されることは実用上まれであり，限られた空間で利得が最大のアンテナを実現することが求められる．
本報告では，遠方界を用いた無線電力伝送用のアレーアンテナの高速な設計法を提案する．提案法は，アレーアンテナの大きさを制約条件とし，最大の指向性利得を有するアンテナを実現することができる．アンテナ素子は八木・宇田アンテナとし，アレーアンテナが 0.2 m × 0.2 m ×0.2

m の立方体に収まるように設計する．設計したアレーアンテナを遠方界無線電力伝送システムの受信アンテナとして用い，1 素子当たりの受信電力を明らかにすることで，整流回路の設計指針を
得る．

2

設計法
設計する八木･宇田アレーアンテナの構造を図 1 に示す．動作周波数は 2.45 GHz であり，波長はおよそ 12cm であるから，放射器，導波器，反射器の素子長はいずれも 6cm 前後となる．一般
的に，反射器は半波長以上の長さとなるから，0.2 m × 0.2 m ×0.2 m の立方体にアレーアンテナ
全体を収めるためには，E 面 (図 1 の yz 面) の八木・宇田アンテナ素子数は高々3 素子が限界であ— 15 —

ベクトル型スーパーコンピュータを用いた八木・宇田アンテナアレーの
高利得化

dx

l6
d5

dy
y
z

l5

l3 l2
d2 d1

l1
a

y
x
Front view

Side viewzx

図 1: 4 × 3 素子の八木・宇田アレーアンテナ.
る．その一方で，H 面 (図 1 の xz 面) の八木・宇田アンテナ素子数は素子間距離次第で自由に変えられるが，素子同士が近接しすぎると開口効率が低下したり，素子間相互結合によってインピー
ダンスの不整合が生じるという問題がある．したがって，限られた寸法内で最大の開口効率が実現でき，素子間相互結合の影響もなるべく低減できる構造が望ましい．そこで本報告は，立方体の両側面と内部にそれぞれ 2 素子，合計 4 素子の八木・宇田アンテナ素子が等間隔で配置された構造を考える．アンテナの素子間隔は 0.2/3 m(≈ 0.67λ) であり，一般的なアレーアンテナの素子
間隔 (0.5λ ∼ 1λ) の範囲に収まる．また，両側面に配置したアンテナの開口は立方体の内部だけでなく外部にも存在するので，立方体の内部にアンテナを配置するという制約を満たしつつ，立方体の側面積から算出される開口を上回る実効開口を持たせることが可能になる．八木・宇田アン
テナ 1 素子あたりの導波器は 4 本，放射器と反射器はそれぞれ 1 本とした．
このような構造を持つ 4 × 3 素子の八木・宇田アレーアンテナは，全ての八木・宇田アンテナ素子が同一であるならば，導波器，放射器および反射器の素子長に関するパラメータが l1 , l2 , l3 , l4 , l5 , l6の 6 つ，それらの素子間距離が d1 , d2 , d3 , d4 , d5 の 5 つある．したがって，4 × 3 素子の八木・宇田アレーアンテナの設計は，これら 11 変数の最適化問題と捉えることができる．このような多変数のアンテナ構造の最適化問題は，アンテナ構造の候補が変数に対して指数関数的に増加するのに加え，アンテナ性能を求めるのに Full-wave の数値解析を何度も繰り返し解く必要があるため，現実的に解くのが非常に困難である．したがって，このような最適化問題を解くには何らかの高速化が必要である．本報告では，最適化アルゴリズムは用いずに総当たりの数値計算を行う一方，
電磁界数値解析はベクトル型スーパーコンピュータ SX-ACE で高速化して実行する．
本アンテナは遠方界を用いた無線電力伝送用のアンテナであるため，最適化のための目的関数は指向性利得とし，その最大化を行うこととする．八木・宇田アンテナ素子の導波器，放射器，反射器の素子長とそれらの素子間隔は，離散的な値の候補から選ぶこととし，素子長は 0.02, 0.03, 0.04, 0.05, 0.06

m の 5 通り，素子間隔は 0.01, 0.02, 0.03 m の 3 通りとした．一般的な八木・宇田アンテナの素子長は半波長前後，素子間隔は 1/4 波長前後なので，このような素子長と素子間隔の候補の与え方は妥当である．その一方で，56 ∗ 35 = 3, 796, 875 とおよそ 400 万通りに上る素子長と素子間隔の組

— 16 —

SENAC Vol. 51, No. 1（2018. 1）

20Eθ on E-plane (yz-plane),
f=2.45 GHz, a=0.001 m,
l1=0.06 m, l4=0.04 m,
l2=l3=l5=l6=0.05 m,
10 d =d =d =d =d =0.03 m
1
2
3
4
5
dx=0.067 m, dy=0.067 m

Directivity [dBi]

My code
FEKO

0

-10

-20
180

120

60

0

θ [deg.]

60

120

180

図 2: 最適化した 4 × 3 素子の八木・宇田アレーアンテナの指向性利得.
み合わせに対応したアレーアンテナを高速に数値解析するため，Richmond のモーメント法 [11],[12] の自己・相互インピーダンス行列の数値計算部分をベクトル化したのに加え，アンテナの電流分布を求める部分にベクトル演算向きの共役勾配法を用いた [13]-[15]．素子長と素子間隔の値は逐次的に更新し，指向性利得が最大となるような素子長・素子間隔を最終的に出力する．

3

数値計算結果最適化した 4 × 3 素子の八木・宇田アレーアンテナの指向性利得を図 2 に示す．数値計算結果の妥当性を示すため，商用の電磁界シミュレータである FEKO で同じモデルを数値解析して得られた指向性利得も示す．指向性利得の最大値は約 17.6 dBi(FEKO)，18.3 dBi(My code) であった．
指向性利得から計算した実効開口効率は約 194%であり，効率が 100%を超えていることから，八木・宇田アンテナを側面に配置して開口を立方体の外にまで広げた効果が見られる．各素子の動作インピーダンスは，実部が 29〜40，虚部が-14〜-18 程度の範囲内に収まった．各素子の動作イ
ンピーダンスがこの範囲に収まっていれば，整流回路を適切に設計することでアンテナ-整流回路間の整合を保つことが可能である．
ベクトル型スーパーコンピュータ SX-ACE を用いて数値計算を実行するのに要した時間は，1ベクトル型スーパーコンピュータを用いた八木・宇田アンテナアレーの高利得化
SCHWARZBECKUSL9143, Pt=30 W,
Gt=5.23 dBi@2.45 GHz— 17 —


(+)0.20.20.2 m3
Gr=18.3 dBi@2.45 GHzr

図 3: 遠方界無線電力伝送システムの一例.
ノード (4 並列) で合計 174,897 秒，ベクトル演算率は 99.3%であった．数値計算 1 回当たりの時間
は 0.04 秒ほどと非常に短いものの，400 万回近い電磁界数値解析を実行したため，設計にはおよそ 50 時間を要した．市販のパソコンだとこの数十倍〜数百倍の時間がかかってしまうので，多変数のアンテナの最適化問題において，ベクトル型スーパーコンピュータは強力な高速化ツールに成り得ることが分かる．

4

受信電力
無線電力伝送システムには，受信側にアンテナと整流回路からなるレクテナを配置する．レク

テナは，アンテナで受信した高周波電力を整流回路で直流電力に変換することができる．整流回
路は非線形素子であるダイオードを含むので，整流回路の整流効率は，周波数だけでなくその入力電力の関数にもなる．したがって，高い整流効率を実現するためには，レクテナが受信する電力を予め見積もっておく必要がある．
そこで実際の遠方界無線電力伝送システムを想定し，設計したアレーアンテナの各素子の受信電力を求める．図 3に，遠方界無線電力伝送システムの一例を示す．この図に示す系は，2017 年

9月に電子情報通信学会の無線電力伝送研究会が主催した無線電力伝送コンテストで実際に用いられたものである．送信アンテナは SCHWARZBECK USL9143 であり，受信アンテナは本報告で設計した 4 × 3 素子の八木・宇田アレーアンテナである．送受信アンテナ間距離は 1m, 送信電
力は 30W である．
フリスの伝達公式から，送受電アンテナ間での不整合を無視すれば，受信電力 Pr は以下の式で
表される．

Pr =

(

λ
4πr

)2

Gr Gt Pt(1)

— 18 —

SENAC Vol. 51, No. 1（2018. 1）

ここで，λ は波長，r は送受信アンテナ間距離，Gr は受信アンテナの絶対利得， Gt は送信アン
テナの絶対利得，Pt は送信電力である．(1) 式に Gr = (101.83 )/12 = 5.63, Gt = 100.523 = 3.33,

Pt = 30, λ = 0.122, r = 1 を代入すると，Pr ≈ 53 mW を得る．この値が 1 素子当たりの受信電力であり，この入力電力を目安として整流回路の設計をすればよいことが分かった．

5

まとめ
本報告では, 遠方界を用いた無線電力伝送システム用の受信アンテナの高速設計法を提案した．

提案手法は，素子間相互結合を厳密に評価しながらアレーアンテナを設計する手法であり，反復法とベクトル型スーパーコンピュータを併用することによって高速化を実現している．一定の容積に収まるという条件下で，八木・宇田アンテナ素子からなる 12 素子アレーの設計を行った結果，200%近い開口効率が得られた．遠方界無線電力伝送システムを構築し，素子当たりの受信電力を見積もり，整流回路設計のための指針を得た．

謝辞
本研究成果は, 東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-ACE を用いて得られたものである. スタッフをはじめ，関係各位に深く感謝する．

参考文献[1] R. Reinhold, “Der empfang elektrischer wellen in der drahtlosen telegraphie”(“The receipt
of electric waves in the wireless telegraphy”), Annalen der Physik IV, 25, 1908, pp. 446-466.
[2] N. Tesla, “The transmission of electrical energy without wires,” Electrical World and Engineer, March 5, 1904.
[3] N. Tesla, US Patent No. 1,119,732, 1914.
[4] A. Kurs, A. Karakis, R. Moﬀatt, J. D. Joannopoulos, P. Fisher, and M. Soljacic,
power transfer via strongly coupled magnetic resonances,

Wireless

Science, vol. 317, no. 5834, pp.

83-86, July 2007.
[5] P.E. Glaser, “Power from the sun: Its future,” Science, vol.162, pp.857-861, Nov. 1968.
[6] P.E. Glaser, “An verview of the solar power satellite option,” IEEE Trans. Microw. Theory
Tech., vol.40, no.6, pp.1230-1238, June 1992.
[7] J. O. McSpadden and J. C. Mankins, “Space solar power programs and microwave wirelesspower transmission technology,” IEEE Microw. Mag., vol. 3, no. 4, pp. 46-57, Dec. 2002.

ベクトル型スーパーコンピュータを用いた八木・宇田アンテナアレーの高利得化— 19 —[8] A.K.M. Baki, K. Hashimoto, N. Shinohara, T. Mitani, and H. Matsumoto, “Isoscelestrapezoidal-distribution edge tapered array antenna with unequal element spacing for solar
power satellite,” IEICE Trans. Commun., vol.E91-B, no.2, pp-527-535, Feb. 2008.
[9] K. Konno, Q. Chen, K. Sawaya, and T. Sezai, “Analysis of huge-scale periodic array antennausing impedance extension method,” IEICE Trans. Commun., vol.E92-B, no.12, pp.38693874, Dec. 2009.
[10] K. Konno, Q. Chen, K. Sawaya, and T. Sezai, “Statistical analysis of huge-scale periodic
array antenna including randomly distributed faulty elements,” IEICE Trans. Commun.,vol.E94-C, no.10, pp.1611-1617, Oct. 2011.
[11] J. H. Richmond and N. H. Geary, “Mutual impedance of nonplanar-skew sinusoidal dipoles,”IEEE Trans. Antennas Propag., vol. 23, no. 3, pp. 412-414, May 1975.
[12] C. W. Chuang, J. H. Richmond, N. Wang and P. H. Pathak, “New expressions for mutual
impedance of nonplanar-skew sinusoidal monopoles,” IEEE Trans. Antennas Propag., vol.
38, no. 2, pp. 275-276, Feb. 1990.
[13] M.R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,”
J. Res. Nat. Bur. Stand., vol.49, no.6, pp.409-436, Dec. 1952.
[14] T.K. Sarkar, K.R. Siarkiewicz, and R.F. Stratton, “Survey of numerical methods for solution of large systems of linear equations for electromagnetic ﬁeld problems,” IEEE Trans.
Antennas Propag., vol.AP-29, no.6, pp.847-856, Nov. 1981.
[15] T.K. Sarkar and S.M. Rao, “The application of the conjugate gradient method for thesolution of electromagnetic scattering from arbitrarily oriented wire antennas,” IEEE Trans.
Antennas Propag., vol.AP-32, no.4, pp.398-403, April 1984.

— 20 —

SENAC Vol. 51, No. 1（2018. 1）［共同研究成果］

乱流の生成とその維持機構に関する研究
伊澤 精一郎：東北大学大学院工学研究科
吉川 穣：宮城県産業技術総合センター李 根燮
：岐阜大学大学院工学研究科
西尾 悠：東北大学大学院工学研究科
福西 祐：東北大学大学院工学研究科

1.

はじめに

私たちは様々な乱流に囲まれて暮らしている．空気や水の流れそのものは直接には見えないが，木枯らしに揺れる木々や工場の煙突からもくもくとはき出される煙，海岸に繰り返し打ち寄せては砕け散る波頭を通して複雑で込み入った流れを日常的に目にしている．また屋内においても，洗濯機やぐつぐつと沸騰した鍋の中で洗濯物やうどんがくるくると踊るさまを通して，激しく波打ち複雑に渦を巻いている流れを知ることができる．このような複雑で不規則な乱流は，激しい変動がもたらす強い混合作用を有するため工学的にも大いに利用されており，冷却や内燃機関における燃焼効率の向上などの直接的な寄与ばかりでなく，自動車の排気ガスや工場の煤煙といった汚染物質の大気拡散を促進する間接的な側面もある．その一方で，管内流や流体中を物体が移動する場合には摩擦抵抗の大幅な増大をもたらす．特に大型の旅客機ともなると，摩擦抵抗の割
合は全抵抗の約半分にも達してしまう[1]．摩擦抵抗を減らし燃費の大幅な改善を図るには，旅客機の後退翼面上に発達する 3 次元境界層の遷移過程を把握する必要があり，その物理を知らなけ
れば空力性能に優れた翼型を設計することは難しい[2]．また，管内を流れる乱流を抑制して圧力損失を回復させようとするのであれば，アクチュエータを駆動させるなどして乱流のエネルギー伝達機構に働きかけ，乱流の制御維持システムを破壊する必要がある．流れを制御するための手法はこれまでにも様々なものが提案されている[3]が，その多くは吹き降ろしのある個所では壁面から噴流を噴射し，逆に吹き上げがある個所では吸い込みをかけるといった対処療法的な方策である．このことは，遷移過程を含む乱流についての知識が未だ十分ではないことを示唆している．
一般に，層流から乱流への遷移の初期段階は線形安定論による流れの安定性問題として扱われ，基本流に対する固有モードの撹乱を求める固有値問題に帰着される[4]．一方，発達した乱流の統計的な性質は，Kolmogorov らの乱流理論によって説明される[5]．しかし，これらの数学的な理論が対象とするのは，振幅の非常に小さい撹乱でありまた Reynolds 数（Re 数）の非常に大きい極限の状態であって，乱流の生い立ちを含む物理的な機構の詳細は未だ詳らかではない．近年のコンピュータ技術の目覚ましい発展とともに，これまで統計的にしか扱うことのできなかった乱流運動を直接解析することが可能になりつつあるが，現在までのところその平均的な描像を得るのが精一杯で詳細な解析まで行うには至っていない．
このような背景のもと，3 次元境界層の不安定化機構の解明と乱流の初生に関わる現象の捕捉という 2 つのテーマについて，東北大学サイバーサイエンスセンターとの共同研究課題として実施する機会を得た．本稿では，その内容について報告する．

2.

3 次元境界層の不安定化機構

2.1 研究背景と問題設定
一般に，音速付近で飛行するジェット機の主翼には，翼に流入する気流速度を低下させ臨界マッハ数に達しにくくする目的で後退翼が採用されている．しかし後退翼においては，外部流の方

— 21 —

乱流の生成とその維持機構に関する研究

向と翼面上の圧力勾配の方向が異なるため，翼の前縁及び後縁付近で横流れが生じて境界層が 3次元化し，乱流遷移が早まるなどの欠点もある．3 次元境界層を対象とした研究はこれまでにも理論や計算，実験など様々な方法で行われてきたが，Re 数の他にも翼型や後退角などパラメータが多いことから，しばしば形状の単純な後退円柱で代用されてきた．特に，一定の角速度で回転する円板上の流れは，流れの方向と最大圧力勾配の方向が異なることから航空機の後退翼と同様に 3 次元化して本質的に同じ遷移過程をたどること，またその速度分布が Navier-Stokes 方程式の厳密解であり相似解となることから，3 次元境界層の代表例として研究が行われてきた．
一定角速度Ωで回転する円板上の速度分布は，
����Ω�� ��� =� �� ��� =
��
��
W(z)
rUθ (z)と無次元化でき，層流境界層厚さ
�
z
�=�
�
rUr(z)が半径方向位置 r によらない自己相似な速度分
布（von Kármán flow）となる．また，円板表面
θからの無次元高さはr� �� = �� =� �

図 1 回転円板上の速度分布で与えられ，壁面垂直方向の速度成分と圧力は������ =
� ���� =���√��と一定の値をとる．ここでρは流体の密度，νは動粘性係数である．このとき，� �
�� = �� =� �

で定義される Re 数は，回転円板中心からの無次元距離に一致する．
これまでこのような回転円板上に形成される流れ場を対象として，様々な実験や理論解析が行われ，3 次元境界層の安定性について調べられてきた．例えば，流れ場の可視化や熱線流速計による計測結果から，遷移領域では横流れ不安定性により周方向に並んだ 28～32 本程度の定在的な縦渦が現れることが知られている[6]．これは，平板表面のわずかな凹凸が流れ場を汚染し，その結果として臨界 Re 数の低い移流不安定性が支配的となって境界層内に含まれていた速度変動が増幅されためである．回転円板上の流れは，このような移流不安定性のほかに全体不安定性も存在することがわかっているが，実験では両者を切り分けて議論することができない．このため，
表面粗さの影響を受けない理論解析や数値計算によって，臨界 Re 数や増幅率，不安定波数を調べる試みが数多くなされてきた．古典的な線形安定性解析では，しばしば並行流近似が用いられ，
流れの局所的な不安定性を調べることで系全体の不安定性（全体不安定性）を推測することが行われる．回転円板流れも自己相似な速度分布を持つことから，並行流近似のもとで流れ場の安定性が調べられてきた[7]が，その妥当性に疑義があること[8][9]に加え，予測された局所的な絶対不安定性と系の全体不安定性の関係についても断片的にしかわかっていない[10][11][12]．一方数値計算では，計算負荷抑制の観点から，周方向に周期境界条件を課すことで扇形の狭い範囲に計算領域を縮小して設定せざるを得ないのが現状である．解像可能な波数が全周分割数の整数倍に限定されるため，周方向の不安定波数（実験で観測された縦渦の本数に対応）の正確な値を知ることができない，異なる波数間の非線形干渉が限定的にしか扱えないなどの問題があった．
そこで本研究では，実際の流れ場にできるだけ近い条件のもとで全体不安定性とその臨界波数について調べるため，東北大学のスーパーコンピュータ SX-ACE を用いて回転円板流れの全周計算を試みた．

— 22 —

SENAC Vol. 51, No. 1（2018. 1）

2.2 計算方法
支配方程式である 3 次元 Navier-Stokes 方程式と連続の式
1 ���� �� � �� � �� � � � �� � ��� �� � � ��� � � � � ������を差分法で解いた．ここで，U = (U,V,W)は基本流の半径方向，周方向，壁面垂直方向の速度成分であり，u = (u,v,w)は各軸方向の速度変動成分である．外力 f は，流入及び流出境界における速度変動の非物理的な反射を抑制するために設けた緩衝領域で作用する減衰力である．圧力場と速度場のカップリングには MAC 法を用い，空間微分項の離散化には式 3 次精度上流差分(Kawamura- Kuwahara スキーム) を，その他の項には 4 次精度中心差分を用いた．時間進行には
2 次精度 Crank-Nicolson 法を用いた．計算領域として Re 数が 220 から 752 までの領域を確保し，流入（Re = 220～249）及び流出境界（Re = 621～752）には速度変動を滑らかに減衰させるための緩衝領域を設けた．計算方法の詳細は既報[13]を参照されたい．
流れ場に導入する撹乱は，全体不安定性が発現するとされる臨界 Re 数 583[14]近傍 570～590の領域において，0.016T（T は回転周期）の間だけ，正味の質量流束が零で最大振幅がその位置の周方向速度に等しくなるような壁面垂直方向速度をランダムに与えた．攪乱を与えた瞬間を時
刻 t = 0 とした．回転円板の回転数は 1rad/s である．計算格子は壁面近傍に集中させた不等間隔格
子とし，格子点数は 533×4,609×63 点とした．計算は SX-ACE の 64 ノードを使用して行い，円板
1 回転あたりおよそ 2 週間の計算時間を要した．計算は約 7.7 回転分行ったので，計算終了までにはおよそ 1 か月半かかったことになる．
2.3 全体不安定性と縦渦構造
図 2 は，撹乱を導入してから 7 回転後の流れ場の様子である．図中の色は，回転円板表面から
z = 1.3 の高さにおける周方向速度変動成分の大きさの対数値 log|v|を表しており，図は緩衝領域を除いた Re = 250～620 までの有効計算領域のみ表示してある．外周寄りの領域には，ほぼ一定の間隔で円板の回転と同期した定在的ならせん状のパターンが見られることから，本計算でも実験で観測された横流れ不安定性による縦渦の存在を確認することができる．なお，ここでは周方向の速度変動成分の絶対値|v|を対数表示しているので，らせん状のパターン 2 本が縦渦 1 本に対応
していることに注意されたい．また，図中のらせんの本数を数えてみると全周で 162 本あることから，都合 81 本の縦渦が形成されていることがわかる．Imayama ら[15]は，本計算と同じ z = 1.3の高さにおいて熱線流速計による詳細な計測を行い，これら定在渦の振幅は周方向に微妙に異なるものの，乱流に遷移する位置は個々の渦の強さではなく Re 数によって規定されることを報告した．図 2 の結果をよく見ると，本計算においても速度変動の強さ，すなわち縦渦の強さが周方向に微妙に変化しており，大きさもまちまちな 15 の小領域単位で変化していることも見て取ることができる．このことは速度変動の周方向の揺らぎが円板の表面粗さとは無関係に起こりうることを示唆しており，波長が一定の変動しか想定することができない安定論や，周期境界を用いる
ことで周方向に強い周期性が要請される従来の数値計算では原理的に捉えることができず，全周計算を行うことで今回初めて明らかになったことのひとつである．
周方向速度変動成分の乱れ強さ vrms を図 3 に示す．Re = 580 付近から導入されたランダムな撹乱は下流へと流下するが，全体不安定領域に達すると増幅に転じ，7 回転後にはほぼ定常状態に達した．本研究では全体不安定性による撹乱の増幅過程に主眼を置いているため，有効計算領域の下流端は完全発達した乱流が現れる Re = 650 よりも上流の Re = 620 の位置に設定している．したがって，定常状態に達していても計算領域内で乱流遷移が完了することはない．なお，回転初期の t = 0.5～1.0 において，Re 数が 550 までの領域で一時的に導入した撹乱とは別の孤立撹乱が増幅しながら流下する様子が見られたが，これは壁面から短時間噴流を噴射した際の圧力変動が上流境界で反射し，緩衝領域で十分に減衰しきれなかったものがノイズとして計算領域にしみだしてしまったためである．その乱れ強さは全体不安定性で増幅する撹乱のそれと比べて小さく，
結果には影響しない．

乱流の生成とその維持機構に関する研究— 23 —

図 4 は，t = 7，Re = 600 における周方向速度変動成分の周波数解析結果である．図から波数 81付近に最大のピークがあり，またその高調波成分にもピークがあることがわかる．Lingwood [7]は，並行流近似にもとづいた線形安定性解析の結果から群速度が零となる位置を絶対不安定性の臨界 Re 数（= 507）として求め，その臨界波数が 68 であることを報告している．最も不安定となる波数は，むろん Re 数に依存した値である．しかし ，撹乱の初期振幅が同じであれば中立安定
点の Re 数が低いものほど早く成長することになるため， 先行研究の多くは線形安定論の枠組みの中で得られたこの臨界波数 68 に着目してきた．これに対して，回転円板全周を計算領域とした
本研究の結果は，従来考えられてきた波数よりも高い波数帯域に全体不安定性により増幅されやすい波数が存在していることを示しており，現象の解明には大規模計算が必要不可欠であることが改めて認識される結果となった．

図 2 周方向速度変動成分で見る流動パターン（z = 1.3）

— 24 —

SENAC Vol. 51, No. 1（2018. 1）

図 3 導入した撹乱の時空間成長（log(vrms), z = 1.3）

6

10

4

10

2

Power

10

0

10

−2

10

−4

10

−6

10

0

200

400

600

800

1000

Azimuthal wavenumber

図 4 周方向速度変動成分のパワースペクトル（z = 1.3）

1200

乱流の生成とその維持機構に関する研究3.

— 25 —

乱流の初生に関わる現象の補足

3.1 研究背景と問題設定
境界層の乱流遷移は，境界層内部に含まれる乱れが増幅され互いに非線形干渉を起こした結果である．その遷移のきっかけとなる乱れの種がどこから来るのかということは，長い間研究者たちの論争の的であった．これは，実際に境界層内で観測される不安定波動の波長や位相速度といったスケールが，乱れの元と考えられた主流に含まれる速度変動や音波などの外乱のそれと必ずしも一致しないためである．この問題は受容性（receptivity）の問題として知られ，理論解析や風洞実験により様々に研究が行われ．多くの知見が蓄積されてきた[16]．現在では，翼の前縁のように速度や圧力の空間変化が激しく流れの非平行性の無視できない領域で波長のマッチングが行
われ，境界層中に乱れの種として取り込まれることが研究者の共通した認識となっている．さらに，壁面曲率が不連続となる箇所でも受容性が高くなること，受容される乱れの振幅は外乱の大きさに比例し，その後の遷移の道筋は初期振幅の大小によって粘性型不安定波動（Tollmien-Schlichting 波，T-S 波）の線形増幅に始まるタイプと撹乱の過渡増幅にともなうストリ
ーク不安定タイプ（バイパス遷移）に大別されること，そしてそれぞれが全く異なる遷移形態をとることもよく知られた事実である．前者の粘性型不安定波動の増幅過程（modal growth）については，すでに述べたように基本流に対する固有モードの撹乱を求める固有値問題に帰着される．
非線形性が極めて強い後者の過程）については，理論的には固有関数の非直交性による過渡増幅機構（non-modal growth）として理解され，縦渦の生成とストリーク構造の形成及び不安定化を経て乱流に遷移する．ストリーク構造は流れ方向に低速領域と高速領域が交互に並んだ縞状の配置をとり，バイパス遷移過程だけでなく乱流境界層の壁面近くに現れる特徴的な流体構造である．
ストリーク構造には，Varicose モード（対称）と Sinuous モード（反対称）と呼ばれる不安定性があるが，通常は Sinuous 不安定性が支配的となってストリークの崩壊が起こり，層流境界層中に孤立した乱流領域（乱流斑点）が次々に生成され乱流遷移が進展する[17][18]．このように乱れの受容から乱流遷移に至る過程について多くの知見が報告されているが，いずれの遷移過程をたどって乱流に遷移するにせよ乱流の誕生は非線形作用が卓越した現象であり，その詳細については詳らかでないことが多い．
そこで本課題では，平板層流境界層中に励起した低速ストリークに対して壁面から短時間噴流を噴射し，下流で乱流に遷移する場合としない場合の流れ場の様子を詳細に比較することで，乱流の初生に関わる構造の特定を試みた．計算は，東北大学の並列コンピュータ LX 406Re-2 を用
い，噴射条件を様々に変えて計算を行った．
3.2 計算方法
支配方程式は 3 次元 Navier-Stokes 方程式と連続の式であり，速度や圧力などの変数は一様流速度及び流入部層流境界層の排除厚さを用いて無次元化して計算を行った．差分法の計算アルゴリ
ズムには MAC 法を用い，時間進行には 2 次精度 Crank-Nicolson 半陰解法，空間微分項の離散化には移流項にのみ 3 次精度上流差分法である Kawamura-Kuwahara スキームを，その他の項には 2次精度中心差分法を用いた．
計算は平板層流境界層中にストリークを励起するための予備計算と，生成されたストリークに対して噴流を噴射して境界層の不安定化を促す主計算に分けて行った．いずれの場合においても，計算領域は線形安定論により予測される不安定領域内に設定し，撹乱が増幅されやすい環境を与えた．座標系の原点は流入部スパン方向中央の壁面上とし，流れ方向を x，壁面垂直方向を y，ス
パン方向を z とした．予備計算の領域の大きさは (x, y, z) ∈ ([0, 400], [0, 30], [-20, 20]) であり，流入境界は平板前縁から 177 離れた位置（排除厚さにもとづく Re 数で 530）にあるものとした．また．計算格子には壁面近傍を密にした不等間隔格子を用い，2,380×114×280 の格子点数を用いて計算を行った．1 ケースあたりの演算時間は，Open MP による 24 並列で 11,200 ステップ計算して約 4.5 日であった．

— 26 —

SENAC Vol. 51, No. 1（2018. 1）Andersson ら[19]の報告によれば，スパン方向の波数が 0.45 となるとき過渡増幅が最大となる．
そこで，噴流噴射位置（x = 100）でこの最適波数の間隔と等しくなるように，その上流 x = 50 の位置に一辺が 1×1×2 の直方体を 4 個スパン方向に 10 おきに並べ，層流境界層中にストリーク構造を励起した．主計算の計算領域はその下流 x = 60 からとし，予備計算の結果を初期速度分布及び流入境界条件として与えた．噴流は 2×2 の小孔から壁面垂直方向に無次元時間 15 の間だけ一様噴射させるものとし，その噴射速度 vjet を変えながら乱流遷移が起こる下限の値を探った．以下では，低速ストリークの中心で噴流を噴射した結果について紹介する．なお，噴射速度が最大の場合であっても，噴射位置と流入境界は十分に離れていることは確認済みである．

3.3 乱流のはじまりに関わる渦構造
物体列から 50 だけ下流 x = 100 の yz 断面における，予備計算で得られた基本流となるストリー
クを含む速度分布を図 5 に示す．物体下流 z = 0 付近には周囲より速度が低下した領域（低速ストリーク）が見られ，またその両側 z = ±1～2 にかけて周囲より速度が速い領域（高速ストリーク）が見られることから，ストリーク構造が形成されていることが確認できる．
図 6 は，この低速ストリークに対して，壁面から一様流速度の 13%及び 16%の強さで噴流を噴射した場合の流れ場の様子を比べたものである．渦構造は速度勾配テンソルの第 2 不変量である
Q 値の等値面で可視化し，局所的な流れ方向の渦度成分の強弱に応じてその表面を着色した．低速ストリークに対して噴流を噴射した直後の t = 30 の時刻では，噴流の強弱によらず噴射孔周囲に首飾り状の渦が，またその下流側には Varicose モードが刺激されて形成されたヘアピン状の渦が見られ，互いによく似た構造が生成されている．どちらのケースにおいても，時間が進むにつれてその上流側に新たなヘアピン渦が生成される様子が見られたが，その後の様子は全く異なったものとなった．噴射速度が小さい 13%のケースでは，生成されたヘアピン渦は粘性により散逸
しながら流れ去ってしまい，後には流れ方向に伸びた縦渦が残されるばかりで乱流に遷移することはなかった．
これに対して噴射速度を一様流の 16%まで大きくすると，初期に形成されるヘアピン渦の位置が高くなり，境界層の速度勾配によりすぐに伸長されて次々に新しいヘアピン渦が誕生する様子

図 5 流れ方向速度分布（x = 100）

— 27 —

乱流の生成とその維持機構に関する研究が見られた．これらのヘアピン渦は，13%のケースと同様下流に流れ去ってしまい，乱流遷移の直接の引き金とはならなかった．しかし，その過程で時刻 t = 200 で見られるような流れ方向に長く伸びた複数の縦渦対が生き残り，初期に形成され壁面近くの低速領域に取り残されていた上流側の縦渦群と干渉し，さらに複雑な構造を生み出す原動力となることがわかった．そこで時刻 t =200 前後の渦構造の変化を細かく見てみると，x = 205 付近に見られる流れ方向に傾いたハの字状の横渦対がその内側の縦渦対とつなぎ替えを起こし，新たにヘアピン渦を形成しながら複雑化が進展することから，この横渦対が乱流のはじまりに関わる鍵となる構造の 1 つであることがわかった．時刻 t = 350 の流れ場を見ると，複雑に絡み合った無数の渦構造が見られ，この時刻では境界層はすでに乱流状態にあるものと考えられる．詳細については，現在解析を行っているところである．

(a) vjet =0.13(b)
図 6 次項に続く

vjet = 0.16

— 28 —

SENAC Vol. 51, No. 1（2018. 1）(a) vjet =0.13

図6

(b)

vjet = 0.16

Q = 0.0015 の等値面で可視した渦構造の比較：t = 30 (1), 100 (2), 200 (3), 275 (4), 350 (5)（局所縦渦成分ωx が正の領域を赤，負の領域を青色付けして表示）

乱流の生成とその維持機構に関する研究— 29 —

謝辞
本研究は、東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．これらの研究成果の一部はすでに学術誌に投稿中であり，また残りの成果についても現在投稿準備を進めているところである．ここに改めて，同センターの関係各位に感謝申し上げたい．
参考文献[1] 中橋 和博: 航空機の空力形状と最適設計，日本流体力学会誌「ながれ」，26 巻 4 号，pp.259‐
265，2007．
[2] 小濱 泰昭, Saric William S. and Hoos Jon A.: 後退翼における境界層の乱流遷移，日本機械学会
論文集 B 編，58 巻 554 号，pp.3053‐3059, 1992．
[3] L.N. Cattafesta III and M. Sheplak: Actuators for Active Flow Control, Annual Review of Fluid
Mechanics, 43, pp. 242-272, 2011.
[4] P.J. Schmid and D.S. Henningson: Stability and Transition in Shear Flows, Springer Science &Business Media, 2012.
[5] P.A. Davidson: Turbulence: An Introduction for Scientists and Engineers, Oxford University Press,
2004.
[6] Y. Kohama: Study on boundary layer transition of a rotating disk, Acta Mechanica, 50, pp.193-199,
1984.
[7] R.J. Lingwood: Absolute instability of the Ekman layer and related rotating flows, Journal of Fluid
Mechanics, 331, pp.405-428, 1997.
[8] N. Itoh: Structure of absolute instability in 3-d boundary layers: Part 1. Mathematical formulation,
Transactions of the Japan Society for Aeronautical and Space Sciences, 44, pp.96-100, 2001a.
[9] N. Itoh: Structure of absolute instability in 3-d boundary layers: Part 2. Application torotating-disk-flow, Transactions of the Japan Society for Aeronautical and Space Sciences. 44, pp.101-105,
2001b.
[10] C. Davies and P.W. Carpenter: Global behaviour corresponding to the absolute instability of therotating-disk boundary layer, Journal of Fluid Mechanics, 486, pp. 287-329, 2003.
[11] B. Pier: Finite-amplitude cross flow vortices, secondary instability and transition in the rotating-diskboundary layer, Journal of Fluid Mechanics, 487, pp. 315-343, 2003.
[12] J.J. Healey: Model for unstable global modes in the rotating-disk boundary layer, Journal of Fluid
Mechanics, 663, pp.148-159, 2010.
[13] K. Lee, Y. Nishio, S. Izawa, and Y. Fukunishi: Effects of location of excitation on the spiral vortices inthe transitional region of a rotating-disk flow, Physics of Fluids, 29, 084106, 2017.
[14] E. Appelquist, P. Schlatter, P. H. Alfredsson, and R. J. Lingwood: On the global nonlinear instabilityof the rotating-disk flow over a finite domain, Journal of Fluid Mechanics, 803, pp.332-355, 2016.
[15] S. Imayama, P.H. Alfredsson and R.J. Lingwood: On the laminar-turbulent transition of the rotating
disk flow: the role of absolute instability, Journal of Fluid Mechanics, 745, pp.132-163, 2014.
[16] Saric, W., Helen L.R. and Edward J.K.: Boundary layer receptivity to freestream disturbances, Annual
Review of Fluid Mechanics, 34, pp.291-319, 2002.
[17] M. Asai, Y. Konishi, Y. Oizumi and M. Nishioka: Growth and breakdown of low-speed streaks leading
to wall turbulence, Journal of Fluid Mechanics, 586, pp.3791-396, 2007.
[18] P. Schlatter, L. Brandt, H.C. de Lange and D.S. Henningson: On streak breakdown in bypass transition,
Physics of Fluids, 20, 101505, 2008.
[19] P. Andersson, M. Berggren, and D.S. Henningson: Optimal disturbances and bypass transition inboundary layers, Physics of Fluids, 11, pp.134-150, 1999.

— 30 —

SENAC Vol. 51, No. 1（2018. 1）［共同研究成果］

スーパーコンピュータによるリアルタイム津波浸水被害予測
越村俊一・阿部孝志：東北大学災害科学国際研究所
撫佐昭裕：東北大学サイバーサイエンスセンター、NEC
村嶋陽一・鈴木崇之・井上拓也：東北大学災害科学国際研究所、国際航業株式会社
太田雄策・日野亮太：東北大学大学院理学研究科
佐藤佳彦：NEC ソリューションイノベータ株式会社
加地正明：株式会社エイツー
小林広明：東北大学大学院情報科学研究科
巨大津波発生直後のきめ細かな津波情報や、迅速な被害情報の推計・把握と配信を通じて被災地を支援し、災害に対するレジリエンス（回復力）の向上と我が国の国土強靱化に資する、リア
ルタイム津波浸水被害予測システムを構築した。スーパーコンピュータ SX-ACE の活用により、沿岸部 10ｍメッシュ分解能でのリアルタイム津波浸水予測（浸水する範囲と深さの予測）を可能とし、浸水域内人口および建物被害予測を、地震発生から 20 分以内を目安に完了することに成功した。本システムは、内閣府の津波浸水被害予測システムとして採用され、2017 年 11 月より運用を開始した。

1.

はじめに

巨大災害後の災害対応や救援活動において最も重要なことの一つは被害の全容把握である。地震や津波災害といった巨大災害の発生直後は、激甚な被害を受けた地域からの情報が断片的となり、被害全容の把握がきわめて困難になるとともに、被災地の救援活動や復旧活動も難航する。
2011 年東北地方太平洋沖地震津波の被災地は広大であり、発災直後には、激甚な被災地がどこにあるかを把握することさえ困難であると同時に、現地調査期間や人的資源の制約により被害全容を把握するにはきわめて長い時間を要した[1]。この教訓を踏まえ、効果的な災害対応、被災地の救援のためには、津波高さの予測だけでなく浸水範囲・被害の迅速かつきめ細かな予測が必要である。
本研究は、最先端のシミュレーション・センシング・ICT を統合し、津波発生直後のきめ細かな津波情報や、迅速な被害情報の推計・把握と配信を通じて被災地を支援し、災害に対するレジリエンス(回復力)の向上とわが国の国土強靱化に資する、リアルタイム津波浸水被害予測システ
ムを構築することを目標とした。本稿では、災害科学国際研究所、サイバーサイエンスセンター、理学研究科との共同および産学連携により、この課題の解決にむけて構築したリアルタイム津波
浸水被害予測システムについて報告する。

2.

リアルタイム津波浸水被害予測

津波の伝播・陸上遡上の予測には、その領域（沖合・沿岸および浅海域の伝播、陸上での遡上）と分解能（空間・時間）に応じて数値モデルを使い分ける必要がある。いずれにせよ、支配方程式を差分法により離散化する方法が一般的である。長波理論の差分法に基づく津波数値計算の場
合、 座標系と支配方程式を、再現する津波の対象（遠地/近地、 外洋伝播／遡上）に応じて適切に選択する必要がある。2011 年東北地方太平洋沖地震津波の場合、日本近海の津波の再現には直交座標系による非線形長波理論式を、太平洋全体への外洋伝播を含めた津波を再現する場合には球面座標系による分散波理論式および線形長波理論式が必要である。いずれにせよ、正確な津波の予測・再現には、方程式系の適切な選択、津波初期水位分布（断層運動による海底地盤変動）、
詳細な海底・陸上地形の情報と計算の分解能、土地利用状況等による陸上の津波抵抗則の適切なモデル化等が重要な要件となる。

スーパーコンピュータによるリアルタイム津波浸水被害予測— 31 —

津波浸水被害の予測をリアルタイムで行うための課題は 3 つに大別される[2]。
1 点目は、津波の発生予測である。津波数値計算の初期条件には、断層破壊の具合的なメカニ
ズムに関連した断層モデルや津波発生時の初期水位分布の情報が必要で、特に地震学・測地学の研究者との連携が必要になる。近年、GEONET をはじめとする衛星測位（GNSS）技術の発展を背景とした新しい地震・地殻変動観測が普及しており、津波発生モデルの精度向上に期待が持てる。また、特に我が国太平洋側には、日本海溝海底地震津波観測網（S-net）や地震・津波観測監
視システム（DONET）をはじめとして、非常に高密度な津波観測網が整備されつつある。また、これらの観測網の充実を背景として、全く新しい津波発生予測構築手法が提案されており、我が国の津波予測技術の革新が進みつつある。本システムでは、地震情報の取得に、気象庁の緊急地
震速報と GEONET の地殻変動観測情報を用いたリアルタイム断層モデル推定法 RAPiD[3]および
国土地理院の REGARD[4]による解を即時に取得し、津波波源モデルの計算過程を自動化した。
2 点目は、津波の浸水予測である。一般的な非線形長波理論に基づく数値モデルでは従来のワークステーションで計算できるが、リアルタイムでの予測は技術的に困難であったが、High
Performance Computing Infrastructure（HPCI）の普及が課題解決の追い風になっている。例えば、
京コンピュータを利用した大規模高速計算の技術による津波計算の実績も蓄積されている[5]。一方、ほとんどのスパコンの運用体制では災害時の即時的な予測を行うための運用を行うことは難しい。筆者らの研究グループは、津波の予測計算の高速化を、東北大学サイバーサイエンスセン
ターのベクトル型スーパーコンピュータ SX-ACE の独自運用（ディザスターモード：地震発生時に所要の計算リソースを即座にアサインする）により実現し、いつ地震が発生してもスパコンのパフォーマンスを確保している。図 1 のように、10 分以内に津波の発生（断層モデル）を予測、
10m メッシュという高分解能の浸水計算を、10 分以内に完了することを具体的な目標として実証に成功している[6]。また、東北大学の津波解析プログラムは、SX-ACE の高いメモリバンド幅およびベクトル化効率の恩恵により、計算効率が非常に高いことも確認することができた。
上記の成果を踏まえ、本システムは、南海トラフ地震の発生を想定して、鹿児島県から静岡県
までの 6,000km の海岸線を予測対象領域として（図 2）、内閣府の津波浸水被害予測システムの構築・運用を開始した。計算領域が広大であることから、内閣府システムとしては津波被害の概略を迅速に把握することを最優先事項として、沿岸部 30m メッシュの分解能で運用を行う。

図 1 東北大のスパコン SX-ACE による 3 時間分の浸水予測計算（10m メッシュ）のパフォーマンス（高知市で実施）[6]。システムの運用には、東北大・大阪大のスパコンを用いて冗長性を確保する。比較しているのは京コンピュータによるパフォーマンス[5]。

— 32 —

SENAC Vol. 51, No. 1（2018. 1）

図 2 リアルタイム津波浸水被害システムの予測領域。（左：予測領域（90m メッシュまで表示）、
右：浸水予測を行う 30m メッシュの領域構成（高知県の例）。
3 点目は、被害の予測である。津波の浸水域は、湾の構造や建物の密度などによって左右される。量的な被害予測を行うためには、木造建築物と鉄筋コンクリート造などの建物がどれくらいの津波で破壊されるのか量的に解析する必要がある。そのために、東日本大震災での被害調査データから津波被害関数を構築した（図 3）。
図 4 に高知県で実証した津波浸水被害予測結果（南海トラフ想定地震 Mw8.7 の例）を示す。国勢調査が綿密に行われている日本では、人口分布および建物の位置と場所を正確に把握できるため、これらの情報を組み込むことで 10m 区画まで細分化した浸水予測結果から建物被害地図の作成が可能になった。建物の位置と場所を正確に把握し、津波被害関数から得られる流失率を求めることで、細分化した浸水予測結果から建物被害の量的な予測が可能になった。

図 3 2011 年東日本大震災における被害データから得られた津波被害関数の例[1]。浸水域内建物被害（流失）棟数の推定に利用する。

スーパーコンピュータによるリアルタイム津波浸水被害予測— 33 —

図 4 リアルタイム津波浸水被害予測システムの出力結果例（左上：津波の浸水開始時間、右上：
浸水深分布、左下：浸水域内人口分布、右下：建物被害（流失）棟数分布）

3.

おわりに

東日本大震災の教訓を踏まえ、スーパーコンピュータ SX-ACE を利用したリアルタイム津波浸水被害予測システムを構築した。本研究は、本学災害科学国際研究所、サイバーサイエンスセン
ター、理学研究科との共同、大阪大学サイバーメディアセンターおよび産学連携により、取り組
んだ研究成果であり、地震情報自動取得と津波の発生予測（10 分以内）、津波伝播・浸水・被害予測（10 分以内）
、結果の図化・配信を 10m メッシュで行うという目標をたて、現在では、さらに改良を進め、地震発生から 3 時間分の津波浸水予測を 5 分程度で行うことが可能になり、東北大学・大阪大学の SX-ACE で 30 分以内の全国津波浸水予測が可能となった。東北大学の津波浸水予測プログラムを SX-ACE に最適化し、京コンピュータよりも優れた計算性能を達成した。この成果が高く評価され、平成 29 年 3 月に内閣府津波浸水被害推計システムとして採用され、平成
29 年 11 月から運用が開始されることとなった。今後は、日本全体での予測領域の拡大や、様々な主体・ユーザのニーズに対応した予測情報・活用技術の研究を進めていく。

謝辞
本研究は、東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた。また、研究にあたっては同センター関係各位に有益な指導と協力をいただいた。基礎研究に際しては、JST の CREST（グラントナンバーJPMJCR1411）の補助を受けた。ここに記して謝意を表する。
参考文献[1] Koshimura, S. and N. Shuto, Response to the 2011 Great East Japan Earthquake and Tsunami disaster,Philosophical Transactions of the Royal Society A, 2015a. doi:10.1098/rsta.2014.0373[2] Koshimura, S., R. Hino, Y. Ohta, H. Kobayashi, A. Musa, Y. Murashima, Real-time tsunami inundation
forecasting and damage estimation method by fusion of real-time crustal deformation monitoring and— 34 —

SENAC Vol. 51, No. 1（2018. 1）

high-performance computing, IUGG General Assembly, 2015.
[3] Ohta, Y., T. Kobayashi, H. Tsushima, S. Miura, R. Hino, T. Takasu, H. Fujimoto, T. Iinuma, K.
Tachibana, T. Demachi, T. Sato, M. Ohzono, N. Umino, Quasi real-time fault model estimation fornear-field tsunami forecasting based on RTK-GPS analysis: Application to the 2011 Tohoku-Oki
earthquake (Mw 9.0), J. Geophys. Res., 117, B02311, 2012. doi:10.1029/2011JB008750.
[4] Kawamoto, S., Y. Ohta, Y. Hiyama, M. Todoriki, T. Nishimura, T. Furuya, Y. Sato, T. Yahagi, K.
Miyagawa, REGARD: A new GNSS-based real-time finite fault modeling system for GEONET, J.
Geophys. Res., 2017. DOI:10.1002/2016JB013485.
[5] Oishi, Y., F. Imamura, and D. Sugawara. “Near-field tsunami inundation forecast using the parallel
TUNAMI-N2 model: Application to the 2011 Tohoku-Oki earthquake combined with source inversions”.
Geophysical Research Letters, 42, doi:10.1002/2014GL062577, 2015.
[6] Musa, A., H. Matsuoka, O. Watanabe, Y. Murashima, S. Koshimura, R. Hino, Y. Ohta, H. Kobayashi, A
Real-Time Tsunami Inundation Forecast System for Tsunami Disaster Prevention and Mitigation，The
International Conference for High Performance Computing, Networking, Strage and Analysis (SC15),
Austin, Texas, Nov. 2015.

— 35 —

SENAC Vol. 51, No. 1（2018. 1）［高速化支援］

第１原理電子構造計算コード
VASP の第１原理電子構造計算コード VASP のSX-ACE
向け最適化
SX-ACE 向け最適化

山下 毅 †1 , 長谷川 正之 2 , 青山修也 2 , 千葉貴則 2 , and 西館 数芽 ‡2
1
2

岩手大学

東北大学

理工学部

情報部情報基盤課

システム創成工学科電気電子通信コース

概要
パラメータを用いずに物質の電子構造を調べる第１原理電子構造計算法は、物性物理学の分野における非常に強力な研究手法であり、最近では理論研究者のみならず広く実験研究者にも使われるようになってきた。本稿では密度汎関数法に基礎をおく第１原理電子構造計算コード VASP のスーパーコンピュータ SX-ACE 向けの最適化と、並列コンピュータ LX 406Re-2 での実行結果の比較について述べる。

1 はじめに物性理論の精密化とスーパーコンピュータの演算性能の飛躍的向上に伴い、密度汎関数法（Density

Function Theory, DFT）[1] に基礎をおく第１原理電子構造計算法は物質科学の分野において頻繁に使用される研究ツールとなってきた。この方法は物質の物理的・化学的特性を調べるための非常に強力な手段を提供してくれる。通常の意味での経験的パラメータを用いないため、たとえば経験的理論が未だ存在していない特殊な環境下における物質の物理的・化学的特性をも予測できる。現在では、

WIEN2k、CASTEP、VASP[2] などのさまざまな先進的な DFT コードが入手可能である。Vienna
Ab-initio Simulation Package (VASP) は、平面波を基底とする DFT の枠組みにおける第１原理電子構
造計算コードである [3, 4, 5]。原子の内殻は、PAW 法（Projector augmented-wave method）で記述
し、DFT ハミルトニアンの対角化に効率的なアルゴリズムが採用されている。ここでは物質科学の分野でユーザーが非常に多い第１原理計算コード VASP に焦点をあてる。
高性能科学技術計算の追求は、近年の並列計算機の隆盛をもたらした。現在用いられている DFT コードの大部分が、多数の計算ノードから構成される大規模並列計算機での実行を前提にプログラムされている。しかしながらメモリの閾値や大規模な並列化におけるオーバーヘッドのために、単純に計算ノード数を増加させてもそれに比例した速度向上を得ることはできない。この問題に対する一つの解は、各計算ノードに搭載されているプロセッサあたりの性能を高めることである [6, 7]。

SX-ACE は NEC 製のベクトル型スーパーコンピュータである [8, 9]。4 コアを有するベクトルプロ
セッサは、コアの演算性能とメモリバンド幅のバランスの取れたシステムとして設計されている。ただし、その高い演算性能を引き出すためには、コンパイラの最適化を促進させるためのコード修正が必要で† yamacta@tohoku.ac.jp‡ nisidate@iwate-u.ac.jp

— 36 —

SENAC Vol. 51, No. 1（2018. 1）ある。本稿では VASP.5.4.1（vasp.5.4.1.24Jun15.tar.gz）を対象とし、SX-ACE 環境への移植作業および、SX-ACE 向けコンパイラに対応したコード修正、およびコードのベクトル高速化について述べる。
また、FFT ライブラリは付属の Furth ライブラリのソースコードをコンパイルして使用した。

2 SX-ACE 用コンパイラ向け移植作業とベクトル高速化
2.1 SX-ACE 用コンパイラの仕様
SX-ACE 用のコンパイラとして、NEC 製の FORTRAN90/SX（ISO/IEC 1539-1:1997 準拠の Fortran90/95 コンパイラ）、NEC Fortran2003（ISC/IEC 1539-1:2004 準拠の Fortran2003 コンパイラ）、
および、C++/SX（ISO/IEC 1539-1:2004 準拠の C,C++ コンパイラ）が利用可能である。
今回移植の対象とした VASP.5.4.1 は一部 Fortran2003 規格を使用したソースコードがあるが、現時点では FORTRAN90/SX コンパイラの自動ベクトル化/最適化機能が Fortran2003 コンパイラよりも優れているため、本稿では前者のコンパイラの利用を仮定する。Fortran2003 機能を利用するためのオプ
ションを指示することにより、FORTRAN90/SX（rev.534）コンパイラでも VASP.5.4.1 をコンパイルすることが可能である。また、C,C++ コンパイラとして C++/SX（rev.112）コンパイラを使用する。

2.2 FORTRAN90/SX 向けのコード修正
2.2.1 変数定義部の修正
FORTRAN90/SX コンパイラでは、変数の定義部で組込関数を利用できないため、リスト 1 に示す 3
つのソースコードについて、変数の定義部では型宣言だけを行い、実行部で値を代入するように修正を行った。
リスト 1 変数定義部[ src / vdwforcefield . F ]（6047 行 目 ） 定 義 部 で 組 込 関 数 の C M P L X を 使 用 で き な い
COMPLEX ( q ) :: zdumm0 = CMPLX (0. _q ) , zdumm1 = CMPLX (1. _q ) , zdummM1 = CMPLX ( -1. _q )（7824 行 目 ） 定 義 部 で 組 込 関 数 の S Q R T を 使 用 で き な い

REAL ( q ) , PARAMETER :: t wo _ ov er _ sq rt PI =2. _q / SQRT ( PI ) , f o u r _ o v e r _ s q r t P I =2* t w o_ ov er _ sq rt P I[ src / bse_te . F ]（137 行 目 ） 定 義 部 で 組 込 関 数 の R E A L を 使 用 で き な い
REAL ( mq ) , PRIVATE , PARAMETER :: HARTREE = 2. _mq * REAL ( Rytoev , mq )[/ src / dmatrix . F ]（1321 行 目 ） 定 義 部 で 組 込 関 数 の S Q R T を 使 用 で き な い
REAL ( q ) , PARAMETER :: SQPi54 =4. _q * SQRT ( PI /5. _q ) , SQPi52 =2. _q * SQRT ( PI /5. _q ) ,mSQPi52 = -2. _q * SQRT ( PI /5. _q ) , SQPi15 =6. _q * SQRT ( PI /15) ,mSQPi15 = -6. _q * SQRT ( PI /15)
関 連 し て DATA文で定義される配列XIXJ2YLMも実行部で値を代入する

2.2.2 WRITE 文の修正
FORTRAN90/SX コンパイラの仕様に基づき、リスト 2 に示す 2 つのソースコードについて、WRITE
文中で/=を使用する行のコメントアウトと、WRITE 文の直後のカンマを削除した。

第１原理電子構造計算コード VASP の SX-ACE 向け最適化— 37 —

リスト 2 WRITE 文[/ src / gridq . F ]（ 3 3 7 行 目 ） W R I T E 文 中 で /= を 使 用 で き な い
WRITE (0 ,*) GRID % MPLWV /= SIZE (C ,1)[ src / minimax . F ]（182 行 目 ） W R I T E 文 の 直 後 の カ ン マ （ , ） が 文 法 エ ラ ー
IF ( INU >= 0 ) WRITE (* , ’ (" Number of grid points forced to " , I4 ) ’) , N N U

2.2.3 変数の二重定義の修正
変数が二重定義されており、FORTRAN90/SX コンパイラではコンパイル時にエラーとなるため、定義部分の片方をコメントアウトした。
リスト 3 変数の二重定義[ src / base . F ]（169 行 目 か ら190 行 目 ） d i p o l が 2 重 定 義 で エ ラ ー に な る の で 、 定 義 部 分 を コ メ ン ト ア ウ ト

2.2.4 外部関数呼出しの修正
関数 volume を呼び出しているが、他ソースコードで定義されていないためコンパイル時にエラーとなる。実際の利用はないようなので呼び出し部分をコメントアウトした。
リスト 4 外部関数呼出[ vasp .5.4.1/ src / subdftd3 . F ]（ 7 2 4 行 目 ） 外 部 関 数 volume を リ ン ク で き な い

実際の利用はないようなので、コメントアウト

2.3 FORTRAN90/SX 向け高速化
2.3.1 関数 rane のインライン展開
SX-ACE 向けの簡易性能解析情報ツールである、FTRACE を用いたプロファイリングにより、関数
rane が多数回呼び出されおり、この呼び出しオーバーヘッドのコスト割合が高いことが分かった。このオーバヘッドを削減するために、関数 rane のインライン展開を試みた。
インライン展開を行うにはコンパイルオプションに-pi を指定することで同一ファイル内のインライン展開が行われる。複数ファイル間のインライン展開を行うには -pi のサブオプション expin に対象ソー
スコード名もしくは対象ソースコードの保存されたディレクトリ名を指定する。-pi オプションの記述内
容についてはリスト 8 に示す。
関数 rane ではインライン展開を阻害する SAVE 属性の利用があったため、変数を新たな module を定義して変数の値を保持させることに変更した。リスト 5 に関数 rane のインライン展開を促進するために修正を行ったソースコードを示す。呼び出し側ではこの module を利用する宣言と変数の初期化を行う
ことで、インライン展開が可能となり、呼び出しオーバヘッドの削減を行うことが出来た。
リスト 5 インライン展開による高速化[ vasp .5.4.1/ src / random . F ]関数 raneの定義において、インライン展開を阻害するSAVE属性の利用をmoduleに変更する。
module i c a l l m を 設 定 し 、 i c a l l の 値 を 保 持 さ せ る 。
同 じ く イ ン ラ イ ン 展 開 を 阻 害 す る DATA文を代入文に変更する。
[ vasp .5.4.1/ src / main . F ]icallの値を保持するために、use

i c a l l m の 記 述 と i c a l l =0 で 初 期 化 。

— 38 —

SENAC Vol. 51, No. 1（2018. 1）

2.3.2 コンパイラ指示行によるベクトル高速化
ソースコード中の直後の DO ループに依存性がないことを示す、インテルコンパイラ向けの指示行（!DIR$ IVDEP）が記述されている箇所に、FORTRAN90/SX コンパイラ向けの同種の指示行（!CDIRNODEP）を挿入することで、コンパイラによる DO ループのベクトル化が行われるようにした（リスト6）。
リスト 6 コンパイル指示行によるベクトル化[ vasp .5.4.1/ src / fft3dlib . F ]! DIR$ I V D E P が 指 定 さ れ て い る ル ー プ 箇 所 （ 5 7 箇 所 ） に! CDIR N O D E P を 記 述 し て ベ ク ト ル 化 す る 。

2.4 SX-ACE 向け Makefile
VASP に付属の Makeﬁle の動作は、.F90 の拡張子を持つソースコードからプリプロセッサがオプションに従って.f90 のソースコードを生成し、このソースコードをコンパイルする、という流れでオブジェク
トファイルを生成する。

SX-ACE 向け高速化では、ファイル間のインライン展開を行うことで、関数・サブルーチン呼び出しの部分が DO ループ内で展開され、ベクトル化が促進されることがある。そのため make 時はコンパイ
ル時間を短縮するために、最適化オプションを最低レベルの-C ssafe として.f90 のソースコードを生成した後に、リスト 8 に示すスクリプトファイルで相互ファイル間のインライン展開の指定と、最大の最適化オプション-C hopt を指定してコンパイルとリンクを行った。
リスト 7 SX-ACE 向け makefile.include と src/makefile[ vasp .5.4.1/ makfile . inclued ]CPP_OPTIONS = - DMPI - DHOST =\ " NECSX \ " - DMPI - DPGF90 - Davoidalloc - DscaLAPACK - Du se _ co ll e ct iv e \- USX - Dpro_loop - DMPI_BLOCK =80000
CPP= sxcpp -P $ ( FUFFIX ) $ ( SUFFIX ) $ ( CPP_OPTIONS )FC= sxmpif90
FCL= sxmpif90
FFLAGS= - f5 - f2003
OFLAG= - Cssafe
BLAS= - lblas
LAPACK= - llapack
BLACS= - lblacsF90init - lblacs - lblacsF90init
SCALAPACK = - lscalapack - lblacsF90init - lblacs - lblacsF90init - llapack - lblas
OBJECTS= fft3dfurth . o fftmpi . o fftmpi_map . o fft3dlib . o
OBJECTS_O1 += fft3dfurth . o fftmpi . o
OBJECTS_O2 += fft3dlib . o
CC_LIB= sxcc
CFLAGS_LIB =FFLAGS_LIB = - Cvopt[ vasp .5.4.1/ src / makefile ] sxmpif90 の 最 適 化 オ プ シ ョ ン と 拡 張 自 由 形 式 の オ プ シ ョ ン に 変 更OFLAG_1 = - Cssafe - f5
OFLAG_2 = - Cssafe - f5
OFLAG_3 = - Cssafe - f5
OFLAG_4 =ﬀt3dlib.f90 と vdw_nl.f90 については、インライン展開のネスト数が 3 を超えた場合、コンパイラの最適化に要するメモリ不足により、コンパイル途中で処理が中断した。また force.f90 と fock.f90 については最適化オプションのレベルによって計算結果が異なることが分かった。このため、これらのファイルに対しては他のファイルとは異なるオプションを指定してコンパイルを行っている。

第１原理電子構造計算コード VASP の SX-ACE 向け最適化— 39 —

リスト 8 sxmake.sh（ vasp_stdをコンパイル・リンク）

cd ./ build / std
rm -f *. o *. F vasp
sxmpif90 - Chopt - f5
sxmpif90 - Chopt - f5
sxmpif90 - Chopt - f5
sxmpif90 - f5 - f2003
sxmpif90 - f5 - f2003- f2003- f2003- f2003-c - pi-c - pi-c - pi-c - pi-c - pi
nest =5
nest =5

nest =5 line =1000 expin =./ ./*. f90
nest =3 line =1000 expin =./ fft3dlib . f90
nest =3 line =1000 expin =./ vdw_nl . f90
line =1000 expin =./ force . f90
line =1000 expin =./ fock . f90

sxmpif90 -o vasp base . o mpi . o smart_ alloca te . o xml . o constant . o jacobi . o main_mpi . o ... \（その他必要なオブジェクトファイル）... fft3dfurth . o fftmpi . o fftmpi_map . o fft3dlib . o main . o \- Llib - ldmy - lscalapack - lblacsF90init - lblacs - lblacsF90init - llapack - lblas - llapack - lblas
cp ./ vasp ../../ bin / vasp_std

3 計算結果の比較
3.1 対象としたベンチマークファイル
ベンチマークファイルとして、スウェーデンのリンショーピン大学（Linköping University）国立
スーパーコンピューティングセンターのペーター博士（Peter Larsson, PhD）が、Cray XC-40（通称Beskow）において VASP のベンチマーク用に設定した系を用いる。*1 対象として GaAsBi-512 と呼ばれる系を使用した。これは、GaAs（格子定数 5.6537 Å）の 1 つの As のみを Bi で置換したものである。この系を表現するインプットファイルを用いて、SX-ACE と当センターの並列コンピュータ LX
R
E5-2695v2）との計算結果を比較した。LX 406Re-2 向けのコンパイルに使用
406Re-2（Intel⃝Xeon
R コンパイラおよびした makefile.include は vasp.5.4.1.24Jun15.tar.gz に同梱のものを使用し、Intel⃝

R
MKL ライブラリ（インテル⃝Parallel
Studio XE 2016）を用いてコンパイルとリンクを行った。

インプットファイルの 1 つである、計算条件のキーワードを指定するファイル、INCAR ファイルの共通パラメータは以下の通りである。
リスト 9 INCAR ファイルの共通パラメータ
SYSTEM = GaAs!
ISTART =0!
ICHARG =2!
PREC = Accurate !
ENCUT =313!
ISPIN =1!
ISMEAR = -5
ALGO = fast
LCHARG =. FALSE .
LWAVE =. FALSE .
LREAL = Auto
NELMIN =1
NELM =20
EDIFF =1 E -4
NBANDS =1536

系の名前始めから計算原子位置から電荷密度を推測
計算精度の設定
カ ッ ト オ フ エ ネ ル ギ ー [ eV ]スピンは考慮しない

3.1.1 SX-ACE での計算結果
SX-ACE は 1 ノードあたり 1CPU と 64GB の主記憶を搭載し、1CPU は 4 コアで構成される。1 ノードの理論演算性能は 276GFLPOS、メモリバンド幅は 256GB/sec である。

SX-ACE での実行時に INCAR に記載した並列化パラメータは以下の通りである。1 ノードあたり 4プロセス、32 ノードを用いて合計 128 プロセスで実行を行った。
*1 https://www.nsc.liu.se/~pla/blog/2015/01/13/vaspstudy-crayxc40/

— 40 —

SENAC Vol. 51, No. 1（2018. 1）

リスト 10 SX-ACE での並列化に関するパラメータ
NCORE =4
KPAR =4
NSIM =1
NPAR =32

SX-ACE の実行時にはゼロ割演算の警告を抑制するために、リスト 11 に示す環境変数をバッチリクエ
ストのスクリプトファイル内に追記した。
リスト 11 実行時環境変数# PBS -v F_ERROPT1 = " 252 ,253 ,0 ,0 ,2 ,2 ,2 ,2 "リスト 12 SX-ACE での実行結果（OUTCAR ファイルの抜粋）
total drift :-0.000019-----------------------------------------------------0.000019

0.000000

FREE ENERGIE OF THE ION - ELECTRON SYSTEM ( eV )--------------------------------------------------free energy
TOTEN =
-2110.33769038 eVenergywithout entropy =-2110.33769038

Total CPU time
User
SystemElapsed

used
time
time
time

( sec ):( sec ):( sec ):( sec ):

energy ( sigma - >0) =-2110.33769038

930.070929.5700.500
933.133（ イ タ レ ー シ ョ ン 回 数 は20 回 ）3.1.2 LX 406Re-2 での計算結果
LX 406Re-2 は 1 ノードあたり 2CPU と 128GB の主記憶を搭載し、1CPU は 12 コアで構成される。
1 ノードの理論演算性能は 460.8GFLOPS、メモリバンド幅は 119.4GB/sec である。
LX 406Re-2 での実行時に INCAR に記載した並列化パラメータは内容は以下の通りである。1 ノード
あたり 24 プロセス、8 ノードを用いて合計 192 プロセスで実行を行った。
リスト 13 LX 406Re-2 での並列化に関するパラメータ
NCORE =24
KPAR =4
NSIM =24
NPAR =2

リスト 14 LX406 での実行結果（OUTCAR ファイルの抜粋）
total drift :-0.000019-----------------------------------------------------0.000019-0.000000

FREE ENERGIE OF THE ION - ELECTRON SYSTEM ( eV )--------------------------------------------------free energy
TOTEN =
-2110.32291858 eVenergywithout entropy =-2110.32291858

Total CPU time
User
SystemElapsed

used
time
time
time

( sec ):( sec ):( sec ):( sec ):

energy ( sigma - >0) =-2110.32291858

1003.195
1001.277
1.919
1010.096（ イ タ レ ー シ ョ ン 回 数 は13 回 ）SX-ACE と LX 406Re-2 でそれぞれ実行した結果、イタレーション回数はそれぞれ 20 回、13 回と差異があったが、最終的な系のエネルギーの値は小数点第一位まで一致する結果が得られた。

第１原理電子構造計算コード VASP の SX-ACE 向け最適化— 41 —

4 VASP 実行オブジェクトの提供
サイバーサイエンスセンターでは VASP のソフトウェアライセンスを保有する利用者に、VASP.5.4.1の SX-ACE 向けバイナリと LX 406Re-2 向けバイナリを提供しています。利用方法については共同利用
支援係（連絡先は本誌の表紙内側を参照）までお問合せください。

5 まとめ
本稿では VASP.5.4.1 の SX-ACE 向け移植方法と高速化について述べた。当センターでは SX-ACE
は最大で 1,024 ノード（2,048 プロセス並列）で、60TB メモリを利用した実行が可能である。SX-ACE
で大規模な分子モデルの解析を行っていただき、利用者の研究ツールとしてご活用いただければ幸いで
ある。

参考文献
[1] P. Hohenberg and W. Kohn. Phys. Rev., 136:B846, 1964.
[2] Belgium Center for Molecular Modeling, Ghent University. Comparing Solid State DFT Codes,
Basis Sets and Potentials. https://molmod.ugent.be/deltacodesdft, 2016.
[3] G. Kresse and J. Hafner. Phys. Rev. B, 47:558, 1993.
[4] G. Kresse and J. Furthmuller. Phys. Rev. B, 54:11169, 1996.
[5] Jürgen Hafner. Journal of Computational Chemistry, 29:2044, 2008.
[6] R. Egawa, S. Momose, K. Komatsu, Y. Isobe, H. Takizawa, A. Musa, and H. Kobayashi.
Early evaluation of the SX-ACE processor. SC14: The international conference for High Performance Computing, Networking, Storage and Analysis, New Orleans, LA., 2014. The extended abstract can be found here: http://sc14.supercomputing.org/sites/all/themes/
sc14/files/archive/tech_poster/poster_files/post196s2-file3.pdf.
[7] S. Momose, T. Hagiwara, Y. Isobe, and H. Takahara. The brand-new vector supercomputer, SXACE. In Proceedings of 29th International Supercomputing Conference (ISC’14), pages 199–214.
Springer, 2014.
[8] 萩原孝 浜口博幸 山信田恒. スーパーコンピュータシステム SX-ACE のハードウェア. SENAC,
48-1:5–14, 2015.
[9] 山下毅 森谷友映 佐々木大輔 齋藤敦子 小野敏 大泉健治 岡部公起 江川隆輔 小林広明. スーパーコン
ピュータシステム SX-ACE の紹介. SENAC, 48-1:39–46, 2015.

— 42 —

SENAC Vol. 51, No. 1（2018. 1）［大学 ICT 推進協議会

2017 年度年次大会論文集より］

大規模科学計算システムにおける利用者プログラムの特性分析
大泉 健治 1), 山下 毅 1), 穂苅 寛光 2), 江川 隆輔 3,4), 滝沢 寛之 3,4), 小林 広明 4,3)
1) 東北大学 情報部情報基盤課
2) 日本電気株式会社 3) 東北大学 サイバーサイエンスセンター
4) 東北大学 大学院情報科学研究科
oizumi@cc.tohoku.ac.jp

Analysis of the characteristics of user programson Supercomputer system
Kenji Oizumi 1), Takeshi Yamashita 1), Hiroaki Hokari 2),
Ryusuke Egawa 3,4), Hiroyuki Takizawa 3,4) , Hiroaki Kobayashi 4,3)
1) Information Infrastructure Division of Information Department, Tohoku Univ.
2) NEC Corporation, 3) Cyberscienece Center, Tohoku Univ.
4) Graduate School of Information Sciences, Tohoku Univ.
概要
東北大学サイバーサイエンスセンターで提供している大規模科学計算システムにおいて、
実利用プログラムの性能状況をとらえるため、また性能改善の指標とするため、利用者プ
ログラムの特性分析を実施した。41 本の利用者プログラムを対象とし、プログラムが要求する B/F(Bytes/Flop)値から、性能のボトルネックとなる箇所を演算器律速、ADB バンド幅
律速およびメモリバンド幅律速に分類し傾向を分析した。その結果、メモリバンド幅律速のプログラムが約 7 割を占めることがわかった。

1 はじめに東北大学サイバーサイエンスセンター（以下、

2 運用状況
2.1 大規模利用課題の動向

本センター）の大規模科学計算システムは、日本

図 1 に 1999 年から本センターで取り組んでい電気株式会社（以下、NEC）製ベクトル型スーパーるセンター独自の共同研究、学際大規模情報基盤

コンピュータ SX-ACE を主力計算機として運用し、全国の研究者へサービスを提供している。また本センターでは 1999 年より、プログラムの高精度共同利用・共同研究拠点（JHPCN）課題および革新的ハイパフォーマンス・コンピューティング・イ

化、大規模化の支援を目的とした共同研究制度を

ンフラ（HPCI）課題採択数の推移を示す。本セン実施している。利用者、計算機科学を専門とする

ター独自の共同研究は恒常的に年 10 課題ほど実センター教員、技術職員、およびベンダー技術者

施されていることに加え、近年では JHPCN、HPCIが連携して利用者プログラムの高速化に取り組んを介した共同研究数が増加している。これは、セでいる。これらの高速化支援活動の取り組みの中で、プログラムの性能分析や高速化に関する個々の事例については多数報告してきた。本稿では、ンターとの共同研究を通して利用者プログラムが高度化・大規模化し、 JHPCN、HPCI 採択課題へと利用者プログラム 41 本について、性能状況をとら

ステップアップした結果だと考えることができ、え性能改善の指標とするため特性分析を行ったの

我々の継続的な高速化支援活動が一定の成果を挙で報告する。

げていることが分かる。

— 43 —

大規模科学計算システムにおける利用者プログラムの特性分析

3.1 SX-ACE のアーキテクチャ
表 1 にスーパーコンピュータ SX-ACE の諸元を、
図 3 に CPU コア構成をそれぞれ示す。1 ノードは 1
つの CPU とメインメモリから構成される。1 CPU
は 4 つのコアを搭載し 256 GFLOPS(64 GFLOPS×4)のベクトル演算性能を持ち、メインメモリは 64
GB、ADB(Assignable Data Buffer)は各コアで 1 MB
図１ 共同研究課題採択数

搭載している。また、SX-ACE では、冗長なメモリ
ア ク セ ス を 抑 止 す る た め に MSHR(Miss Status

2.2 利用ノード数の動向

Handling Register)を新たに実装している。以上

図 2 に、SX-ACE の利用ノード数と利用ノード時のような構成により、高メモリバンド幅を必要と間の分布を示す。SX-ACE は、2015 年から運用を開するプログラムで高い実効性能を達成することを始し、利用者はプログラムに応じて最大 1,024 ノ

特徴としている。しかしながら、CPU の演算性能

ードの計算資源を最大 1 ヶ月まで利用できる運用が向上している一方で、半導体チップへの実装にを行っている。動向分析の対象は、比較的新規利は制約があることからメモリバンド幅の性能向上

用が多い層であるセンターの共同研究課題の利用は難しくなっている。このためメモリバンド幅と者とした。

演算性能の比である B/F(Bytes/Flop)値は減少傾

利用ノード数は、システム導入当初は 1 ノード向にある。図 4 は本センターで提供してきたスーおよび 32 ノードまでの利用がほとんどであった

パーコンピュータの B/F 値を示している。1998 年

が、期間を追うごとに 33-256 ノードを利用した割に導入した SX-4 は B/F 性能値が 8 であるが、2015

合が高くなっており、利用者プログラムの並列化

年に導入した現行機種 SX-ACE の B/F 性能値は 1 でが進んでいる。これは、これまで我々が取り組んあり、代を重ねるごとに減少していることがわかできた利用者プログラムへの高速化支援によるとる。

ころが大きいと考えられる。また、利用者プログ
ラムの大規模並列化促進を目的に、利用ノード数が多いほど割安となる単価設定を行っていることも、高並列化を促進する一因と考えられる。

図 2 利用ノード数の割合（共同研究課題）

3 利用者プログラムの特性分析
本章では、演算性能とメモリアクセス性能の面からプログラム性能のボトルネックとなる箇所について、利用者プログラムを対象にして分析する。
性能分析にあたり、対象とするスーパーコンピュ
ータ SX-ACE の諸元について述べる。

表 1 SX-ACE の諸元— 44 —

SENAC Vol. 51, No. 1（2018. 1）の B/F 性能値は 1 となる。一方、ADB と各コアの間はそれぞれ 256 GB/sec のバンド幅性能を備えているため、CPU 全体での B/F 性能値は 4 となる。
表 2 は Prog B/F と Actual B/F の 2 つの指標を、
SX-ACE の B/F 性能値とそれぞれ比較し性能決定要因を分類したもので、図 6 は横軸を Prog B/F、縦軸を Actual B/F とし、要因ごとに色分けして示している。青色領域は演算器律速、緑色領域は ADB
バンド幅律速、赤色領域はメモリバンド幅律速と分類される。
図 3 CPU コア構成

図 5 SX-ACE のメモリアクセス
表 2 性能決定要因の分類表

図 4 歴代スーパーコンピュータの B/F 値
3.2 Prog B/F と Actual B/F
B/F 値は 1 演算(Flop)に必要なデータ量(Bytes)を示す指標である。プログラムによって要求される B/F 値は異なり、実行する計算機の B/F 性能値
より小さければ、演算器へ十分なバンド幅でのデータ供給が可能であり、CPU の演算性能を引き出すことができる。逆に大きければ、演算性能に見合ったバンド幅でのデータ供給をできず実効性能が低下する要因となる。プログラムの要求する B/F
値と計算機の B/F 性能値を比較することで、性能を律速している箇所が CPU の演算性能にあるか、メモリのデータ転送性能にあるかを推定することができる。

図 6 性能決定要因の分類図

プログラムが要求する B/F 値は Prog B/F と定義
され、プログラム上（ベクトルロード命令）のメ

3.3 分析結果

モリアクセス要素数をカウントすることで算出さ本センターの SX-ACE で実行されている利用者れる。一方、Actual B/F は、実際にメモリから CPU

プログラム 41 本を分析対象とした。表 3 は利用分に転送されたデータ量をもとに算出される（図 5）。

野とそのプログラム数である。

また、SX-ACE の B/F 性能値はつぎのとおりであ

まず、図 7 に各プログラムの Prog B/F と Actual

る。メモリと CPU の間は 256 GB/sec のメモリバン

B/F の測定値をグラフにて示す。Prog B/F を降順

ド幅性能を備える。4 コアすべてが同時にメモリに並べている。Prog B/F は最大 9.2、最小 0.1、

アクセスした場合、各コアで利用可能なメモリバ

平均 3.2 であり、その分散は大きいものの、全体

ンド幅は平均 64 GB/sec となり、メモリと CPU 間的に高いメモリバンド幅を必要とするプログラム

大規模科学計算システムにおける利用者プログラムの特性分析が多いといえる。一方、Actual B/F は最大 10.3、— 45 —

メモリから読み込まれる余分なデータ量が増える

最小 0.0、平均 2.4 であった。Prog B/F と比較すために、Actual B/F が Prog B/F よりも大きくなると、75 %のプログラムで、Prog B/F よりも Actual

るといえる。

B/F の方が小さいことが明らかになった。これは、

表 3 利用分野とプログラム数

ADB や MSHR が効果的に機能することによって、無駄・冗長なメモリアクセスを削減できているためと考えられる。残り 25%のプログラムでは、Actual
B/F の方が Prog B/F よりも大きな値を示している。
これは、ADB のデータ管理に関係していると考えられる。データはライン単位で ADB に読み込まれるために、実際にはアクセスされない余分なデー
タも読み込まれる。ADB にヒットしない場合は、

図 7 Prog B/F と Actual B/F の測定値
次に、性能決定要因の分類にもとづき利用者プ

メモリアクセス数が大きくなっていると考えられログラムを分類した結果の分布を、図 8 に示す。

る（橙色領域）。後者は、ADB の効果を確認できる

最も多いのはメモリバンド幅律速（赤色領域）の領域で、Actual B/F 値が相対的に小さくなるほどプログラム 28 本で、全体の 70%を占める。次いで、ADB のヒット率が高く効果が大きい（黄色領域）。

演算器律速（青色領域）が 11 本で 25%、残り 2 本

全体を分析した結果、メモリバンド幅律速のプが ADB バンド幅律速（緑色領域）となった。

ログラムが 70%と大部分を占め、本センターの利

演算器律速領域（青色領域）は、メモリバンド

用者プログラムにはメモリバンド幅性能を求める

幅および ADB バンド幅ともにデータ転送性能は十プログラムが多いことがわかった。しかし、ベク

分余裕があり、理論的には CPU のピーク演算性能

トル化率やベクトル長等、改善の余地が残されてを引き出せている。ADB バンド幅律速領域（緑色いるプログラムも存在していると考えられるため、領域）は、無駄・冗長なメモリアクセスが少なく、個々のプログラムのさらなる分析が必要である。

ADB が効果的に機能したデータ転送ができている。

また、Prog B/F と Actual B/F の 2 つの値から、メモリバンド幅律速のプログラムは、さらにプログラム性能のボトルネックとなる要因を比較

Prog B/F <=Actual B/F の 9 本、Prog B/F>Actual

的簡単に見積もることが可能なため、プログラム

B/F の 19 本に分けられる（図 9）。前者は、Actualの初期性能分析や性能改善の指標として活用でき

B/F が相対的に大きい領域で、ストライドアクセ

ることがわかった。

スやリストアクセス等でメモリロードが頻発し、

— 46 —

SENAC Vol. 51, No. 1（2018. 1）

参考文献

図 8 利用者プログラムの分布

図 9 メモリバンド幅律速領域

4 まとめ
本稿では、本センターの利用者プログラムを対象に Prog B/F と Actual B/F の測定を行い、本センターにおける利用者プログラムの特性を分析した。これらの分析結果により、利用者プログラムの求める性能が定量的に明らかになり、メモリバ
ンド幅律速のプログラムが多いことがわかった。
今後も引き続き動向を分析し、プログラムチュー
ニングのツールとして、また利用者のニーズに応えられるシステム設計、運用設計に活用していきたい。

[1] R.Egawa, K.Komatsu,S.Momose, Y.Isobe,
H.Musa, H.Takizawa, H.Kobayashi,Potentialof a modern vector supercomputer forpractical applications: performance
evaluation of SX-ACE, The Journal ofSupercomputing, Sep. 2017, Volume73, Issue
9, pp 3948-3976[2] 萩原孝 浜口博幸 山信田恒、スーパーコンピ
ュータ SX-ACE のハードウェア、東北大学サイバーサイエンスセンター大規模科学計算シ
ステム広報(SENAC) Vol.48,No.1、pp.5-14、2015．
[3] 工藤淑裕 横谷雄司、SX-ACE でのプログラミ
ング（ベクトル化偏）、東北大学サイバーサ
イエンスセンター大規模科学計算システム広報(SENAC) Vol.48,No.1、pp.15-32、2015．
[4] 東北大学サイバーサイエンスセンター、高速化推進研究活動報告 第 6 号、2015．
[5] 小野敏 他、スーパーコンピュータ SX-ACE の運用について、大学 ICT 推進協議会 2015 年度
年次大会．

— 47 —

SENAC Vol. 51, No. 1（2018. 1）［大学 ICT 推進協議会

2017 年度年次大会論文集より］

反応・相変化を伴う多分散系混相流シミュレーションコードの最適化
佐々木季広 2)，磯部

洋子 2)，

弘貴 3)，渡部広吾輝 3)，志村啓 3)，奥野

亜紀子 3)，江川隆輔 4), 5)，滝沢寛之 4), 5)，小林笠原
松尾大輔 1)，加藤航平 3)，
広明 5), 4)1) 東北大学 情報部情報基盤課
2) 日本電気株式会社
3) 慶應義塾大学
4) 東北大学 サイバーサイエンスセンター
5) 東北大学 大学院情報科学研究科
d-sasaki@cc.tohoku.ac.jp

Performance Optimization of Large-scale Simulation of Polydisperse
Multiphase Flow with Reaction and Retastability
Daisuke Sasaki1), Toshihiro Kato2), Yoko Isobe2), Hirotaka Kasahara3), Hiroaki Watanabe3),
Kei Shimura3), Kohei Okuno3), Akiko Matsuo3), Ryusuke Egawa4), 5),
Hiroyuki Takizawa4), 5), Hiroaki Kobayashi 5), 4)
1) Information Infrastructure Division of Information Department, Tohoku Univ.
2) NEC Corporation.
3) Keio University.
4) Cyberscience Center, Tohoku Univ.
5) Graduate School of Information Sciences, Tohoku Univ.
概要
東北大学サイバーサイエンスセンターでは，大規模科学計算システムを最大限活用する
ために計算科学者と本センターの計算機科学者が連携しながら，プログラムの高速化技法の研究・開発に取り組んでいる．本稿では，慶應義塾大学の研究グループが研究・開発している反応・相変化を伴う多分散系混相流シミュレーションコードのベクトル型スーパーコンピュータ SX-ACE における高速化について紹介する.スーパーコンピュータを利用することで，ワークステーションでの実行よりも圧倒的に大規模なモデルを高速に解析することが可能となり，粉塵爆発現象に対する安全防護技術の発展に貢献することが期待される.

図1

スーパーコンピュータ SX-ACE

— 48 —

1

はじめに
東北大学サイバーサイエンスセンター(以下，本センター)は，全国共同利用機関として大規模科学

SENAC Vol. 51, No. 1（2018. 1）

リバンド幅をそのコアで占有して利用でき，4Bytes/FLOP での実行が可能となっている．SX-ACE
は高ベクトル演算性能と，高データ供給性能により，高い実行性能を達成することができる．

計算システムの運用と，本システムを最大限に活本センターの SX-ACE は 2,560 ノードで構成さ用可能なコードの高速化技法や新しいシミュレー

れている．各ノード間は最大 4 GB/sec×2(双方向)ション技術の研究・開発に取り組んでいる．本セで接続された 2 段ファットツリーネットワークでンターでは 1997 年からさまざまな分野における
実アプリケーションの高速化支援を計算科学分野構成され，システム全体の性能は 706.6 TFLOPS となっている．利用者には，1 ジョブあたり最大
1,024 ノードの大規模な実行環境を提供している．

の利用者との共同研究を通じて行っている[1]．また，センター独自の共同研究に加え，全国の情報基盤センター等と連携して，学際大規模情報基盤

表 1 SX-ACE の諸元CPUあたりコア数
理論最大演算性能
最大ベクトル演算性能
ADB
メモリバンド幅(コア-ADB間）
メモリバンド幅(CPU-メモリ間）

4個
276GFLOPS
256GGLOPS
1,024KB×4
256GB/s×4
256GB/s

ノードあたり

CPU数
メモリ容量

1個
64GB

システムあたり

ノード数
CPU数
メモリ容量
理論最大演算性能
最大ベクトル演算性能

2,560ノード
2,560個160TB
706.6TFLOPS
655.4TFLOPS

共同利用・共同研究拠点(JHPCN)や革新的ハイパフォーマンス・コンピューティング・インフラ(HPCI)を構成し，多様なニーズに応える計算環境の提供も行っている．以上のように本センターでは，利用者である計算科学者と本センターの計算機科学の専門家が密接に連携しながら，科学・工学の恒常的な進歩を支える共同研究，利用者コードの高速化支援活動を推進している．
本稿では，慶應義塾大学の研究グループが研究・開発している反応相変化を伴う多分散系混相
流シミュレーションコードの高速化について紹介する．

3 反応・相変化を伴う多分散系混相流シ
ミュレーションコード
本コードでシミュレーションを行っている粉2

スーパーコンピュータ SX-ACE

塵爆発とは，可燃性の固体粒子が一定条件下で空気中に浮遊している状態で，何らかの現象により本センターに導入されているスーパーコンピュ

発火し爆発する現象のことである．工場等のもの

ータ SX-ACE を図１に，諸元を表 1 に示す．SX-ACE

つくり現場において，大規模な爆発を伴う火災がの CPU はこれまでのベクトルプロセッサと同様に，
高ベクトル演算性能と高メモリバンド幅を継承し発生する危険がある.この爆発を伴う火災は，工場内で発生した大量の粉塵が浮遊・堆積し，静電気などにより発生した火花に引火することが原因でている．また，CPU は 4 コアで構成され最大ベク発生する.また，この火災で生じる衝撃波により堆

トル演算性能は 256 GFLOPS である．各コアには容積していた粉塵が舞い上がりさらに大規模な火災量が 1,024 KB の ADB(Assignable Data Buffer)がに発展する．

搭 載 さ れ ， ADB と コ ア 間 の メ モ リ バ ン ド 幅 は
256GB/sec を有する．データが ADB 経由でアクセスされる場合は， データ転送速度と計算速度の比

慶應義塾大学の研究グループでは，粉塵爆発現象について，粉塵からの可燃性ガスの揮発及び炭
素個体の燃焼や気相燃焼，揮発ガス・既燃ガス・空気の混合を考慮しシミュレーションを行っていである B/F 値が 4 Bytes/FLOP となり，メモリ負荷る．また，爆発に伴って生じる衝撃波，爆風によの高いアプリケーションでも高い実行効率で実行る粉塵の巻き上げについて解明するため，シミュ可能となっている．また，単一コアで演算を実行

レーションコードの開発に取り組んでいる．

する場合には，CPU-メモリ間の 256 GB/sec のメモ

反応・相変化を伴う多分散系混相流シミュレーションコードの最適化

4

コードの高速化

4.1

SX-ACEにおける初期性能— 49 —

慶應義塾大学の研究グループが開発したプログ
ラムの初期性能解析結果を表 2 に示す．性能解析には SX-ACE の簡易性能解析ツール FTRACE を用いた．また，測定データのグリッドサイズは 1,901
×161 である．
初期性能解析の結果に基づき，全体の実行時間に対して約 90%を占めているサブルーチン A～C に

図 2 サブルーチン A,B のオリジナルコード

ついて高速化を行うことを検討する．また，将来の大規模並列実行を見据え，コードの並列化も検両サブルーチンの最内ループは 3 重ループの最内ループとなっているため，最内ループの外側の討する．

ループでベクトル化することで性能の改善が期待
表 2 初期性能解析できる．今回は図 3 のように unroll 指示行を挿入

ベクトル ベクトル化率
長subroutine A
8039.993
91.75
21.9
subroutine B
2314.420
82.60
18.8
subroutine C
1164.311
13.23
12.3
other subroutines 1850.370
total13397.172
88.92
26.2

し，最内ループを展開，最内ループの外側のルーサブルーチン名 実行時間

プでベクトル化することで性能改善を試みる．

以下では，本コードに適用した高速化について
述べる．
4.2

サブルーチン A，B

表 2 からサブルーチン A，B はベクトル長が非常に短いことが分かる．ベクトルプロセッサは複数

図 3 サブルーチン A,B の高速化後のコードのデータに対して一括で演算を行うことができ，最内ループの反復(イテレーション)回数であるベ

高速化後の性能を表 3 に示す．指示行を挿入す

クトル長が長いほど演算の効率が高くなる．

ることでベクトル長がサブルーチン A で 21.9 から

SX-ACE ではベクトル長を 256 の時に演算効率が最225.9 に，サブルーチン B で 18.8 から 250.9 まで大になる．

大きく改善され，サブルーチン A で約 44 倍，サブサブルーチン A，B の内部を詳細に解析したとこ

ルーチン B で約 17 倍の性能向上を得た．

ろ，両サブルーチンは同様のループ構造となっており，非効率なベクトル化が行われている箇所があった(図 2)．このループの反復回数は 10 であり，ベクトル化されているものの，SX-ACE の理想のベ
クトル長より非常に短く，ベクトル化の恩恵を受けていない．

表 3 高速化の効果

実行時間 ベクトル ベクトル 加速
化率
長
率
[sec]
before 8039.993
91.75
21.9
44.17
subroutine A
after
182.020
99.18
225.9
before 2314.420
82.60
18.8
subroutine B
16.96
after
136.431
99.97
250.9

— 50 —

4.3

SENAC Vol. 51, No. 1（2018. 1）

サブルーチン C

表 2 からサブルーチン C はベクトル化率が低
く，ベクトル長が短いことが分かる．サブルー
チン C の内部を詳細に解析したところ，ループ内に依存関係があり，最内ループがベクトル化
されていないことが分かった(図 4）．また，最内ループは外側のループ毎にループ長が異なり，かつそのループ長は短い．外側ループのループ長が長いため，外側ループを最内ループに移動することで性能改善が期待できる．

図 5 サブルーチン C のベクトル化後のコード

図 4 サブルーチン C のオリジナルコード
そこで，ループの一部を分割し，最内ループのループ長のうち最大になるものを検索する．次に，最内ループの終端を検索した最大ループ長に変更
し，ループ間の依存関係を解消する．最大ループ長に満たないループはマスクにより処理を行わないようにしている．依存関係が解消したため，ループ交換を行い，外側のループを最内側ループに移動することでループ長を長く確保し，ベクトル化率を高めた(図 5)．さらに，ベクトルプロセッサに搭載されている ADB を有効に活用できるように最内ループのブロック化を行う． SX-ACE は１

図 6 サブルーチン C の高速化後のコード命令で 256 要素を同時に演算できるため，ブロッ
ク長は 256 とした(図 6）．

高速化後の性能を表 4 に示す．高速化によりベクトル化率が 12.23%から 98.91%に改善し，ベクトル長が 12.3 から 151.3 に伸びることで，サブルー
チン単体で約 6 倍の性能改善を得ることができた．

表 4 サブルーチン C の高速化の効果
実行時間 ベクトル ベクトル 加速率
化率
長
[sec]
before 1164.311
13.23
12.3
6.11
after
190.607
98.91
151.3

反応・相変化を伴う多分散系混相流シミュレーションコードの最適化— 51 —

以上の高速化を適用後のコード全体の性能を表5 に示す．シングルコア実行で高速化前に比較し約 17 倍の性能改善を得ることができた．

表 5 高速化後の性能

before
after

実行時間
[sec]
13397.172
790.884

加速率16.94

図 7 並列化性能
4.4

並列化

本コードはベクトル化を行うことでシングルコ

5

まとめ

アでの性能を高めてきた．さらなる性能の向上の

本稿では反応・相変化を伴う多分散系混相流シ

ためには，ネットワークに接続される複数のノー

ミュレーションコードに，本センターが今までにドを利用し，多数のコアで計算をする並列化を検取り組んできた高速化支援で得られた高速化技法

討する必要がある．並列化にはノード内のコア毎を適用し高速化を行った事例について紹介した．

に プ ロ セ ス を 処 理 す る MPI(Message Passing

本センターの提供しているスーパーコンピュータ

Interface)並列を行う Flat MPI を用いた．

SX-ACE の性能を最大限に引き出すためには，ベク

本コードの計算領域は 1,901×161 の 2 次元空間

トル化率を高く保つ必要がある．また，SX-ACE はで，この空間中に粒子が存在する．この 2 次元空

従前のベクトル型スーパーコンピュータより，さ間中の熱や圧力等の計算後，粒子の相互作用の計

らにノードを多数接続することで性能の向上を図算を行う．並列化をするにあたり，計算領域を分る構成となっている．従って，今後の高速化支援

割し分割後の各領域に MPI のプロセスを割り当て，においては，シングルコアでの高速化に加え，複各領域の熱や圧力等と存在する粒子に対して相互数のノードのコアを利用する並列化を行うことで作用の計算を行うように並列化を行った．

性能を引き出していくことが重要となる．

図 7 に並列化後の性能を示す．2 次元分割を行

本コードは，高速化により数時間程度を要する

うと計算の初期で粒子が存在しない領域が多数存

シミュレーションを数分程度で計算することが可在し，プロセス間のインバランスが大きくなるた

能となった．今後，さらなる大規模なモデルでシ

め，計算領域を x 方向(1,901)での 1 次元分割とす

ミュレーションを行うことで安全防護技術の発展

る．高速化後のシングルコア実行に比較し，64 プに貢献していきたい．

ロセス実行で約 9 倍の性能改善を得るにとどまった．実行後，時間が経過すると，各領域に存在する粒子数に偏りが発生し，プロセス間のインバランスにより，性能向上が抑制されていると考えらえる．今後，さらなる高速化のため，多数の粒子が存在する領域をさらに細かく分割することや，ノード内の並列化に OpenMP を用いた Hybrid MPIなどの検討が必要である．

参考文献[1] 東北大学サイバーサイエンスセンター， 高速化推進研究活動報告，第 6 号，p8-12，(2015)．

— 52 —

SENAC Vol. 51, No. 1（2018. 1）［大学 ICT 推進協議会

2017 年度年次大会論文集より］

デジタル時代の教育・研究を支える基盤としてのeduroam と次世代ホットスポット
後藤英昭 1) ,中村素典 2) ,曽根秀昭 1)

1) 東北大学サイバーサイエンスセンター
2) 国立情報学研究所
hgot@cc.tohoku.ac.jp

eduroam and the Next Generation Hotspot as the Infrastructure
for the Research and Education in the Digital Age
Hideaki Goto1) , Motonori Nakamura2) , Hideaki Sone1)
1) Cyberscience Center, Tohoku University
2) National Institute of Informatics

概要
学術系の無線 LAN ローミング基盤である eduroam は，世界 89 か国，国内 189 機関 (2017 年 9 月現在, 前年末より +28 機関) に成長した．初等・中等教育機関や，市街地の公衆無線 LAN への導入も進められており，デジタル時代の教育・研究を支える基盤としての重要性が高まってきている．本報告では，世界規模の巨大な通信事業として，また，教育インフラのひとつとして eduroam を捉え，現在の開発状況を紹介する．世界で展開が進んでいる
次世代ホットスポット (NGH) と eduroam を連携する仕組みを開発し，City Wi-Fi Roaming トライアルに参加することで，その有効性を確認した．
200

1 はじめにeduroam (エデュローム) は，世界の大学や研究所等において，キャンパス無線 LAN の相互利用を実現する．日本では，eduroam JP の名称で 2006 年に実証実験的なサービスが始まり，2016 年 4 月からは国立# of institutions

国 際 的 な 学 術 無 線 LAN ロ ー ミ ン グ 基 盤 で あ る

175eduroam JP

150
125
100
75
50
25
0
'09.12 '10.12 '11.12 '12.12 '13.12 '14.12 '15.12 '16.12 '17.09

情報学研究所 (NII) の正式サービスとして運用されて

図 1 eduroam JP 参加機関数の推移

いる [1]．

eduroam は当初，主に大学等でキャンパス無線 LANの相互利用を実現するためのシステムとして開発・導入が進められていた．近年では，[2, 3] で既報のように，市街地でのサービス提供が増加し，初等・中等教育機関への導入も進められるようになり，社会全体における教育・研究を支える基盤としての様相を呈してきている．

ルに参加したので，概要と成果を示す．

2 教育・研究の基盤としての無線 LAN シス
テム2017 年 9 月時点で，eduroam は世界 89 か国 (地域)

本報告では，教育インフラのひとつとして eduroamに普及し，国内では前年末より 28 機関増となる 189を捉え，現在の開発状況を紹介する．特に，世界で機関が参加するに至った．図 1 に国内の参加機関数の展開が進んでいる次世代ホットスポット (NGH) と推移を示す．国内には約 1,200 の高等教育機関がある

eduroam を連携させることによって，ネットワーク利ことから，普及率で見ると 16% 弱であるが，7 大学や用環境が大幅に改善される可能性がある．その実現を

主要私立大学にはほぼ導入されており，eduroam は重目指して，2017 年の City Wi-Fi Roaming トライア

要な教育・研究基盤として認識されていると言える．

デジタル時代の教育・研究を支える基盤としての eduroam と次世代
ホットスポット— 53 —

2015 年度年次大会 [2] 等で既報のとおり，大学のり，電気や水道のように，行政が提供すべきインフラ

キャンパス内に留まらず，市街地や博物館等の施設，としての側面が打ち出されている．ニューヨークの空港・鉄道駅，大学関係者の出入りの多い病院などでLinkNYC [5] では，高速無線 LAN 基地局と大型デジの eduroam サービスの提供も，世界各地でよく見ら

タルサイネージ, 公衆電話, USB 充電端子を備え，ギれるようになった．特に，北欧各地の空港や，欧州原

ガビット対応のバックボーン回線を持つキオスク端末

子核研究機構 (CERN) の最寄りのジュネーブ空港は，が多数，主要な通りに並べられている．LinkNYC は，

Twitter や Facebook でもしばしば話題に上るほど好ニューヨーク市政府が中心となって，様々な企業が参評である．国内でも，世界的に早い 2010 年に，関東地画して作られているが，観光客ばかりではなく，地元区の市街地で eduroam サービスが開始された [4]．現

市民へのネットワーク提供も主な目的のひとつとして在もカフェや貸会議施設等で利用可能となっており，

掲げられている．特に，自宅にネットワーク回線を引

教職員・学生から好評を得ている．

けない層にネットワークアクセス手段を提供すること

キャンパス外での eduroam サービスは，単に利用場の福祉的効果も，期待されている．

所の増加のみならず，学習・研究活動にも影響を与え

欧州では，2020 年までに EU 全域の町村にフリーる可能性が高い．例えば，興味ある文献などを出張先

Wi-Fi を展開するプロジェクト，WiFi4EU が立ち上で見つけた際に，個人契約の携帯電話のテザリングでがった [6]．このプロジェクトでは，居住地や収入の額

大きなファイルをダウンロードすることは，学生はもに関わらずネットワーク利用を可能とすることがうた

ちろんのこと，研究者でも躊躇することが多いと思わわれており，これにより公共サービスへのアクセス手れる．思い立った時にすぐに学術リソースにアクセス段を確保するとともに，デジタル時代の社会福祉に貢できる環境を提供することは，教育・研究の支援の観点献することが志向されている．

でも有用なものと考えられる．学外での調査等の多い

以上のように，公衆無線 LAN についても，公共性研究者にとっては，特に有効であろう．収集したデーの高いインフラとしての重要性が高まり，学校に限ら

タを最寄りのカフェで整理し，大学の計算サーバを用ない教育活動・文化活動への貢献は小さくないと考えいてその場で解析するといったことも可能となる．

られる．

なお，市街地やカフェ等店舗では無償で利用できる
公衆無線 LAN，通称「フリー Wi-Fi」が提供されていることがある．現在のフリー Wi-Fi は暗号化のない接3 次世代ホットスポットと eduroam の連携
3.1 フリー Wi-Fi のセキュア化動向

続方式が一般的であり，偽基地局に誘導されやすいと利用者が無償利用できる，いわゆる「フリー Wi-Fi」

いうセキュリティ上の問題が大きい．最近では，盗聴は，観光地を中心に国内でも広く展開されている．しを防ぐための VPN サービスを契約し，容易に利用でかしながら，前章に述べたように，現在の多くのフきるようになったが，偽基地局への接続は避けられな

リー Wi-Fi には盗聴と偽基地局の問題がある．

いため，セキュリティの観点で十分な対策ではない．

幾つかの国では，フリー Wi-Fi においても，暗号化

一方，eduroam が採用している 1X 認証では，正規ののない SSID(Service Set Identiﬁer) に加えて，1X 認

運用者が設置した基地局でなければ利用者認証を成功証用の SSID が併設されていることがある．例えば，させることができず，認証に失敗すると端末は IP ア

シンガポールには Wireless@SG というフリー Wi-Fi

ドレスの取得さえできないことから，偽基地局への対サービスがあり，暗号化無しの SSID=“Wireless@SG”

策も可能である．すなわち，市街地における eduroamに加えて，1X 認証対応の SSID=“Wireless@SGx”がサービスは，端末を偽基地局から守り，暗号化により吹かれている．このサービスでは，通信事業者 4 社が無線通信を保護するといった，セキュリティ面での利認証連携しており，一か所で取得したアカウントがど点も大きい．

の事業者の基地局でも利用できる．ニューヨークやサ学術利用に限らない，一般市民向けの公衆無線 LAN

ンフランシスコ，サンノゼなどの都市では，1X 認証をについて見ると，世界の様々な都市において，フリー

ベースとした新しいシステム「Passpoint」(次節参照)

Wi-Fi が提供されている．従来，フリー Wi-Fi は買が導入されている．

い物客や観光客が主な対象であったが，最近の幾つ

一方，日本国内のフリー Wi-Fi は，暗号化無しのかのプロジェクトでは，住民へのネットワークアク

SSID にキャプティブポータルを組み合わせたシステ

セスの提供も目的のひとつに積極的に挙げられておムが一般的であり，セキュア化が進んでいない状況に

— 54 —

SENAC Vol. 51, No. 1（2018. 1）

ある．

(NGHSIG)」を発足させた [7]．

3.2 次世代ホットスポット (NGH) の動向

3.3

NGH 連携による市街地 eduroam サービスの拡充

世界の Wi-Fi 業界の動向を見ると，2014 年頃より北eduroam では，セキュリティ確保の観点から，ウェブ

米を中心に，次世代ホットスポット (NGH, Next Gen-

認証の利用が禁止されている．また，利便性確保の観

eration Hotspot) の 導 入 が 活 発 化 し て き た ．NGH

点から，原則として SSID=“eduroam” を使うことがは，Wi-Fi の業界団体である Wireless Broadband Al-

求められている．これらの制約により，ウェブ認証を

liance (WBA) と Wi-Fi Alliance が共同で推進してき基本とする現行の公衆無線 LAN システムと eduroam

た規格である．NGH では，端末の認証と接続を自動を認証連携させることは，不可能であった．もし将来，化し，認証情報を業者間で相互利用 (ローミング) でき

市街地のフリー Wi-Fi に Hotspot 2.0(NGH) が導入るようにすることで，公衆無線 LAN の利便性を携帯されると，セキュリティ及び利便性の問題が解決され，電話並みに高め，また，携帯電話とのシームレスな共

eduroam サービスの市街地展開が容易になる可能性存を図ることが志向されている．端末の自動認証・接がある．

続には，IEEE802.11u(以下 802.11u) と IEEE802.1X

eduroam における従来の市街地 (キャンパス外) 展

を組み合わせた，Hotspot 2.0 と呼ばれる仕様が用いら

開では，SSID として “eduroam” を基地局に追加し，れる．Hotspot 2.0 とオンラインサインアップ (OSU)

必要に応じてアクセス回線も用意する必要があった．

などの機能をまとめて，Wi-Fi Alliance の認証を取得

公衆無線 LAN の基地局は，既に複数の事業者の SSID

したものは，Passpoint と呼ばれる．サンフランシスが相乗りしていることが多く，追加の SSID を吹く

コなどのフリー Wi-Fi や，無線 LAN 事業者が提供す余裕がない，高額な運用経費が必要になるといったる公衆無線 LAN サービスなどで，Passpoint が導入さ問題があり，eduroam サービスの展開は容易ではなれている．ニューヨークの LinkNYC でも，Passpoint

かった．

が併設されており，一度の利用者登録のみで，以降はHotspot 2.0 の場合，SSID は任意となるが，802.11u

セキュアなフリー Wi-Fi が自動接続で利用できる．

の部分で eduroam 用のプロファイルが必要となる．

Hotspot 2.0 を利用した自動接続において，後半の802.11u では，ローミング・コンソーシアムを示す部分は eduroam が採用している 1X 認証と同じもの

NAI realm (NAI: Network Access Identiﬁer) が必要である．前半の 802.11u の部分で自動接続が実現されとなるが，GÉANT(eduroam の国際運用の中心機関)

る．従来の 1X 認証では接続先の基地局の SSID を端において，“eduroam.org” が既に予約されている．た

末に登録しておく必要があるが，Hotspot 2.0 ではこ

だし，認証方式などの細かい設定は確立されたものがれに代わって，利用者が契約している通信事業者のプ

なく，開発や検証が必要である．

ロファイルが登録される．基地局側にも通信事業者のフリー Wi-Fi は，個々の店舗や商店街組合，自治体

プロファイルが登録されており，SSID を問わず，合致など，基地局設置場所のオーナーの負担や，広告収入するプロファイルが存在する基地局に対して，1X 認証

によって運営されており，ビジネスが局所的に閉じてによる接続が試行される仕組みである．フリー Wi-Fi

いることが多い．そのため，認証連携により無線 LANで用いられる SSID は，しばしば場所のオーナーなどシステムを利用させてもらうのに，経費がそれほどかの広告を兼ねているため，統一されたものを使うこと

からなくて済む可能性が高い．実際に，次章で述べるは難しい．Hotspot 2.0 では，SSID に依存しない接続

City Wi-Fi Roaming トライアルでは，都市のフリーを実現したことにより，異なる場所でも端末の再設定

Wi-Fi が世界各国の携帯電話会社に無償でサービスを無しにシームレスな利用が可能となる．

提供するモデルが模索されている．フリー Wi-Fi を提日本国内では，2016 年時点で，携帯電話や公衆無供する側にとっては，利用者が匿名よりも，個人に紐線 LAN の事業者に NGH を推進する動きが見られな付けされたアカウントの方が，ネットワーク提供の責かった．また，フリー Wi-Fi のセキュリティ問題が指任を問われにくいという利点がある．また，既存のア

摘されているにも関わらず，セキュアな接続手段の導カウントを利用できるため，サインアップや本人確認

入が進んでいなかった．このような状況の下，国内のの手間が省けることも利点になる．

公衆無線 LAN のセキュア化と NGH 導入を推進する目的で，本稿第一著者が発起人・幹事となって，2017
年 1 月に「セキュア公衆無線 LAN ローミング研究会

— 55 —

デジタル時代の教育・研究を支える基盤としての eduroam と次世代
ホットスポット

4 City Wi-Fi Roaming トライアルにおける実証実験
4.1 WBA City Wi-Fi Roaming
City Wi-Fi Roaming は，WBA が主催し，世界各地の都市で提供される公衆無線 LAN を NGH 基盤で結ぶことにより，ひとつのアカウントで相互利用 (ローミング) できる環境を構築しようとする，世界規模のトライアル (プロジェクト) である．2016 年にニューヨーク，サンフランシスコ，サンノゼ，シンガポール
を結ぶ第 1 回が開催された．第 2 回目となる 2017 年
は，6 月 20 日の World Wi-Fi Day を起点として，8月末までの期間，20 程度の都市を結んで開催された(図 2) [8]．本トライアルは，世界の通信事業者と共に認証連携の実証実験を行う，大変良い機会であると考えられた．
c
図 2 WBA City Wi-Fi Roaming (WirelessBroadband Alliance / World Wi-Fi Day 2017)

我々は，セキュア公衆無線 LAN ローミング研究会として，主に以下の目的をもってトライアルに参加した [9]．

DNS-based realms• WBA 及び世界の通信事業者との接点を作り，今XXroamOpName-based realms
anyroam

WBA/City Wi-Fi, etc.

後の様々な活動の足がかりとする．

• 研究会で構築中の NGH テストベッドシステムをSP

携の実証実験や，Hotspot 2.0 の技術に関する情報収集・交換を行う．

• 国 内 の 事 業 者 に 参 加 を 呼 び か け ，Hotspot

Tohoku University

JP hub

世界の NGH 基盤に接続することにより，認証連
Shopping Mall
Wi-Fi,
Seaport Wi-Fi,
etc.

IdP/SP

IdP

IdP/SP
IdP/SP

2.0/NGH の啓発と接続支援を行う．

Basic
connections
Optional/example
connections(bilateral)

図 3 JP hub と NGH テストベッド• 多数のローミング・コンソーシアムを結んで，世界規模のローミング基盤を実現するための，イン
ターローミングの認証連携アーキテクチャを開発する．eduroam と govroam を扱える仕組みを考案し，世界の NGH 基盤の上で実証実験システムを構築，評価する．

4.2 NGH テストベッドと JP hub
研究会では，NGH の国内基盤を構築するために，認証連携の中心となる JP hub を立ち上げ，これに幾つかの機関のサーバを接続して，NGH テストベッドを構築した．システムの概要を図 3 に示す．
執筆時点の仕様では，NGH テストベッドに接続されたすべての SP(Service Provider, 無線アクセスを提供する機関) と IdP(Identity Provider, アカウント
情報を提供し認証処理を行う機関) は，相互利用が可能となっている．

4.3

eduroam-NGH 連携の技術的課題

通常，ひとつの通信事業者はひとつないし少数のレ
ルム (realm) を利用しており，認証連携では，ユーザ名の後ろに@とレルム名を付けたものがユーザ ID として用いられる．認証連携ネットワークは，RADIUS
プロキシが相互接続されたものであり，利用者端末から送られてきた認証要求をレルムを頼りにルーティン
グ (経路制御) して，利用者の IdP(ホーム IdP と呼ばれる) まで転送する．
事業者間のローミングにおいて，SP ではレルムを見ることによって，ローミング契約のある事業者の利用者かどうかを判断できる．例えば，端末から送られてきたユーザ ID に boingo.com というレルムが含まれていれば，米 Boingo 社の利用者であることが判る．
レルムが domru.ru であれば，露 ER Telecom 社だ

— 56 —

SENAC Vol. 51, No. 1（2018. 1）と判る．一方，eduroam では，レルムを見ただけでNGH hub operators

NGH
hub

は eduroam の利用者であると判断するのが難しい．

eduroam は巨大なローミングコンソーシアムであり，

NGH
hub

NGH
hub

世界中の参加機関のレルムが存在するが，機械的に判断できる共通の文字列がレルムに含まれているわけで

UK

NL

BE

US

JP

はない．例えば，SP で tohoku.ac.jp というレルムを見たとき，ac.jp が日本の大学関係のレルムだという事前知識がなければ，認証要求を eduroam のプロキ

All AuthN requests with .jp
realms except the known
operators’.

Federation-level
Hub layer

UK
NL
BE
JP
UK
NL
BE
US
UK
NL
BE
US
JP

シに転送するという判断ができない．

図4

eduroam をベースに開発された，政府関係機関向け

XXroamgovroam

インターローミング・アーキテクチャのローミング基盤である govroam [3] でも，同様の問題がある．膨大な数のレルムとコンソーシアムの対応
表をすべての SP が持ち，常に最新の情報にアップデー

トに含まれる Operator-Name 属性値を調べることにトすることは，現実的ではない．そのため，eduroam

よって，ローミング契約のある SP かどうかを判断し，や govroam を NGH 基盤に接続しようとすると，認認証処理の可否を決めることが可能である．

証要求を正しいコンソーシアムのプロキシに転送する

4.4 実証実験

ための，効率的な仕組みが必要となる．

2017 年 5 月頭に City Wi-Fi Roaming のキックオ

そ こ で ，以 上 の 問 題 に 対 処 で き る よ う な ，イ ン

フの遠隔ミーティングがあり，参加機関・都市の募集

ターローミングを実現するアーキテクチャを開発し，と NGH 基盤への接続作業が開始された．世界の 4 つ

TNC17 Mobility Day で発表した [10]. 図 4 にシステの NGH hub オペレータのうち，前年のプロジェクト

ムの全体構成を示す．各国には eduroam や govroam，

で実績のある GlobalReach Technology 社 (以下，GR

その他のローミング基盤 (図中の XXroam) が存在し，社) をパートナーに選び，GR 社の hub に JP hub をそれぞれに世界規模で接続されたコンソーシアムが接続して，これにより世界の NGH 基盤への接続を果あると仮定する．eduroam や govroam のレルムは，たした．

DNS ドメイン名に対応するように設計されているこ

最終的には，40 程度の携帯電話会社が 3GPP ネッとから，ccTLD (country code Top Level Domain)

トワークを介した EAP-SIM 認証を提供し，約 20 のごとにルーティングを分散処理することを考える．各都市が参加表明を行った．我々は国内の NGH 基盤と国に代表となるプロキシ (以下，代表プロキシ) を 1 台

eduroam の，ふたつのコンソーシアムについて，こ用意し，他のコンソーシアムと同様に，世界で水平にれらのアカウントを供給する立場に加えて，参加する

認証連携した Hub layer を設ける．代表プロキシは，事業者を募って国内数か所に基地局を立てる計画を立当該国内のすべてのコンソーシアムとレルムの対応表てた．

を持ち，適当な間隔で更新しているものとする．国内
では図 3 の JP hub がこの役割を担う．

当初，末尾が .jp のすべてのレルムについて，認証要求を JP hub に集約できるように，NGH hub オペ

例えば，日本国内では，すべての ac.jp レルムはレータに依頼した．研究会では nghsig.jp というレル

eduroam に所属しているとみなせるため，該当する機ムを使用するが，これも .jp に含まれる．ただし，既関のレルムをひとつの正規表現でまとめることができ知の事業者のレルムや，eduroam/govroam が含まれる．go.jp については，eduroam の対象機関と政府関ることのない co.jp, ne.jp などのレルムは除外した．

係機関が混在しているため，個別のレルムを対応表に

しかしながら，NGH hub オペレータ側で正規表現を登録する必要がある (注: 日本はまだ govroam に加盟

用いたレルム・マッチングがすぐには実装できないとしていない)．また，kek.jp のように，汎用 JP ドメイ

いうことで，以下の 3 種類を固定レルムとして転送しンを使用している所は，執筆時点で 10 機関である．

てもらうように計画変更せざるを得なかった．

このアーキテクチャでは，レルムの対応表が各国に分散されるため，スケーラビリティの問題が軽減される．SP 側で利用者のコンソーシアムを知ることは依然として困難であるが，IdP 側は認証要求のパケッ• tohoku.eduroam.jp(東北大学の eduroam 実験用アカウント)• tohoku.nghsig.jp(NGHSIG, 東北大学用アカウント)

デジタル時代の教育・研究を支える基盤としての eduroam と次世代
ホットスポット• wwd.nghsig.jp— 57 —

し，Passpoint プロファイルを作成するシステムを期(NGHSIG, World Wi-Fi Day 配布用 (デモ) アカ

間内に用意できなかったことから，RADIUS のテス

ウント)

トコマンドや，端末上で 1X 認証の手動設定を利用し基地局については，まず東北大学サイバーサイエン
スセンターの実験室に Hotspot 2.0 対応のものを一基設置して，屋外からも利用できるようにした．一般に，大学等の既存の eduroam 基地局及びネットワーク回
線は，教育・研究目的の利用に限定されていることから，一般利用者のローミング利用を受け入れることは困難である．そのため，上流のネットワーク回線には

NII より実験用として借りたものを用いた．基地局のベンダより協力の申し出があり，追加で東京にも基地局を設置することができた．この他，7 月上旬に 3 か所程度の基地局設置を予定したが，基地局の入手や現場の調整が大幅に遅れたことから，一般に公開できたのは 8 月に入ってからであった．最終的には，仙台 1，
東京 2, 京都 2 の計 5 か所でサービスを提供できたが，京都の 1 か所は旅館であり，eduroam を提供する国内初の宿泊施設となった．

4.5 接続試験と評価
最初に，GR 社が提供する Passpoint OSU (Online

Sign-Up) を利用して Passpoint プロファイルを取得
し，iPhone や Android スマートフォンから無線 LANが利用できるかどうかを試した．iPhone では概ね良好に自動接続が可能であったが，Android 6 系はプロ
ファイルが読み込めなかったり，プロファイルが入っても自動接続しなかったりなど，トラブルが多発した．
基地局のファームウェアにもバグが多いことが判明した．デバイスの実装や相性の問題が多いことは WBA
でも把握しており，執筆時点でもテストと改良が続けられている．
スマートフォンを用いた EAP-SIM 認証について
は，AT&T 等の SIM を用いて，正常に認証・接続が可能なことを確認した．City Wi-Fi Roaming の参加都市では，観光客が自国の携帯電話の SIM を使い，現地で利用登録なしに，無償かつセキュアに公衆無線 LANを利用できる．これは，観光客誘致の視点でも，ネッ
トワークの不正利用防止の観点でも，大きな利点であると考えられている．

Wi-Fi プロバイダとの認証連携については，著者が契約している Boingo の Passpoint プロファイルを利用して，国内の基地局で正常に認証・接続できることを確認した．
国外からの接続試験は，以下のように行った．ただて，認証連携の部分のみを評価した．まず，GR 社からの認証テストが正常なことを確認した．次に，ER

Telecom の担当者にテストを依頼して，ロシアからも正常に認証が通ることを確認した．
海外のフリー Wi-Fi での接続試験を行うために，参加都市に在住あるいは訪問予定の人をボランティア
として募集したが，期間内に適切な人が現れなかった．そのため，著者自らが旅行を計画し，イギリスのバーミンガム及びリーズにて，接続試験を行った．
もし Passpoint が利用できるなら，利用者が現地のSSID を知らなくても自動接続されるはずである．しかし，今回はプロファイルを用意できなかったため，

City Wi-Fi Roaming に対応した基地局を現地で探し出すのに苦労した．どの事業者が提供しているフリー Wi-Fi が該当するかも不明で，ウェブでもこの情報を得られなかったことから，それらしい商店街などを歩き回って探すしかなかった．結局，両市ともSSID=“#Passpoint” を吹いている基地局が該当のものであると判明した．その後，iPhone と Android の両端末で，手動で 1X 認証の設定を行い，自前のアカ
ウントで正常に認証・接続できることを確認した．

NGH hub オペレータにおいて，すべての .jp レル
ムを転送できるようにするため，正規表現に対応したシステムが開発中である．この開発が完了すれば，国内のすべての eduroam アカウントが利用できるようになる．また，日本以外の国でも，代表プロキシを設置すれば，容易に eduroam-NGH 連携が実現できるようになる．一方，現行の eduroam のルールでは，Hotspot 2.0 対応システムの実装方法や，eduroam 以
外の SSID の利用についての規定がないため，今後，GÉANT の eduroam コミュニティでも調整を進める
必要がある．

5 むすび
eduroam は世界各地において市街地や空港・駅，病院，博物館などに導入が進み，キャンパス無線 LANの枠を越えて，社会における教育・研究を支える基盤として認識されるようになってきた．一方，近年では公衆無線 LAN のセキュア化が進められるようになり，
既に一部のフリー Wi-Fi には 1X 認証や次世代ホットスポット (NGH) が導入されている．本報告では，セキュア化されたフリー Wi-Fi と eduroam を連携する

— 58 —

SENAC Vol. 51, No. 1（2018. 1）

ことで，教育・研究基盤がさらに強化される可能性を

2017/06/press20170627-01.html

示した．また，その実現に向けて参加した City Wi-Fi(2017 年 9 月 25 日参照)

Roaming トライアルの概要と，実証実験の成果を紹介した．

Hotspot 2.0 仕様による自動接続については，端末[10] Hideaki Goto, “Inter-roaming architecture foreduroam/govroam/public WLAN integrated
services,” The TNC17 Networking Conference,

側・基地局側とも，まだ実装上の不具合が多い．そのMobility Day, 2017.

改善動向を見ながら，eduroam-NGH 連携方式を洗練

https://tnc17.geant.org/core/event/23

し，実用化していくことが，今後の課題である．

(2017 年 9 月 25 日参照)

本研究の一部は，平成 29 年度国立情報学研究所公募型共同研究の助成を受けた．

参考文献[1] eduroam JP: https://www.eduroam.jp/(2017 年 9 月 25 日参照)
[2] 後藤英昭, 中村素典, 曽根秀昭, “キャンパス無線
eduroam の国内外の最新動向,” 大学 ICT 推進協
議会 2015 年度年次大会 論文集 3E3-2, 2015.

[3] 後藤英昭, 中村素典, 曽根秀昭, “キャンパス無線
eduroam と関連サービスの最新動向,” 大学 ICT推進協議会 2016 年度年次大会 論文集 WE25,

2016.
[4] “ライブドアと NII、学術無線 LAN ローミング
基 盤 の 共 同 実 験 を 開 始,” INTERNET Watch,2010.3.9.
http://internet.watch.impress.co.jp/
docs/news/353536.html(2017 年 9 月 25 日参照)
[5] LinkNYC: https://www.link.nyc/(2017 年 9 月 25 日参照)
[6] WiFi4EU | Free Wi-Fi for Europeans:
https://ec.europa.eu/

digital-single-market/en/policies/
wifi4eu-free-wi-fi-europeans(2017 年 9 月 25 日参照)
[7] セキュア公衆無線 LAN ローミング研究会 (NGHSIG): http://nghsig.jp/(2017 年 9 月 25 日参照)
[8] City Wi-Fi Roaming: http://worldwifiday.
com/city-wi-fi-roaming/(2017 年 9 月 25 日参照)
[9] “City Wi-Fi Roaming 大学としては世界初参加世界中でつながる次世代ホットスポット (NGH)と eduroam の連携を推進 -,” 東北大学プレスリ
リース, 2017.6.27.

https://www.tohoku.ac.jp/japanese/

— 59 —

SENAC Vol. 51, No. 1（2018. 1）

[報 告]
［報 告］

SC17 報 告東北大学サイバーサイエンスセンター・スーパーコンピューティング研究部
江川隆輔

2017 年 11 月 12 日〜17 日に米国コロラド州デンバー市のコロラドコンベンションセ
ン タ ー に お い て ， International Conference for High Performance Computing,
Networking, Storage and Analysis (SC17)が開催されました．SC は高性能計算とその関連技術に関する国際会議で，今年度で 27 回目を迎えます．
今年は世界各国から 12,753
人の研究者，システム開発者・ベンダーが集い，当該分野に関する最新の成果発表や技術紹介があり，将来のエクサフロップスレベルのシステムとアプリケーション開発に関して活発な議論がされました．サイバーサイエンスセンターは，本学の高性能計算に関する研究成果の発信を目的に，流体科学研究所，金属材料研究所，メディカルメガバン
クと連携し展示会に参加しました．展示ブースでは，本センターのシステム紹介，将来の高性能システムのための研究開発に関する展示・成果発表を行いました．また，JST
CREST の展示ブースに於いてもスーパーコンピューティング研究部の研究成果の一部を発表しました．例年同様，展示ブースに多くの訪問者を迎え，広く当センターの活動をアピールしただけでなく，大変有意義な議論と情報交換を行う事ができました．

東北大学展示ブース

会場にそびえ立つ熊— 60 —

SENAC Vol. 51, No. 1（2018. 1）[Web 版大規模科学計算システムニュースより]
大規模科学計算システムニュースに掲載された記事の一部を転載しています。 http://www.ss.cc.tohoku.ac.jp/tayori/

利用負担金額の表示コマンドについて (No.251)

本センター大規模科学計算システムでは、利用者の利用額とプロジェクトごとに集計した負担額、請求情報を表示するためのコマンドとして ukakin, pkakin があります。また、利用者のジャーナル
情報とプロジェクトごとに集計したジャーナル情報を CSV 形式で出力するコマンド ulist, plist があります。これらのコマンドは、並列コンピュータ(front.cc.tohoku.ac.jp)にログインして使用し
ます。
コマンド名

機

能

ukakin

利用者ごとの利用額を各システム、月ごとに表示

pkakin

プロジェクトごとに集計した負担額、請求情報を表示

ulist

利用者ごとのジャーナルを CSV 形式で出力

plist

プロジェクトごとに集計したジャーナルを CSV 形式で出力

いずれも、前日までご利用いただいた金額を表示します。コマンド使用例は大規模科学計算システ
ムウェブページをご覧ください。
負担金の確認http://www.ss.cc.tohoku.ac.jp/utilize/academic.html#負担金の確認(共同利用支援係)

— 61 —

SENAC Vol. 51, No. 1（2018. 1）― SENAC 執筆要項 ―１．お寄せいただきたい投稿内容
サイバーサイエンスセンターでは、研究者・技術者・学生等の方々からの原稿を募集しております。
以下の内容で募集しておりますので、皆さまのご投稿をお待ちしております。なお、一般投稿いただ
いた方には、謝礼として負担金の一部を免除いたします。
・一般利用者の方々が関心をもたれる事項に関する論説・センターの計算機を利用して行った研究論文の概要
・プログラミングの実例と解説・センターに対する意見、要望
・利用者相互の情報交換

２．執筆にあたってご注意いただく事項(1)原稿は横書きです。
(2)術語以外は、「常用漢字」を用い、かなは「現代かなづかい」を用いるものとします。
(3)学術あるいは技術に関する原稿の場合、200 字～400 字程度のアブストラクトをつけてください。
(4)参考文献は通し番号を付し末尾に一括記載し、本文中の該当箇所に引用番号を記入ください。
・雑誌：著者,タイトル,雑誌名,巻,号,ページ,発行年
・書籍：著者,書名,ページ,発行所,発行年
３．原稿の提出方法
原稿のファイル形式は Word を標準としますが、PDF での提出も可能です。サイズ*は以下を参照してください。ファイルは電子メールで提出してください。
－Word の場合－
・用紙サイズ：A4
・余白：上＝30mm 下＝25mm 左右＝25mm 綴じ代＝0
・標準の文字数（45 文字 47 行）＜文字サイズ等の目安＞
・表題＝ゴシック体 14pt 中央 ・副題＝明朝体 12pt 中央
・氏名＝明朝体 10.5pt 中央
・所属＝明朝体 10.5pt 中央
・本文＝明朝体 10.5pt
・章・見出し番号＝ゴシック体 11pt～12pt*余白サイズ、文字数、文字サイズは目安とお考えください。
４．その他(1)執筆者には、希望により本誌*（10 部以内の希望部数）と本誌 PDF 版を進呈します。
*2014 年末で、別刷の進呈は終了しました。
(2)一般投稿を頂いた方には謝礼として、負担金の一部を免除いたします。免除額は概ね 1 ページ
1 万円を目安とします。詳細は共同利用支援係までお問い合わせください。
(3)投稿予定の原稿が 15 ページを超す場合は共同利用支援係まで前もってご連絡ください。
(4)初回の校正は、執筆者が行って、誤植の防止をはかるものとします。
(5)原稿の提出先は次のとおりです。
東北大学サイバーサイエンスセンター内 情報部情報基盤課共同利用支援係
e-mail uketuke@cc.tohoku.ac.jp
TEL022-795-3406

— 62 —

SENAC Vol. 51, No. 1（2018. 1）

スタッフ便り新年，明けましておめでとうございます．サイバーサイエンスセンターがある青葉山にも冬が訪れ，時折白一面の雪化粧が見られるようになりました．そんな寒い中，サイバーサイエンスセンターのスーパーコンピュータは計算需要が逼迫するこの時期を乗り越えるべく，必死でがんばってくれています．本年も日本の科学技術・産業・社会を支える重要な基盤として，皆様の研究開発活動を支援できるよう取り組んでまいります．どうぞよろしくお願いいたします．
仙台市営地下鉄東西線が開通して早 2 年強経ちますが，「青葉山」駅に近いサイバーサイエンスセ
ンター周りも大きく変わってきています．サイバーサイエンスセンター向かいの青葉山新キャンパスに農学部・農学研究科などが移転したり，新しい建物が今も建設されています．人の流れも駅中心に変わり，駅からキャンパスバスが運行されています．青葉山キャンパスにお越しの際には，ぜひサイ
バーサイエンスセンターにもお立ち寄りください．スーパーコンピュータなどの見学が随時できます．
皆様のお越しをお待ちしております．(K.K)
青葉山キャンパスには，平成 10 年 4 月から 2 年間，薬学部に勤務していたことがあります．その頃の私は，体重も増えてきたので，昼には薬学部内にあったテニスコートでテニスをしたり，青葉山
キャンパス内を歩いたりしていたように記憶しています．行先は，宮教大の向こうのゴルフ練習場
だったり，時には薬用植物園の中だったり，いろいろ足を延ばしていました．当時はクマの出没など全然なかったので，何の心配もなく季節を楽しみながら歩いていました．
平成 29 年 7 月に異動してくるとキャンパスの様子は大きく様変わりをしていました．サイバー
サイエンスセンターの屋上から見る景色は新キャンパスを一望し，いくつもの新しい建物が立ち並んでいました．建物内に入ってみると大きな空間があったり，凝ったデザインのテーブルや椅子が置かれていたりしています．
地下鉄も便利で，東京へ出張しようとすると片平キャンパスから仙台駅に歩くよりも，青葉山駅から地下鉄を利用した方が早く，汗もかかずに仙台駅に着いてしまいます．
ただ，ちょっとだけ残念だなあと感じていることがあります．青葉山コモンズ（農学分館）は，す
ごく立派なんだけれども，もっと学生さんが利用してくれればいいのに～．就職した後に始めたスポーツですが，若いころはよく野球をやっていたので，やっぱりグランドがあるといいなあ～．(T.O)

情報部情報基盤課スタッフ採用のお知らせ[採用]
2017.12.18

畑山詩乃（総務係・派遣職員）

SENAC 編集部会
滝沢寛之 水木敬明 後藤英昭 江川隆輔
伊藤昭彦 吉田貴子 大泉健治 小野 敏
斉藤くみ子

編集・発行

印
開発中の青葉山キャンパス

刷

平成 30 年 1 月発行
東北大学
サイバーサイエンスセンター
仙台市青葉区荒巻字青葉 6-3
郵便番号 980-8578
東北大学生活協同組合
プリントコープ

索
SENAC Vol.50 (No.1～No.4)
項 目引
2017.1～2017.10

巻頭言
気象予測とコンピュータ
岩崎俊樹
スーパーコンピューティング研究部教授に着任して滝沢寛之
共同研究成果DNA鎖切断の反応動力学シミュレーションとAVS/Expressを用いた三次元可視化
菱沼直樹他
ベクトル型スーパーコンピュータを用いた大規模リフレクトアレーの高利得化
今野佳祐他
アルゴン熱プラズマジェットによる鉄ナノ粒子の大量創製プロセスの３次元シミュレーション
茂田 正哉
数値シミュレーションから見出された余効すべり
伝播速度と摩擦特性の関係
有吉慶介他
有機-金属界面での電子準位接続の精密決定に向けた GW space-time コードの高速化
柳澤将
民間航空機開発における大規模 CFD 解析の適用
前田一郎
周期構造の電磁界散乱特性解析
有馬卓司他
次世代低騒音タイヤ開発に向けた高精度流体音響
解析
藤井孝藏他

No.

p.

項 目

1

1

報告
SC16 報告1＜計算科学・計算機科学人材育成のためのスーパーコンピュータ無償提供利用報告＞ 東北大学大

2

1

3

1

13

1

19

2

3

3

1

3
4

8
1

4

6

研究成果
HPGMG-FV を用いた SX-ACE の性能評価
江川隆輔他
大規模科学計算システム
SSH アクセス認証鍵生成サーバの利用方法
アプリケーションサービスの紹介
三次元可視化システムの紹介
非線形構造解析プログラム MSC.Marc の紹介

3

2
2
2
4

25

10
17
36
14

1

25

1

31

1

37

1

41

2

5

3

19

紹介
サイバーサイエンスセンター高性能計算技術開発（NEC）共同研究部門のご紹介

2

40

退職のご挨拶
退職にあたって

2

53

お知らせ
平成 29 年度講習会サイバーサイエンスセンター
講習会のご案内サイバーサイエンスセンター講習会(夏季、秋季開催分)のご案内

千葉実

学院データサイエンスプログラムにおけるプ
ログラミングトレーニング演習第 13 回シナジー研究会報告瀨川悦生
第25回高性能シミュレーションに関するワークシ
ョップ(WSSP)を開催しました小林広明
サイバーサイエンスセンターセミナー報告

p.

1

45

2

44

2
2

46
48

3

23

3

28

3

31

3

33

4

29

4

30

4

31

4

35

3

21

1
2
2
2

45
55
55
56

2
2

56
58

2

58

2
2
2

59
60
62

2
2

63
63

2
3
3
4

64
34
34
36

＜計算科学・計算機科学人材育成のためのスーパーコンピュータ無償提供利用報告＞ 情報科学研究

科 超高速情報処理論利用報告
滝沢寛之＜プレスリリース報告＞東北大・大阪大・NEC・国際航業・エイツー世界初、地震発生から 30 分以内にスーパーコンピュータを用いて津波浸水被害を推計するシステムが内閣府「津波浸水被害推定システム」として採用
City Wi-Fi Roaming 大学としては世界初参加- 世界中でつながる次世代ホットスポット(NGH)とeduroam の連携を推進 ＜成果報告＞日本気象協会が本センタースーパーコ
ンピュータ（NEC SX-ACE）を活用して「熱中症ゼロへ」プロジェクト『熱中症セルフチェック』を新たに開発
小林広明センター長特別補佐が「情報化促進貢献個人等表彰」文部科学大臣賞を受賞しました
JHPCN 学際大規模情報基盤共同利用・共同研究
拠点第 9 回シンポジウム報告
江川隆輔＜岩手大学グループ利用制度を利用して＞ 高速計算サービスの機関利用
中西貴裕
サイバーサイエンスセンターオープンキャンパ
ス報告
小松一彦
利用相談室便り
平成 29 度の利用相談について

大学 ICT 推進協議会 2016 年度 年次大会論文集より転載『銅酸化物の有効モデルに対する揺らぎ交換近似』コードの SX-ACE 向け最適化
山下毅他
スーパーコンピュータ SX-ACE の省エネルギー運用
森谷友映他
キャンパス無線 eduroam と関連サービスの最新動向
後藤英昭他
東北大学自営光ファイバーの経年劣化調査
七尾晶士他

小松一彦

No.

「Web 版大規模科学計算システムニュース」より
コンパイラのバージョンアップについて（No.232）
利用負担金額の表示コマンドについて(No.234)
研究成果リスト提出のお願い（No.235）
コンパイラのバージョンアップについて（No.238）
バッチリクエストの最大経過時間の指定について（No.239）
Gaussian16 の提供開始について（No.239）
数式処理プログラム Mathematica のバージョン
アップについて（No.239）
科学技術計算言語MATLABのバージョンアップについて（No.239）
平成 29 年度利用負担金について（No.239）
平成 29 年度共同研究について（No.240）
計算科学・計算機科学人材育成のためのスーパ
ーコンピュータ無償提供制度について（No.240）
民間企業利用サービスについて （No.240）
大規模科学計算システムの機関（部局）単位での利用について（No.240）
利用負担金額の表示コマンドについて(No.243)
コンパイラのバージョンアップについて（No.244）
利用負担金額の表示コマンドについて(No.248)ࢩࢫࢸ୍࣒ぴィ⟬ᶵࢩࢫࢸ࣒ᶵ ✀ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ

6;$&(୪ิࢥࣥࣆ࣮ࣗࢱ

/;5H

ࢧ࣮ࣂ࡜࣍ࢫࢺྡ
ࣇࣟࣥࢺ࢚ࣥࢻࢧ࣮ࣂIURQWFFWRKRNXDFMS

66+ ࢔ࢡࢭࢫㄆド㘽⏕ᡂࢧ࣮ࣂ

NH\FFWRKRNXDFMSࢧ࣮ࣅࢫ᫬㛫฼⏝ࢩࢫࢸ࣒ྡ➼฼⏝᫬㛫ᖏ

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ

㐃

⥆

㐠

㌿

୪ิࢥࣥࣆ࣮ࣗࢱ

㐃

⥆

㐠

㌿

ࢧ ࣮ ࣂ

㐃

⥆

㐠

㌿

ྍど໬ᶵჾᐊ

ᖹ᪥ 㹼

㤋ෆ฼⏝ᖹ᪥ 㹼

ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢱ㸦6;$&(㸧ࡢ฼⏝ᙧែ࡜ไ㝈್฼⏝ᙧែ฼⏝ࣀ࣮ࢻᩘ ᐇ⾜᫬㛫ไ㝈ͤ㸦⤒㐣᫬㛫㸧
つᐃ್㸸 㐌㛫㏻ᖖࠥ↓ᩱ ᫬㛫

ࠥ ᫬㛫

ࠥ ᫬㛫

ࢹࣂࢵࢢ᭱኱್㸸 ࣨ᭶࣓ࣔࣜࢧ࢖ࢬไ㝈 T ࢜ࣉࢩࣙࣥ E ࢜ࣉࢩࣙࣥ

*%ࣀ࣮ࢻᩘ

฼⏝ࣀ࣮ࢻᩘ

V[*%
*%ࣀ࣮ࢻᩘI
GHEXJ

฼⏝ࣀ࣮ࢻᩘ

ͤࣀ࣮ࢻ௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣ03,ࡢ฼⏝ࡀᚲ⏝୪ิࢥࣥࣆ࣮ࣗࢱ㸦/;5H㸧ࡢ฼⏝ᙧែ࡜ไ㝈್฼⏝ᙧែ฼⏝ࣀ࣮ࢻᩘ
ͤ

㏻ᖖ

ࠥ

࢔ࣉࣜ
ࢣ࣮ࢩࣙࣥᐇ⾜᫬㛫ไ㝈㸦⤒㐣᫬㛫㸧つᐃ್㸸 㐌㛫᭱኱್㸸 ࣨ᭶࡞ࡋ࣓ࣔࣜࢧ࢖ࢬไ㝈
*%ࣀ࣮ࢻᩘ

㸦 ࢥ࢔ࡲ࡛㸧฼⏝ࣀ࣮ࢻᩘ
O[*% ᫬㛫
*%
㸦&38 ᫬㛫ྜィ㸧ͤࣀ࣮ࢻ௨ୖࢆ฼⏝ࡋࡓ୪ิᐇ⾜࡟ࡣ03,ࡢ฼⏝ࡀᚲ⏝఍ヰᆺT ࢜ࣉࢩࣙࣥ E ࢜ࣉࢩࣙࣥ

D