巻頭言

Vol.16, No.1 号の発刊に当たって京都大学学術情報メディアセンター

牛島 省

本号では、京都大学学術情報メディアセンターにおいて実施された平成 28 年度の共同「若手・女性研究者奨励枠」において研究報告について特集いたします．平成 28 年度は，27 件，
「プログラム高度化共同研究」において 2 件の共同研究が実施されました．
学術情報メディアセンターでは，スーパーコンピュータ共同研究制度として，「若手・女性研究者奨励枠」，「大規模計算支援枠」，「プログラム高度化共同研究」という３種類の研究支
援事業を実施しています．課題募集は年度開始前後に行われ，追加募集も実施される場合があります．詳細は，http://www.iimc.kyoto-u.ac.jp/ja/services/comp/support/kyodo/をご覧下さい．応募された課題の採否は，スーパーコンピュータ共同研究企画委員会にて審査されます．
これらの共同研究制度のうち，「若手・女性研究者奨励枠」の応募資格は，(1)40 歳未満の若手研究者 (学生を含む, 性別は問わない)，あるいは，(2) 女性研究者 (年齢は問わない)，と平成 27 年度より拡大されており，スーパーコンピュータを利用することで学術的にインパクトがある成果を創出できる課題に対して，計算機利用負担金の全額または一部を本センターが負担しています．さらに，「若手・女性研究者奨励枠」は，平成 28 年度より JHPCN (学際大規模情報基盤共同利用・共同研究拠点) の活動の一環となり，優れた
課題は，JHPCN の萌芽型共同研究へ推薦されることになりました．推薦された課題は，
JHPCN シンポジウムにてポスター発表を行うことが可能となるなどの特典がありますので，この共同研究制度をご活用ください．本号に掲載された 27 件の「若手・女性研究者
奨励枠」の研究成果報告では，各種の流体および流体構造連成問題，破壊力学，波動問題，
有機分子材料，量子化学，分子動力学計算，MRI データ解析，力学系軌道問題，カオス力
学系，メモリ一貫性モデルなど多様な研究が行われ，全国 15 大学にわたる若手研究者との共同研究となりました．
また，「プログラム高度化共同研究」では，京都大学のスーパーコンピュータをグループ
コースまたは専用クラスタコースでご利用の研究グループ (JHPCN にて京都大学のスー
パーコンピュータを利用する研究者も含む) を対象に，プログラムの制御構造・データ構造の改良による性能チューニングや並列化手法の改良，問題分割・負荷分散方式などの改良による大規模計算プログラムの高度化・高性能化を支援しています．平成 28 年度は，地盤の変形・破壊を扱う有限要素解析と，津波浮遊物が関係する流体構造解析のプログラム高度化支援を行いました．一般固有値解析における数値計算ライブラリの利用方法やハイブ
リッド並列計算法，また DEM と FEM をカップリングした解法に対する並列化など，様々な高度化・高性能化が行われており，計算効率を向上させる際の参考になると思います．
今後も皆様の研究，教育にご活用いただけるようにセンター教職員も尽力していきます
ので，ご利用・ご支援のほど、よろしくお願い申し上げます．

1
京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 1

スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告

種々の波動問題に対する積分核に依存しない演算子積分時間領域高速多重極境界要素法の開発
斎藤 隆泰 ∗∗∗

1

群馬大学 大学院理工学府 環境創生部門

はじめに法は, Lubich により提案された, 畳み込み積分を精度良く, 安定に計算する方法である. 本研究では, こ

本研究では, 種々の波動問題に対する積分核に依存の演算子積分法を時間領域境界積分方程式に現れるしない演算子積分時間領域高速多重極境界要素法を畳み込み積分に適用することで, 時間領域境界要素

開発することを行う.

法における, 時間増分が小さい場合の数値不安定性

境界要素法 (BEM: Boundary Element Method)の問題を解決することを行った. 演算子積分法を用

は, 波動問題に対する有効な解析手法として発展していた場合, 境界積分方程式における積分核は, 時間領

きた. しかしながら, その時間領域解法には, いくつ

域ではなく Laplace 変換領域でのものとなる. 通常,か欠点が存在する. まず, 一般的に時間領域境界要素

例えば粘弾性波動問題や, 飽和多孔質弾性波動問題の法は, 陰解法であるため, 大規模問題の解析は難しい

ように弾性波速度が周波数に依存するような問題は,点が挙げられる. すなわち, 計算効率の問題である.

時間領域で閉じた形の基本解を求めることはできな

また, 時間増分が小さい場合に数値解が不安定になる

いが, Laplace 変換領域の基本解は容易に求めること

ことも欠点として知られている. さらに, 時間領域境ができる. すなわち, 演算子積分法の利用により, 時

界要素法は, 基本解を用いて計算をする必要があるた

間増分が小さい場合の数値不安定性の問題と, 時間領

め, 基本解を閉じた形で求めることができない問題に

域境界要素法の適用範囲の問題を改善することがで対しては, 通常, 適用することができない. このこと

きる.

は, 時間領域境界要素法は限られた問題しか適用できない適用範囲の問題を持っていることとなる. このような欠点が原因となり, 時間領域境界要素法を用いた大規模波動解析はこれまでほとんど行われていな
いのが現状である. そこで, 本研究では, これら 3 つ

3

積分核に依存しない高速多重極法の適用の欠点を克服し, 境界要素法の適用範囲を拡大した新

前節で述べた演算子積分法を適用したとしても, 結しい時間領域境界要素法の開発を目指す. 以下では,局の所, 異方性弾性波動問題のように Laplace 変換領

これら 3 つの欠点についてどのような方法を用いて域で閉じた基本解を得ることができないような問題対応したかについてまとめ, 簡単な数値解析例を示もある. また, 演算子積分法を適用したとしても, 時

すと共に, 今後の課題について述べる.

間領域境界要素法の欠点である計算効率の問題を克服することにはならない. そこで, 本研究では, これ

2演算子積分法の適用
まず, 本研究では, 時間領域境界要素法の時間方向の離散化に演算子積分法 (CQM: Convolution Quadra-ture Method)[1] を適用することを行う. 演算子積分

2 − Academic Center for Computing and Media Studies, Kyoto University

らに対応するために, 積分核に依存しない高速多重極法 (KIFMM: Kernel Independent Fast Multipole

Method) を適用することを考える. 積分核に依存しない高速多重極法は, Ying ら [2] より提案された比較的新しい高速多重極法であり, 擬似的な多重極展開係数を作り出した後, 通常の高速多重極アルゴリズ

incident wave

incident wave

図 2: 空洞群周辺の弾性波動場.

図 1: 解析モデル.

ムの下で計算を高速化する方法である. その際, 積分
核の多重極展開を直接に必要としないため, 先に述べたように異方性弾性波動問題のような閉じた基本解が求まらず, 多重極展開が難しい問題に対しても高速多重極法の適用が可能であると考えられる. また,その際, 従来の高速多重極法で必要となった変形ベッセル関数等の複雑な関数を取り扱う必要も生じない

容量の制限により解析することが難しい問題であるが, 本手法を用いることにより, 単一ノードの計算機でも解析が可能となっている. 図 2 は, 無次元化時刻

ct/a = 72.0 における解析結果であり, 入射波の波面が丁度左端の空洞に当る時刻を ct/a = 0.0 としている. 入射波の波面が空洞群を通過し, それに伴う多重散乱波が生じている様子を見て取ることができる.

ことから, 本研究が将来有効な大規模波動解析手法と成り得る可能性を秘めている.

5

まとめと今後の課題
積分核に依存しない時間領域高速多重極境界要素

4

数値解析例

法を開発した. 今後は超音波非破壊評価法に関する

以下, 簡単な数値解析例を示す. 社会基盤構造物の種々の弾性波動問題へ適用する予定である.

維持管理において, 構造物を壊すことなくそれらの健全度を評価し (非破壊評価法), 健全度に応じて適切な謝辞

補修・補強等を施すことは重要である. その非破壊評
価法において, 現場で最も広く利用されている方法は

本研究は平成 28 年度スーパーコンピュータ共同研

超音波非破壊評価法であるが, 超音波は固体中では弾究制度（若手・女性研究者奨励枠）の支援により行性波の性質を示すため, 超音波非破壊評価法の数値シ

われ, 現在も, 平成 29 年度学際大規模情報基盤共同ミュレーションは, 弾性波動シミュレーションへと帰利用・共同研究拠点公募型共同研究 (jh170045, 研究着される. そこで, ここでは, 種々の波動問題の中で最代表者：斎藤隆泰), にて引き続き京都大学の支援をも基礎的はスカラー波動問題 (弾性波動問題における受けている.

SH 波動問題) に対する本手法による解析結果を示す.
図 1 は 2 次元無限弾性体中の半径 a の 500 個の空洞群による入射スカラー波の散乱問題の解析モデルを示している. 入射波は x1 軸の正の向きに波速 c = 1で伝搬する平面波とし, 各空洞は, 72 境界要素でモデル化されている. 全時間ステップ数は 128 であることから, この問題の自由度は 4608000(72×500×128)である. 通常の時間領域境界要素法では, 当然, 記憶

参考文献[1] Lubich, C.: Convolution quadrature and discretized operational calculus I, Numer. Math., 52,
pp. 129-145, (1998).
[2] Ying, L., Biros, G. and Zorin, D.: A kernelindependent adaptive fast multipole algorithm in
two and three dimensions, J.Comput. Phys., 196,
pp. 591-626, (2004).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 3

光機能性有機分子材料の構造と電子構造の解明による有機系太陽電池の素子性能向上
梅山有和
京都大学大学院工学研究科分子工学専攻

1 緒言(1/1)溶液中の吸収スペクトルを Figure 2 に示す。
ZnPTD の吸収スペクトルは Soret 帯、Q 帯ともに

ポルフィリンは 400~450 nm に Soret 帯と呼ば

LDD1 よりも長波長シフト・ブロード化しており、れる強い吸収と 550~600 nm に Q 帯と呼ばれる

可視・紫外領域で光捕集能が改善されていること

中程度の吸収をもち、色素増感太陽電池の増感色がわかった。

素として有望である。
最近 Diau らは 800nm 程度
までの光捕集能をもつ LDD1 というポルフィリ
ンダイマーを用いることで変換効率η = 8.8%を報告している。本研究では、この優れた光捕集能をもつポルフィリンダイマーのさらなる光捕集能向上を目指し、ジアリールアミノ基を meso 位に導入したポルフィリンダイマーZnPTD を設計した(Figure 1)。

Figure 2. UV-vis absorption spectra of ZnPTD (solidline) and LDD1 (dashed line) in ethanol/toluene (v/v =1/1).

色素増感太陽電池の作製と評価を行ったところ、
ZnPTD において JSC = 16.0 mA cm–2, VOC =0.626 V, FF= 0.615, η = 6.2 %、LDD1 において

JSC = 19.0 mA cm–2, VOC = 0.633 V, FF= 0.680,
η = 8.2 %という結果が得られた。光捕集能の改善
により、820nm 以降においては ZnPTD の方がFigure 1. Molecular structures of porphyrins.

LDD1 よりも高い IPCE 値(外部量子収率)を示している。しかし、より短波長の領域で ZnPTD の2 結果と考察方が低い IPCE 値となっているために、結果とし2.1 新規ポルフィリン色素の光吸収特性と色素増感太陽電池への応用

この理由としては、2,6-ビス（ドデシルオキシ）ZnPTD および LDD1 のエタノール/トルエン

4 − Academic Center for Computing and Media Studies, Kyoto Universityて LDD1 よりも小さい JSC 値になったといえる。
フェニル基の数がLDD1 では4 つであるのに対して、ZnPTD では 2 つしかないためダイマー間の会合の影響が長鎖アルコキシ基の少ない ZnPTD致する。また、吸着部位の LUMO の電子密度はで大きくなっていることが考えられる。

色素の励起状態と酸化チタンの 3d 軌道の電子カップリングに影響を与えるため、ZnPTD でも

2.2 理論計算

LDD1 同様、吸着部位の LUMO の電子密度が高い

ポルフィリンの最安定化構造およびそのフロン
ティア軌道の電子構造について知見を得るために、
密度汎関数法（DFT）による理論計算（B3LYP）を行った。基底関数系に 6-31G(d)を選択し、Gaussian 09 プログラムを用いた。Figure 3 に
ZnPTDT および LDD1 の最高被占軌道（HOMO）と最低空軌道（LUMO）における電子密度分布を示す。ZnPTD および LDD1 の HOMO-LUMO ギャ
ップを比較すると ZnPTD (1.84 eV) < LDD1 (1.91
eV)となった。ZnPTD では LDD1 と比較して小さな HOMO-LUMO ギャップとなっており吸収の長波長シフトが予想されるが、これは実験結果と一ことから酸化チタンへの良好な吸着と電子注入が期待できる。一方、HOMO の電子密度は吸着部位から遠いアミノ基に主に広がっており、電荷再結合の抑制も期待される。このように、本計算結果
は実験結果の理論的解釈の一助となった点で意義がある。

3 参考文献[1] H. Imahori, T. Umeyama and S. Ito, “Large

π-Aromatic Molecules as Potential Sensitizers in

Dye-Sensitized Solar Cells”, Acc. Chem. Res., 42,
1809-1818 (2009).

[2] H. Imahori, T. Umeyama, K. Kurotobi and Y.
Takano,“Self-Assembling

Porphyrins

andPhthalocyanines for Photoinduced Charge Separation

and Charge Transport”, Chem. Commun., 48,
4032-4045 (2012).

[3] T. Higashino, K. Kawamoto, K. Sugiura, Y.

Fujimori, Y. Tsuji, K. Kurotobi, S. Ito and H. Imahori,“A Effects of Bulky Substituents of Push-Pull

Porphyrins on Photovoltaic Properties of Dye-

Sensitized Solar Cells”, ACS Appl. Mater. Interfaces, 8,15379-15390 (2016).

[4] T. Higashino, K. Sugiura, Y. Tsuji, S. Nimura, S.
Ito and H. Imahori, “A Push-Pull Porphyrin Dimer with

Multiple Electron-Donating Groups for Dye-SensitizedSolar Cells: Excellent Light-Harvesting in NearInfrared Region”, Chem. Lett., 45, 1126-1128 (2016).

Figure 3. Selected molecular orbital diagrams for (a)ZnPTD and (b) LDD1 obtained by DFT calculations
with B3LYP/6-31G(d). To simplify the calculations,alkyl chains on the diarylamino groups were replacedwith methyl ones, and dodecyloxy groups on thephenyl groups were replaced with methoxy ones.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 5

乱流相分離制御による自己組織化構造の能動的選択
岡本幸也
大阪大学大学院基礎工学研究科

本研究では、流動、伝熱、物質移動の移動現象と相分離現象とを相互に作用させ、熱交換器の伝熱効率を、相分離流体を利用して促進する手法開発を目指し、強制対流伝熱下における相分離現象の二次元数値解析を行った。結果から、相分離が生じた壁、出口付近で、相分離に起因する速度変動が生じ、局所的に伝熱が促進されることが観察出来た。以上から、あらかじめシミュレーションによって伝熱が促進される位置を知ることは熱交換器設計において重要であると考えら
れる。

1

緒言

2

解析手法

高温で完全に混じり合っている多相流体は、急計算対象とした系は、実験[1]で用いられたアセ

冷などによって熱力学的に不安定な状態となるとトンとヘキサデカンからなる非圧縮性の二成分混

各成分の相へと分離する。この相分離現象はスピ合流体であり、初期状態として、完全に混合されノーダル分解と呼ばれ、特に、相分離の際にネッた状態(成分比率が、50%-50%)を想定した。相分

トワーク状の自己組織化構造が形成されることが離現象には、二成分の界面をオーダーパラメータ

知られている。加えて、局所的な物性の変化や、によって拡散的に再現する手法であるフェーズフ

Korteweg 効果に起因する速度変動が知られておィールド法を用いた。本研究では、オーダーパラ

り、これらの特徴から、石油の増進回収法や、多メータとして、アセトンのモル分率 xac を用いた。

孔質材料形成、食料品生産まで幅広いプロセスでこのオーダーパラメータを用いて相分離挙動を表の応用が期待されており、自己組織化構造形成と現する時間発展方程式は以下の Cahn-Hilliard 式、流れ場との相互作用の理解が重要となっている。
近年、相分離の工学的応用として、熱交換器の作動流体として、相分離流体を利用し、伝熱効率を促進する研究がなされている[1]。しかし、バルク










1


 

の温度変化の測定が主であり、相分離流体による
局所的な伝熱促進効果の解明はなされていない。

となる。ここで、t は時刻、xj は座標、ui は速度、また、実用化のためには、実験を行う前段階としRe は Reynolds 数、Sc は Schmidt 数である。化学て、シミュレーションを行い、最適な溶液選択やポテンシャルは、正則溶液モデルを拡張した以下

設計指針を得ることで、効率の良いプロセス開発の One-parameter Margules Model[1]を用いた。

が期待される。そこで本研究では、強制対流伝熱下における相分離現象の二次元数値計算を行い、相分離による伝熱促進効果の局所的なメカニズムについて検討を行った。

  log12

1 2    2 

ここで、ξは界面幅であり、分子動力学シミュ

6 − Academic Center for Computing and Media Studies, Kyoto Universityレーション[2]から、1 nm ほどであると考えられて場合の局所 Nusselt(Nu)数の分布を Fig. 2 に示す。

いるが、本研究では、格子解像度の制限から、10Kays et al. によって求められている、一様壁温のcm とした。

円管内層流熱伝達の局所 Nu 数の理論値も合わせ-3

流体場の支配方程式は、以下の連続式とNavier-Stokes 式、て示す[3]。横軸は、温度境界層の解析に用いられる流れ方向の無次元長さであるグレツ数とした。
結果から、相分離が無い場合の局所 Nu 数は、
円管内層流伝熱の局所 Nu 数の理論値である約 3.7


 0,


に漸近していることが分かる。これは、設定した計算条件において、十分温度境界層が発達してい







   1 

2

ることを示している。相分離がある場合の局所 Nu 

数は、出口付近で約 1.09 倍増大しており、相分離によって伝熱が促進されたことがわかる。

2 

2 2
   

相分離による伝熱促進効果の局所的なメカニズムを検討するため、時間 t = 0.2 s における壁垂である。ここで、P は圧力、ξn は、界面幅

直方向の速度成分の分布を Fig. 3 に、濃度分布を(Navier-Stokes 式中では 1 nm に設定)、R は気体定

Fig. 4 に示す。結果から、壁、出口付近に速度変数、Tc は相分離臨界温度(アセトン、 ヘキサデカ動が分布しており、Fig. 2 の局所 Nu 数が増大してン 系 で は 301 K) 、 M は モ ル 質 量 で あ る 。

いる領域と一致していることから、この速度変動

Navier-Stokes 式中の右辺第三項は、ポテンシャル

が伝熱を促進したと考えられる。一方、濃度分布の勾配を駆動力として流体が受ける力を表現してから、単相状態で流入した混合流体は、壁から相おり、相分離と流体場をカップリングしている分離が開始し、出口に向かって成長する様子が観

対流伝熱の支配方程式は以下のエネルギー式、察された。これは、Fig. 3 で速度変動が分布している領域と一致しており、相分離によって速度変








4



1

1 




動が生じ、伝熱を促進したと考えられる。特に、

 

相分離の成長している出口付近での速度変動が大

きく、相分離の成長段階に応じて、伝熱促進効果が変化していると言える。

4

結言である。ここで、T は温度、Pr はプラントル数、
Cp は定圧比熱である。右辺第二項は相分離に伴っ

強制対流伝熱下における相分離現象の二次元数て生じる混合エンタルピーの温度場への寄与を表

値解析を行い、計算結果から、相分離を伴う多成している。時間積分には、一次精度オイラー陽解

分混合流体は、熱交換器内で相分離を生じる際に、法を用いた。

速度変動を生じ、伝熱を促進することが分かった。

計算領域は、直径 1 mm、 長さ 5 mm の二次元

また、促進効果の強さは、相分離の成長段階に応チャンネルとし、境界条件も合わせて Fig. 1 に示じており、実際の熱交換器の設計の際には、伝熱す。格子点数は、壁垂直方向に 100 点、流れ方向が強く促進される位置を知ることが重要であるとに 500 点とした。計算条件は、Re = 0.8、Pr = 13.8,

言える。

Sc = 200 と仮定して行った。

3

結果と考察
相分離による強制対流伝熱への影響を検討する

参考文献

ため、相分離がある場合(Re = 0.8, Sc = 200)と無い

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 7

1. Ulmann A., Maevski K. and Brauner N.,“Enhancement of forced and freeconvection heat transfer rates by inducing
liquid-liquid phase separation of apar tial ly- mis cible eq ual- den sity bin ary
liquid system”, Int. J. Heat Mass.
Tr a n s f e r, 70 , pp . 36 3 - 37 7 (2 0 1 4 ) .
2. Gyula I. T. and Bjon K., “Phase field

Fig. 2 Distribution of local Nu at t = 0.2 s, Re = 0.8, Sc

modeling of spinodal decomposition in the= 200: Red line is with phase separation, Blue line isoil/water/asphaltene system”., Phys.

without phase separation and Black line is theoreticalvalue for tube.

Chem. Chem. Phys., 17, 20259 (2015).
3. W. M. Kays, Convective Heat and Masstransfer, McGRAW-HILL BOOKCOMPANY (1966).

5 mm= Random= (y)= 0 gradient

1 mm
=

wall,

= No-flux,= 0,

Fig. 1 Schematic of calculation domain and boundary conditions.

Fig. 3 Snapshot of concentration field at t = 0.2 s.

Fig. 4 Snapshot of vertical velocity field at t = 0.2 s.

8 − Academic Center for Computing and Media Studies, Kyoto University= 0 gradient= Constant flux
= Constant flux
= 0 gradient=0Singlet Fission 相互作用の量子化学計算
鬼頭(西岡)宏任、佐藤竜馬、重田育照
筑波大学 計算科学研究センター

1 INTRODUCTION

間自己無撞着(CASSCF)法に基づき、SF に関与する電子配置(図 1)を適切に考慮して、SF 相互作用

Singlet Fission(SF: 一重項励起子分裂)は、光を従来よりも簡便に計算する方法を示す。[3]吸収によって一重項励起状態(S1)になった分子が、隣接する基底状態(S0)の分子にその励起エネルギーの一部を無輻射過程で移動することによって、三重項励起状態(T1)の２つの分子が生成される現象である[1]。SF 自体は、1960 年代には既に、テトラセンやペンタセンなどのアセン分子結晶で起こることが観測されていた。しかし、SF を活用することで、単一接合太陽電池の持つ約 34％の太陽エネルギー利用効率限界(Shockley-Queisser
limit)を約 46％まで改善する可能性が理論的に示されて以降[2]、ここ 10 年間で注目が集まり、盛んに研究が行われるようになった[1]。

図 1

SF は通常、以下の過程で説明される。

S0 + S1 ↔ 1 (T1T1 ) → T1 + T1
光吸収直後の S1 が、隣接した２つの分子が作考慮する電子配置と SF 経路

2 METHODS

る”correlated triplet pair” 1(T1T1)へ、スピン許容以下の手順による 2 段階のユニタリー変換を用

遷移で変換し、その後スピン・デコヒーレンスに

いて SF を記述するための(準)透熱基底ハミルト

よって２つの三重項励起子が生成される。S0+S1

ニアンを作成する[3]。

から 1(T1T1)への変換は、分子結晶のモフォロジー

A)

分子結晶から、二量体モデルを作成に依存して、数十フェムト秒から数百ピコ秒と大B)

State–
二量体のGS からS5 までの6 状態を、
averaged(SA-)CASSCF(4,4)/6-31G(d)計算

きく異なった反応時間が観測される。
本研究では、SF 反応時間を決定する因子である S0+S1 と 1(T1T1)間の電子的相互作用(SF 相互作用)を求める、新しい計算手法を開発した[3]。通C)

手順(B)の最適化活性軌道をユニタリー変換
し、フラグメント局在軌道を作成

D)

手 順 (C) の 局 在 活 性 化 軌 道 を 用 い て常、材料系の有機分子励起状態計算では、計算コ

CAS-CI(4,4)を実行し、配置間相互作用(CI)

ストと精度のバランスから、一電子励起配置間相

係数 C A とエネルギーを保存

互作用(CIS)法や時間依存密度汎関数(TDDFT)法が用いられる。しかし、CIS や TDDFT は一電子励起状態を記述する手法であり、1(T1T1)のような多電子励起は考慮できない。そこで、完全活性空E)

手順(C)の局在活性化軌道を用いて S0S1 や
S1S0、電荷移動励起子 AC, CA、1(T1T1)などの波動関数を ORMAS-CI 法[4]から求め、その CI 係数 C D,ref を保存

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 9F)G)

特異値分解から C D,ref を再現するような C Aの実ユニタリー変換行列 U を決定
手順(D)のエネルギーを対角成分に持つ(準)は、良く一致していることが分かる。SF-CIS 計

断熱基底ハミルトニアンを行列 U からユニの一部が無視されている。また SF-CIS 計算は、

タリー変換して、透熱基底ハミルトニアン
H D を計算
上記の手順で作成された H D を基に、Dexter 機
構によるダイレクトな経路(図 1 青矢印)と、超交換機構による AC と CA をヴァーチャルに使った経路(図 1 赤矢印)を考慮し、式(1)を用いて SF 相
互作用計算を実行した。
V = S1S0 Ĥ D T1T1 −−D

スピン汚染の影響も含んでいる。一方、我々のCASSCF 計算は動的電子相関が考慮されないが、SF-CIS ではその一部が取り込まれている。これ
らの理由から、ペンタセンや DTT に対して、計算結果の違いが生じたと考えられる。
原因の一つに挙げた動的相関は H D の対角成分

D

差を取る際にそのエラーは相殺され、SF 相互作

CA Ĥ CA − Etun

用に与える影響は小さくなると予想される。そこ

D

S1S0 Ĥ AC AC Ĥ T1T1
D

AC Ĥ D AC − Etun

(1)

ここで Etun は、S0S1 と 1(T1T1)のエネルギー(H Dの S0S1 と 1(T1T1)の対角成分)平均値に設定する。
CI 計 算 は 、 GAMESS に 実 装 さ れ て い る
determinant CI コードを用いた。

で、Van Voorhis らが同様の二量体モデルに対して制約 DFT(CDFT)-CI 法を用いた結果[6]と比較する。CDFT-CI 計算では、ペンタセンに対しては 84meV、DTT に対しては 6.2meV の SF 相互作用を持つ[6]。計算条件の違いはあるが、動的相関が十分に考慮された CDFT-CI 計算が、Hsu らに較べて、我々の結果と良く一致している。この3 RESULTS & DISCUSSION
Hsu らは、Spin-flip CIS(SF-CIS)と fragment
spin difference(FSD)法を組み合わせて、
ペンタセ
ンとテトラセン、及びそれらの誘導体の分子結晶から二量体モデルを作成し、SF 相互作用の第一原理計算を実行した[5]。これらのアセン分子結晶は、構造だけでなく、SF 反応時間も実験的に分かっている[6]。そこで Hsu らの二量体モデルに我々の手法を適用し、その結果を比較した(表 1)。
表 1

算では、1(T1T1)を記述するためのスピン電子配置に大きな影響を与えるが、式(1)の分母で Etun と

S1S0 Ĥ CA CA Ĥ T1T1

D

セン(DTT)以外では、我々の結果と Hsu らの結果

SF 相互作用(meV)の計算結果
Dimer Model

Ours

Hsu’s

Pentacene

50.19

14.75

TIPS-P

23.52

43.83

DTP

6.12

8.21

DBTP

0.082

0.054

DBP

0.001

0.023

DPP

0.072

0.087

Tetracene

17.09

11.38

Ruburene

0.014

0.012

DTT

2.686

20.7

表 1 から、ペンタセンと 6,13-di(2’-thienyl)テトラ

10 − Academic Center for Computing and Media Studies, Kyoto University

ことから、スピン汚染やアンバランスな電子配置
が、Hsu らの結果を悪くしている可能性が高い。
今後、開発した手法を多参照摂動理論に適用できるように拡張することで、SF 相互作用への動的相関の影響をより精密に調べる予定である。

4 ACKNOWLDGEMENTS
本研究では、京都大学学術情報メディアセンタ
ーのスーパーコンピュータ共同研究制度(若手奨励枠)を活用させて頂きました。
[1]M.B. Smith and J. Michl, Chem. Rev. 110,
6891 (2010); Annu. Rev. Phys. Chem. 64, 361(2013).
[2]M.C. Hanna and A.J. Nozik, J. Appl. Phys.
100, 074510 (2006).
[3]H. Kitoh-Nishioka, R. Sato, and Y. Shigeta(in preparation).
[4] J. Ivanic, J. Chem. Phys. 119, 9364
(2003).
[5] C.-H. Yang and C.-P. Hsu, J. Phys. Chem.

Lett. 6, 1925 (2015).
[6]S.R. Yost, et al. Nature Chem. 6, 492(2014).

高解像度大気海洋波浪結合モデルを用いた海面物理モデルの現地観測結果に基づく検証
二宮 順一
金沢大学理工学域環境デザイン学系

1 はじめに

2 研究手法海洋の運動量や熱エネルギーは主に大気との間で交換されており，大気海洋間の相互作用を正しく理解することは海洋の流動や水温変化等の理解に直結する重要な課題である．海面での運動量輸
送はバルク式として，

2.1 数値モデル
数値モデルには大気海洋波浪結合モデルとして開 発 さ れ た COAWST (A Coupled Ocean
Atmosphere߬ ൌ ‫ܥ‬஽ ߩ‫ݑ‬ଶのように表現される．ここで，CD は海面抵抗係数，

ρは大気密度，u は 10m 高度での風速が用いられる．
海面抵抗係数は観測結果から一定値ではなく，風速や海面の状態との相関を示すことが知られているが，実用的な観点から海面抵抗係数は風速の関数として定義された形が主流となっている．

Wave

Sediment

Transport

Modeling System) を用いる．これは，全てオープンソースの気象モデル WRF，海洋モデル
ROMS，波浪モデル SWAN を MCT で結合したモデルである．

2.2 計算の概要
海面抵抗係数のモデルとして，本研究では風速

海面抵抗係数に関連する研究は古くから行われ依存の式（Charnock），波形勾配の式（Taylor &ている．
一定値から始まり，
風速の一次関数近似，

Yelland），波齢の式（Drennan）の 3 式を用いる．

海面の状態との関係性が議論されるようになってこれらの式は海面粗度を与えるものであるが，海きた．一方で，近年まで 25m/s 以下程度の風速域

面粗度は海面抵抗係数と 1 対 1 の関係になっていまでしかカバーしていなかったが，最近では高風る．

速下での観測結果が得られた．高風速下では，海2.3 現地観測結果

面抵抗係数は飽和して一定値もしくは風速の増加
とともに減少する傾向が示された．海面抵抗係数
は風速の一次関数として一定値まで増加し飽和する式や，一定値から減少する式，海面の状態の依存性のため，海面勾配や波齢の関数などの定式化がなされている．しかし，提案はされるものの検証することは全く進んでいない．そこで本研究では，大気海洋波浪結合モデルに様々な海面抵抗係数のモデルを適用して，高風速を伴う台風計算に検証に用いる現地観測結果は，アメリカと台湾の共同観測プロジェクト結果である ITOP
（Impact of Typhoons on the Ocean in the
Pacific）を用いる．このプロジェクトでは，台湾
東沖の台風を対象として，航空機やドロップゾン
デを利用した気象観測，海中固定式ブイや漂流式
フロート，ドリフターを用いた海洋観測が行われた．

どのような影響を与えるのか検討を行った．

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 11

図 2 瞬間値から台風中心付近 100km の範囲から抽出した風速と海面抵抗係数の関係
図 1 海面抵抗係数を変更した台風再現計算の結ルでは海面抵抗係数を低風速域でやや大きく推定

果 (a)台風経路（Best Track と各計算結果

し，30m/s を超えるような高風速域では海面状態

および観測機器の投下位置や経路を示に依存したモデルが大きく推定する傾向がある．

す．
），(b)中心気圧の時間変化，(c)最大風速
の時間変化

3 結果

3.2 海洋場に対する海面抵抗係数の影響
海洋場に対する海面抵抗係数の影響として，台風経路に最も近く，水温変化が大きいドリフター

気象場および海洋場に対する海面抵抗係数の影響について以下にまとめる．

3.1 気象場に対する海面抵抗係数の影響
海面抵抗係数を変更した台風の再現計算結果を
図 1 海面抵抗係数を変更した台風再現計算の結果 (a)台風経路（Best Track と各計算結果および観測機器の投下位置や経路を示す．
），(b)中心気圧の時間変化，(c)最大風速の時間変化図 1 に示す．
経路，中心気圧，最大風速の変化をそれぞれ示している．
経路，
最大風速に対する影響は小さいが，中心気圧は 1010hPa を基準として気圧深度
70hPa に対して，最大で約 10hPa の大きな差が見られた．台風の再現がまだ十分でなく，海面抵
抗係数モデルの正しさについて議論できないが，海面抵抗係数は台風の発達にも関連する重要なパラメータであることがわかる．
高風速域における海面抵抗係数と風速との関係
を図 2 に示す．風速の一次関数で示される直線的な分布と二次関数的な分布，広がった分布と三種の分布が確認できる．分布から明らかなように，
波浪依存のモデルでは同一の風速で様々な値を取りうる．また，台風下において，風速依存のモデ

12 − Academic Center for Computing and Media Studies, Kyoto University

による観測結果を用いて比較を行った．台風通過に伴って，通常は cold wake という表層付近の水温低下が生じるが，風速依存モデルによる結果では，この cold wake がやや弱く再現された．これは，30m/s を超えるような強風速下での強混合において，海面抵抗係数を小さめに推定して海洋への運動量輸送を比較的小さく見積もったためと考えられる．一方で，波浪依存モデルでは水温低下を過大に推定しており，これらの海面抵抗係数の分布幅の中に収まる分布があるものと考えられる．

4 まとめ
風速または海面の状態に依存する海面抵抗係数のモデルを適用した大気海洋波浪結合モデルを用いて，強風を伴う台風の計算を行った．
その結果，
低風速では風速依存のモデルが，高風速では海面
状態依存のモデルが海面抵抗係数を大きく推定する傾向が見られた．
また，
大気側への影響として，
本実験では海面状態依存のモデルが台風強度を強く推定した．海洋側への影響として，cold wakeの再現に注目すると，海面状態依存のモデルでは水温低下が過大に再現された．

iPS 細胞懸濁培養の最適化を目的とした撹拌槽内粒子挙動解析矢野雅貴
大阪大学大学院基礎工学研究科物質創成専攻化学工学領域

本研究では、iPS 細胞の三次元的な浮遊培養の実用化を目指し、撹拌槽内粒子挙動の数値解析を行った。細胞培養に悪影響を及ばすせん断応力と衝突頻度の二つの要素を求めて撹拌方法ごとに比較し、よりよい撹拌方法の知見を得た。その結果、振とう撹拌は底面に粒子がたまることで衝突頻度が増大することが分かった。この欠点を改善するため、一定時間ごとに撹拌方向を変更
する正逆交互回転のしんとう撹拌を開発した。実際にこの撹拌方法は通常よりも粒子を浮遊させることができると分かった。その理由は流体に上下方向の速度が発生したからである。

1 緒言( f vf )    ( f v f v f )t

iPS 細胞培養の産業化に向けて現在の熟練者による静置培養からプラントで行う自動化された懸濁培養にスケールアップすることが強く望まれている。
しかし、実用化のためには撹拌方法や速度、槽形状など多くのパラメータを最適化する必要があり、これを実験のみで行うのは困難である。そこで数値解析が最適化のための重要なツールとなりうる。一般に多能性幹細胞はせん断応力が作用する状況下で分化し  、細胞塊が底面中央部に集積することで不均一な大きさになる。そこで本研
究では、撹拌槽内の粒子挙動を解析する手法を開発し、一方向および正逆交互回転と往復運動の撹拌において細胞塊が受けるせん断応力や底面に接地する細胞塊数を算出し最適な撹拌方法についての検討を行った。


2 解析手法流動の支配方程式は連続式と以下に示す
Navier-Stokes 式である。U は容器自体が動く速度 p    (  f v f )   f g - ms


dU
 F
dt

回転運動
U   Rs sin t ,  Rs cos t , 0 
往復運動

 (1)

(2)

U   Rs sin t , 0, 0 
気液界面形状の解析には VOF 法を用いた。また
iPS 細胞塊を粒子とみなし、その挙動解析には以下の Newton の運動方程式を用いた。 


dv s
dU

(3)
 FD  Fg  ms
dt
dt
ここで添え字 f  s はそれぞれ流体と粒子を表しており、F は表面張力、F は流体抵抗、Fg はms

D

粒子に加わる浮力と重力、 g は重力加速度、m は粒子の質量、p は圧力、t は時間、v は速度、 は
粘度、 は密度、Rs は回転半径（12.5 mm）、 は

角速度（6.28 rad/s）である。解析領域は直径42 mm、高さ43.2 mm の円筒容器で撹拌速度は60 rpm とした。また、壁、撹拌翼と粒子とは完全弾性衝突すると仮定した。粒子数は192個、直径は0.2 mm、粒子と流体の比重は1.04とし、初期条件として粒子を等間隔に配置した。

ベクトルであり、撹拌方法によって異なる。なお
回転運動に関しては通常の一方向撹拌のみでなく、一秒ごとに回転方向が逆転する正逆交互撹拌に関してもを行った。

3 結果と考察
Figure 1 (a) に往復運動、一方向および正逆交互
回転運動時のせん断応力の平均値の時間変化を示す。回転運動の正逆交互撹拌、一方向撹拌、往復

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 13

運動の順に粒子に作用するせん断応力は小さくな   (a)った。また回転運動の正逆交互撹拌、往復運動ではせん断応力の値は大きく振動した。回転方向逆
転時に大きなせん断応力が作用することから回転方向を逆転する際には速度変化をつけて徐々に行う必要があると考える。Figure 1 (b) に底面に接地している粒子数の時間変化を示す。回転運動の正逆交互撹拌、一方向撹拌、往復運動の順に底面に接地する粒子数が減少した。Figure 2 に中心軸を通る底面から 10 mm までの断面における上下方
向の時間平均速度を示す。(a)より(b)の中心部の上  ( b)向き速度が大きいことから、正逆交互撹拌の方が底面に接地する粒子数を減少させることができたと考えられる。

4 結言
iPS 細胞の懸濁培養実用化へ向けて、懸濁培養時の細胞塊挙動を解析する手法を開発し、
往復運動、一方向および正逆交互回転運動撹拌におけるせん
断応力、底面に接地する細胞塊数を求めることができた。また底面に接地する細胞塊数を減らすた
めには中心部の上向きの流れが重要であり、今後そのような特徴がより強い撹拌方法を開発する必要があることがわかった。

 
Figure 1. Time variation of(a) average shear stress and(b) number of colonies on the bottom.

謝辞
 本研究は一部国立研究開発法人日本医療研究開発機構（AMED）の「再生医療の産業化に向け+た細胞製造・加工システムの開発」プロジェクト
および MEXT 科研費 15H04173の助成を受けた。
参考文献

[1] Adamo, L. et al., Nature, 459, 1131-1135(2009).
[2] Yamamoto, T. et al., J. Chem. Eng. Jpn., inpress.

14 − Academic Center for Computing and Media Studies, Kyoto University

Figure 2. Time averaged of horizontal speed at sectional viewat 30 seconds.

大規模津波による流出油の拡散挙動予測高木 洋平
横浜国立大学大学院システムの創生部門

1 はじめにSTOC による津波計算では太平洋を含む広域の計算メッシュから大阪湾石油コンビナート地域の将来の発生確率が高い南海トラフ地震では，東狭域の計算メッシュを作成し，計算格子を数段階

日本大震災と同様に，津波によるコンビナート地にネスティングすることによって解析を行ってい域からの油流出が懸念されている．大阪湾石油コる．大阪湾における大規模な油流出は図 1 に示す

ンビナート地域では，小規模なタンクが集積して北港地域及び堺泉北地域であり，これらの地域でいる北港地域と，大型のタンクが多数存在する堺

津波越流の可能性があるかどうかを津波計算の結

泉北地域が主な油流出による被害が予想される地

果より確認を行った．初期流出開始位置とする石域である．本研究では，この二地域から南海トラ

油タンク配置箇所での水位の時間変化を確認する

フ地震によって最大流出量の油が流出すると仮定と，各ポイントで海抜高度が異なるが，津波発生

し津波によってどれくらいの範囲に油が拡散するから 100 分〜200 分の間に越流が起こっているこ

かを数値シミュレーションによって検討した結果とがわかった．

を述べる．

北港地域及び堺泉北地域における想定される最大流出油は，大阪府石油コンビナート防災計画の2 解析手法

見積もりによれば，北港地域で 27227 kL，堺泉
北地域で 4838 kL がそれぞれ総量として推定さ本数値解析では広域の津波による油流出の挙動れている．北港地域でタンクを集約した 8 点，堺を解析するが，実際の解析ではまずはじめに南海

泉北地域で 4 点の流出開始ポイントから津波の越トラフを震源とする広域の津波計算を行い，得ら流によって流出・拡散する数値シミュレーション

れた津波計算データを用いてタンク群からの油流を実施した．計算には STOC と同じく港湾空港技

出挙動を計算する．

術研究所が開発した STOC-OIL(2)を利用した．

津波計算は港湾空港技術研究所が開発した

STOC-OIL では流出した油を仮想的な粒子に置

STOC（Storm Surge and Tsunami Simulator in

き換え，津波による移流や油粒子群の拡散を考慮

Oceans and Costal Areas）(1)を利用して行った．

しながら設定した粒子それぞれの運動軌跡を

STOC は静水圧近似を用いた広域の津波シミュレ

STOC による津波計算データを入力して追跡する．

ータであり，地形情報や津波発生源での初期水位

本解析では，前述の大阪府石油コンビナート防災

または断層の移動時間を入力することによって，計画の見積もりから，北港・堺泉北のそれぞれの計算領域内の各地点・各時刻における津波高さ・地域に5000粒子ずつ設定し，追跡計算を行った．

速度を出力として得ることができる．
本解析では，

仮想粒子一つあたりの体積は，北港地域で 5.42

地形情報については内閣府の南海トラフの巨大地

m3，堺泉北地域で 0.97 m3 となっている．油の流震モデル検討会が公開しているデータを利用し，出開始は，地震動による津波発生から各地域まで

津波発生の初期条件である断層モデルに関しては，

津波第一波が到達するまでの時間であり，北港地

同検討会が報告している「三重県南部沖～徳島県域で津波発生から 126 分後，堺泉北地域で 137 分

沖」と「足摺岬沖」に「大すべり域＋超大すべり

後とした．

域」を 2 箇所設定したケース 11 を利用した．

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 15

3 計算結果

出事故の再現計算, 土木学会論文集 B2, 70,
pp.15-30, 2014.

図 2 に北港地域，堺泉北地域での津波到達後 70
分経過後の時間における油拡散結果を示す．北港
地域では小型の石油タンクが多く初期流出量も多
いため，陸側・海側の両方に多量の油が徐々に拡散してゆく様子が計算された．陸側においては，追跡計算をさらに継続して流出開始から 90 分が経過すると，住宅やオフィスが密集する難波及び梅田まで油が拡散しており，東日本大震災時の気仙沼地域のように津波火災が発生すると，甚大な被害を被ることが推定される．海側に流出・拡散する油粒子も多く，長期的な油流出が続くと海流・潮流の影響によって大阪湾全体に油が拡散し，大阪湾に停留する船舶や神戸地域も被害を受けることが予想される．一方，堺泉北地域での油流出範囲に注目すると，北港と同じく陸側の住宅街と海側の両方に油粒子が拡散していることがわかる．
陸側への流出は南北に走る上町断層で止まってお図 1 大阪湾における石油コンビナート地域

り，津波発生直後の避難としてはこの断層よりも内陸に避難することが重要である．

4 まとめ
本研究では広域津波シミュレーションと油粒子追跡計算を組み合わせた解析を，南海トラフ地震
によって想定される大阪湾石油コンビナート地域
で津波被害に対して行った．北港地域，堺泉北地域のタンク密集地域からの最大流出量見積もりを初期条件として流出油の津波による拡散範囲を数値予測した結果，両地域において市街地への多量の油拡散が発生し、津波火災の危険性があることがわかった．今後はシミュレーションの精度向上
や火災シミュレーションへの連結を検討し，災害予測・防災技術の発展に貢献したい．

5 引用文献(1) 冨田孝史・柿沼太郎：海水流動の 3 次元性を考慮した高潮・津波数値シミュ レータ
STOC の開発と津波解析への応用, 港湾空港技研研究所報告, 44, pp.83-98, 2005.
(2) 松崎義孝・藤田勇：海水面における流出油の拡散・移流に関する数値計算法の開発と油流

16 − Academic Center for Computing and Media Studies, Kyoto University

図 2 大阪湾における越流発生後 70 分での油拡
散範囲: (上) 北港，(下) 堺泉北

視聴覚モダリティ間デコーディングによる感覚間協応のメカニズムの解明
金谷翔子
京都大学 大学院人間・環境学研究科 共生人間学専攻 認知・行動科学講座

1 背景と問題た図形（鋭角図形）二つと、曲線で構成された図形（曲線図形）二つを用意した。また、聴覚刺激

人間の脳においては、視覚、聴覚、触覚といっとして、ソフトウェアシンセサイザーによって音た異なる感覚モダリティの情報は、異なる場所で、響特徴を複雑に操作した電子音を多数作成し、そ異なる方法によって表現される。しかし、我々はの中から上述の四つの視覚刺激のそれぞれに最も時に、別々の感覚モダリティの情報の間に何らか合うと感じられるもの四つを、多数の学生に対しの共通性を感じることがある。例えば、楽器の絵て行ったアンケートの結果に基づいて選択した。

や写真（視覚情報）と音（聴覚情報）はいずれも、

図１は、これらの計 8 種類の視覚刺激、および計

例えば「バイオリン」といった共通の表象に結び

8 種類の聴覚刺激のスペクトログラムを表す。

つく。また、このような経験的な結びつきがない
場合でも、例えば高い音には明るい色、低い音には暗い色が合うというように、本来は無関係ない
複数の特徴の間に非恣意的な対応が感じられることがあり、この現象は感覚間協応と呼ばれている[1]。複数感覚モダリティ間の共通性が、脳のどのような働きによって表現されているのかを調べる
ことは、感覚モダリティに依存しない高次知覚、認知を支えるメカニズムの理解につながると考えられる[2]。本研究では、fMRI を用いて実験参加者が様々な視聴覚刺激を観察している際の脳活動を計測し、このデータに MVPA（Multivariate

図 1 実験に用いた刺激

Pattern Analysis）解析によるモダリティ間デコ(a)風鈴のイラストと音、(b)絃楽器のイラストと音、(c)ーディングを適用することで、視聴覚情報の共通

鋭角図形と適合する音、(d)曲線図形と適合する音。聴性に関わる神経基盤の検討を行った。

覚刺激のスペクトログラムは横軸が時間、縦軸が周波数、色がパワーを表す。

2 実験方法これらの視覚刺激、および聴覚刺激を一つずつ、男女合わせて 16 名の成人が実験に参加した。

モニタおよびイヤフォンによって呈示し、これを経験に基づく視聴覚情報の対応について検討する

実験参加者が視聴している際の脳活動を fMRI にため、風鈴および弦楽器のイラスト（二つずつ）、

よって計測した。全ての刺激を一回ずつ呈示する

風鈴および弦楽器の音（二つずつ）を用意した。

系列をブロックとし、実験全体で 10 ブロックを

感覚間協応については、
特によく報告されている、繰り返した。fMRI から得られる生の脳活動デー

視覚刺激の形状と聴覚刺激の音響特徴の対応に着タに様々な下処理を加えた上で、以下の手順で目した。まず、視覚刺激として、鋭角で構成され

MVPA 解析を行った。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 17

3 解析方法
ある特定のブロック（e.g. ブロック 1）を除いた全てのブロックで取得した脳画像データを用いて、ある二カテゴリの刺激（e.g. 鋭角図形と曲線図形）を観察している際の脳活動を判別するように、
コンピュータプログラムに機械学習させる（デ
コーダの作成）。その上で、学習時に除外したブロック（e.g. ブロック 1）で取得した脳画像データ
をデコーダに入力し、学習させた二カテゴリの弁別が正しく行えるかどうか（e.g. 実際に鋭角/曲線図形を観察していた際の脳活動データが、正しく
鋭角/曲線図形と判定されるかどうか）をテストする。さらに、テスト時に用いるブロックを入れ替えて、これを 10 回繰り返し（e.g. ブロック 1、ブロック 2、…ブロック 10）、全ブロックにおける平均正答率を求める。
このような学習とテストを、
視覚情報と聴覚情報の両方を表現していると考えられる様々な小領域(Region of Interest, ROI)に

図 2 MVPA 平均正答率

対して行い、特に正答率の高くなる脳領域を探し

上段は風鈴と絃楽器の弁別、下段は鋭角図形/音と曲線

た。ROI として、全視覚刺激に対する平均脳活動

図形/音の弁別における正答率。LH は左半球、RH は右と、全聴覚刺激に対する平均脳活動がどちらも一半球の結果を表す。

定以上の高い値を示した、楔部(Cuneus)、楔前部(Precuneus)、内側前頭回(Medial Frontal Gyrus,

中心後回が多くの条件下で高い正答率を示した

MFG)、中心後回(Post Central Gyrus, PCG)を設ことから、この領域において視聴覚情報の共通性定した(図 2)。特に正答率の高くなる領域においてが比較的頑健に表現されている可能性がある。まは、当該の二カテゴリの刺激を識別するための情

た、弁別させる刺激の種類（風鈴/絃楽器、鋭角/報が表現されていると考えられる。

曲線）によって異なる半球の中心後回が高い正答率を示し、かつ楔部や内側前頭回の正答率も刺激

4 結果

によって異なっていたことから、経験的に得られる視聴覚情報間の結びつきと、感覚間協応のよう本研究では、
抽象刺激については鋭角か曲線か、
具体刺激については風鈴か絃楽器かを弁別させるな抽象的な視聴覚情報間の結びつきでは、異なるメカニズムが関与している可能性がある。

ように学習を行った。また、学習時には視覚(V)
刺激観察時の脳活動、テスト時には聴覚(A)刺激観

参考文献

察時の脳活動を用いるV-A モダリティ間デコーデ

1. Parise, C.V. and Spence, C., Audiovisualィングと、学習時には聴覚(A)刺激観察時の脳活動、

crossmodal correspondences and sound

テスト時には視覚(V)刺激観察時の脳活動を用い

symbolism: a study using the implicit

る A-V デコーディングの両方を試みた。上述のassociation

ROI における MVPA の平均正答率を図 2 に示す。

Research, 2012. 220(3-4): p. 319-333.

test.

Experimental

Brain

チャンスレベル(50%)よりも有意に高い正答率を

2. Man, K., et al., Sight and sound converge to示したのは、右中心後回(V-A、風鈴/絃楽器)、右

form modality-invariant representations in

楔部と右中心後回(A-V、風鈴/絃楽器)、右内側前temporoparietal

頭回と左中心後回(V-A、鋭角/曲線)であった。

Neuroscience, 2012. 32(47): p. 16629-16636.

18 − Academic Center for Computing and Media Studies, Kyoto University

cortex.

Journalof‫࠵ࡵٴߘۿ‬ՃՕஞ⎬⎋⎒ ⏌ ࡵ࠵ՃՕஞ⎯ѽॾॵ෱ƄխງƄ (λֹ੒Ӂ੒Ӂώ࢙Ӂ٥֋ѕ)
yoshi.saiki@r.hit-u.ac.jpҞ์1⎬⎋⎒⏌ Siegel ball ⏓ ૛⎺⎠ゎ ⎔ ⎔ ⎨␍␥ ␝ ⏻ α ⎩ β ⏓ ろ

α = 2 cos θ cos φ − cos2 φ, β = cos θ + i sin θ ⎩ ⎃⏊ ⏐⎚⎔⎩ ⎬ ⎚⏌ ゎ Hα,β ⎯೑ஓ୅⎰ろ y∗ = cos θ ± (cos θ − cos φ)ඕۣ⎨⎰ろ ‫࠵ࡵٴߘۿ‬ՃՕஞ⎬ ⎋ ⎒ ⏌ Ѿୃॅ⎯‫ۿ‬ਙ‫ظ‬‫ޔ‬൚⎩ ⎞ ⎯஦ड़⏓ ෗⏊ ⎌⎬ ⎘ ⎠ ゎ ⎞ ⎯Ҟ์⏓ ೾৔⏤␊ ␯⎨⎃ ⏌ ⎍ろ y∗ = cos φ ⎯൘⎬ મෳ⎘ ⎧ µ = cos φ + i sin φࠃਉ⎯ Siegel ball ⎬ ԥ⎿ ⏍⏌ ࡵ࠵ՃՕஞ⏓ ຶ⎬ ⎩ ⎤ ⎧ ࣅ ⎩ ⎚⏌ ⎩ ろ ⎔ ⎯೑ஓ୅⎬ ⎋⎒ ⏌ ␠ ⏯ ␏۱ໆ⎯‫چ‬ยખ⎰ βµ
ѻ⎚⏌ ゎ⎩ βµ ⎨⎃ ⏋ ろ ⎞ ⎯ഫҺ⎰ θ + φ and θ − φ ( mod 2π)⎨ ⎃ ⏌ ゎ ␍␥ ␝ ⏻ ⏓ (θ = 0.664, φ = 2.032) ⎩ ঎଑⎘ ろ࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ה‬2

ࡾՃખ (x0 , y0 ) = (−0.500 + 0.126i, −0.387 − 0.163i)⎩ ⎘ ⎠ ⎩ ⎎ ⎬ ண⏊ ⏍⏌ Օஞ⏓ ۪‫ ⏌ ⎚ݺ‬ゎ Օஞ⏓ ࠄϳ␌⍬⏯ ␑ ⎯⏤␧⏰ ⍬␅ ଑๲⎰ろ ⏤␧⏰ ⍬␅ ଯ⎫ ຢӁ‫ذ‬⎘ ⎠ ⏃ ⎯ ⏓ य 1 ⎬ ߟ⎚ゎ ܳय⎰ x ⎯ ߹ೱ⎩ ֚ೱ⏓ ࠄ⎯Օஞ (xn ) ⎬ ਻⎘ ⎧ ଑զ⎖ ⏍⏌ ъ⏊ ⎌⎯ຟ f ⎬ ਻⎘ ⎧ ろϳ⎘ ⎠ ⏃ ⎯ ⎨ ろ ϓय⎰ઢ⏍�ߗ ԝ ܽ ವ (rn , rn−2 ) ⏓ ࠄ−1⎞ ⎯␌⍬⏯ ␑ ച‫ ה‬B N (f )(x) = ΣN
n=0 f (xn )/N ⎰ろ Օஞ�ૠ N → ∞ ⎨‫׸‬ԝച‫ ה‬f dµ ⎬࠴ਖ⎚⏌ ⎩ ⎅⎇ ⏃ ⎯⎨⎃=(Rexn + 0.4)2 + Imxn 2 ⎨⎃ ⏌ ゎ ⎞ ⏍ ⎟ ⏍ ろ (Rex, Imx) = (−0.4, 0) ⎫ ⏊ ⎴ ⎬ϳ ⎘ ⎠ ゎ ⎔ ⎔ ⎨ rn

⏌ . ੫⎘ ろ µ ⎰ຢӁ‫⎯ذ‬೑ബҷ๻ਗ୞⎨⎃⏌ . ѱゖ ⎰ろ ࡵ

(rn , rn−2 ) = (0.14, 0.145) ⎌ ⏊ ‫ ظ‬ਗ ⎘ ⎠ Һ ୞ ബ ш ⏓࠵ՃՕஞ⎬ ୀ⏐⏌ үࠦ⎯ຟ⎯‫ۿ ⏓ޔظ‬ਙш⎖ ⎜⏌ ⎔ ⎩ ⎍य 2(ܳ ࣘ)ろ (ϓ ࣘ) ⎬ ߟ ⎘ ⎠ ゎ ⎔ ⏍ ⏊ ⎯ Һ ୞ ബ ш ⎯⎨⎎ ⏌ ࠡ൚⏓ ଚΕ⎘ ⎠ [1, 2].

␌ ⍬ ⏯ ␑ ച ‫ ⏀ ࡚ ⏓ ה‬೒ ⎎ ␌ ⍬ ⏯ ␑ ച ‫( ה‬WB ⎫ ⏊࡚⏀ԟॅ w : R → [0, ∞) ⏓ ෌‫ځ‬Ѿಘഇя௾⎨⎌⎥ (0, 1)⎴ ⎬ WB2 ) ⎬ ⏉ ⎤ ⎧ ‫ ⎩ ⏌ ⎚ޔظ‬ろ N = 106 ⎯ ⏦ ⍬⎨ w > 0ろ ⎞ ⏍Ιҙ⎨ w = 0 ⎩ ⎫ ⏌ ⏃ ⎯⎩ ⎚⏌ . Ιч⎨⏼ ⍬ ⎯ Օஞૠ⎯ ࣟ‫ ⎬܄‬4 ఝ३୞⎯ Ԥ⏂ ‫ܴڬ‬ଟ୞⎯ ‫ۿ‬⎰ろ ஦⎬

wp (t) :=��exp p −1 p ,t (1−t)0,३୞⎨ Ѿୃॅ⏓ ஔ଑⎚ ⏌ ⎔ ⎩ ⎍ ࡥ๝⎠ ゎ ⎳ ⎩ ⎥ ⏂ ⎯for t ∈ (0, 1)

Ѿୃॅ⎰

(1)

⎨ ⎃ ⏋ろ ⎔ ⏍⎰ρ1=
(θ + φ)/2π
0.42908172657574982523291062605229878ろ ⏃ ⎇for t ∈/ (0, 1).

⎩ ⎚⏌ . ⎔ ⏍⏓ ๆ⎅⎧ f ⎯࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ ה‬WBN
�N −1
⏓ WBN
p (f )(x) :=n=0 ŵn,N f (xn ) ⎩ ଑զ⎚⏌ . ੫⎘ ろ

൘ ⎰ ろ⎨ ⎃ ⏋ ろ ρ1w (n/N )p⎨⎃⏌ . ຢӁ‫ ⎩ذ‬f ⎍ C ∞ ⎨ろ wp ⏓ŵn,N = ΣN −1
w(j/N )

ࣘ⎨ߟ⎘ ⎠ड़߸⏓ ⏃ ⎥ C ∞ ࡚⏀ԟॅ⎩ ⎘ ⎧ ろ ۖ⎬ࠚӺ⎯щ଑⏓ ю⎉ ⏌ ⎩ ろ ௥Τ⎯ߢীॅ m ⎬ ਻⎘ ⎧ ろ ⎃ ⏌ Cm > 0���⎍ਣ‫ ⎧ ⎘ݚ‬ろ �WBN (f )(x) − f dµ� ≤ Cm N −m ⎨⎃⏌ ⎔⎅⏌ ⏃ ⎯⎯ろ ⎞ ⎯ভ‫⎨آ⎯܄ٌآ‬ԟ໓⎦⎅⎧ ⎅⏌ ⎔ ⎩ ⎬⎩ ⎍ߟ⎖ ⏍⏌ [2, 3].

શΤ⎘ ⎠ ⎅ [5]ゎ ⎫ ⎋ろ ⎔ ⎯३୞⎯Ѿୃॅ⏓ ࡓ๝⎯␌⍬⏯
␑ ച‫ ה‬B N ⎨ண⏌ ⎠⏂⎬⎰ろ N 30 ଟ୞⎯Օஞૠ⏓ ์⎚⏌ ゎ೾৔⏤␊ ␯ ࠃਉ

4

Ushiki (2016) [4] ⏓ ‫⎬۪ދ‬೾৔⏤␊ ␯ ࠃਉ Hα,β : C2 →C

2((α, β) ∈ C )
� � �
�
x
y
Hα,β
=
y
β(y 2 + α) − β 2 x,

0.64680568872546264456474361434588417
− ρ2
=
2φ/2π
=

⎬ ண⏊ ⏍⎠ ௖⎥⎯Ѿୃॅ⎰೑ஓ୅⎯‫چ‬ยખ⎬ ԟ໓⎦⎅⎧

p

2

=
λ

0.64680568872546264456474361434589842 ⎬ ⎞ ⏍ ⎟⏍ 31 ‫⎯ن‬३୞⎨λઠ⎘ ⎧ ⎅⏌ ゎ ⎔ ⎯ٌї⎌⏊ ⏐⎌⏌ ⏉ ⎇j=03

0.42908172657574982523291062605231270⎿⎩⏂
ඕۣ⎨⎰ , ⏤␧⏰⍬␅ ଯ⎫ ຢӁ‫ ⏌ ⎒⎋⎬ذ‬үࠦച‫ה‬ຟ⏓‫ ⎬⏂ ⎠ ⏌⎚ޔظ‬ๆ⎅⏊ ⏍⏌ ␌⍬⏯ ␑ ച‫࠵ࡵ ⏓ޔظה‬ՃՕ(2)ஞ⎯ࣟ‫ۿ ⎬܄‬ਙш⎖ ⎜⏌ ⎔ ⎩ ⎍⎨⎎ ⏌ ࡚⏀೒⎎ ␌⍬⏯ ␑

1
京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 19ച‫ ⏓ה‬ๆ⎅⎧ ろ ೾৔⏤␊ ␯ ࠃਉ⎯バ ߘ‫࠵ࡵٴ‬ՃՕஞ⏓ バ
ࠦຳ⎯バ ߘ‫ٴ‬෣⎬ ࠄϳ⎘ ⎠ ⎩ ⎎ ⎯Ѿୃॅ⏓ ‫ ⎠ ⎘ޔظ‬ゎ ⎞
0.2

0.25

⎘ ⎧ ろ ண⏊ ⏍⎠ バ ⎥⎯Ѿୃॅ⎩ ࠃਉ⎯೑ஓ୅⎯‫چ‬ยખ⎩

0.2⎯ԝ⎯ԟ‫ ⎠ ⎘ݺ۪ ⏓ؙ‬ゎ

0.1

0.15

0

0.1

-0.1

5

0.05

-0.2
-0.65 -0.55 -0.45 -0.35 -0.25

0

0

ࠌߤ
ඕۣ⎨ๆ⎅⎠೾৔⏤␊ ␯ ࠃਉ⎬ԟ⎘ ⎧ ろ ֣୙੒Ӂ⎯ϔ೜

0.05 0.1 0.15 0.2 0.25

࡚ᆘগ१⎬ ࣸゖ ⎩ ⎕ָߟ⎅⎠ ⎡⎅⎠ ゎ ӽࠌऎ⎘ ࣘ⎓⎿ ⎚ゎय 1:(a) (Rex, Imx) ച෣⎹ ⎯ Օஞ⎯ ࠄϳ.

⎙ Օ ஞ ⎯ (rn , rn−2 ) ച ෣ ⎹ ⎯ ࠄ ϳ ゎ ⎠ ⎡ ⎘ ろ(Rexn + 0.4)2 + Imxn 2 .

(b) ஔ
rn

=

‫۪ދ‬ഓ٤
[1] S. Das, C. B. Dock, Y. Saiki, M. Salgado-Flores,
E. Sander, J. Wu and J. A. Yorke. Measuring
quasiperiodicity, Europhysics Letters 114, 40005.

1
lift of angle difference

lift of angle difference

1
0.8
0.6
0.4
0.2
0

0

10

-15

0.4

0

2000 4000 6000 8000 10000
N

10[5] S. Das, Y. Saiki, E. Sander and J. A. Yorke. Solving
the Babylonian Problem of quasiperiodic rotation

10-15
10-20
10-25

10-30

10-30
104

N

105

106

10-35
104

(a) ࠃ ਉ ඡ ⎯ (Rex, Imx)

rates, preprint, arXiv:1706.02595.

105

=

106

N

107

(−0.4, 0) ⎌⏊ ⏀ ⎠ Һ ୞ ബ ш (य 1(a) ⎬ ਻ Ы).

(b) ࠃ ਉ ඡ ⎯(rn , rn−2 ) = (0.14, 0.145) ⎌ ⏊ ⏀ ⎠ Һ୞ബш (य 1(b)⎬ ਻Ы).

[4] S. Ushiki. Abundance of Siegel balls in a familyof Hénon maps, 2016, http://www.math.h.kyotou.ac.jp/˜ushiki/papers/AbundanceArt.pdf

WB
WB2

-10

10-25

य 2:

arXiv:1506.06810.

0.2

10-5

WB
WB2

10-20

10-35
103

[3] S. Das and J. A. Yorke. Super convergence of ergodic averages for quasiperiodic orbits, preprint,

0.6

0

error

error

10

publication in Nonlinearity, arXiv:1601.06051.

0.82000 4000 6000 8000 10000
N

10-5
-10

[2] S. Das, Y. Saiki, E. Sander and J. A. Yorke. Quantitative quasiperiodicity, Nonlinearity , accepted for(c) (a) ⎯ Һ୞ബш⎯ ࡚⏀ ೒⎎ ␌ ⍬ ⏯ ␑ ച= (θ +‫( ה‬WB and WB2 ) ⎯ ρ1 ⎹ ⎯ ࠴ ਖ (ρ1φ)/2π = 0.42908172657574982523291062605229878)ゎ
(d) (b) ⎯Һ୞ബш⎯࡚⏀ ೒⎎ ␌ ⍬ ⏯ ␑ ച‫( ה‬WB and
WB2 ) ⎯ ρ1 − ρ2 ⎹ ⎯ ࠴ਖ (ρ1 − ρ2 = 2φ/2π =0.64680568872546264456474361434589842)ゎ

20 − Academic Center for Computing and Media Studies, Kyoto University

メモリ一貫性モデルを考慮したモデル検査器の開発
安部 達也
千葉工業大学人工知能・ソフトウェア技術研究センター
abet@stair.center

1

はじめに

をこれまで行うことができなかった。

メモリ一貫性モデルとは、複数のスレッドが共有メモリを利用する際の取り決めである。この取り決めが厳しいものである場合、例えば、あるスレッド
が共有メモリを利用している際、他のスレッドがその利用を待たされるといったものである場合は並列処理の効率が下がってしまう。このため、現在、多くのコンピュータアーキテクチャやプログラミング言
語がゆるいメモリ一貫性モデルを採用している。しかし、メモリ一貫性モデルがゆるすぎると、複数のスレッドによる共有メモリの誤操作により共有メモ
リ上のデータの整合性が崩れてしまう可能性があり、結果として正しい計算が行われないといったことがメモリ一貫性モデルを入力にとれるモデル検査器の開発にはいわゆる状態爆発の悪化という問題がある。そもそもモデル検査は、入力されるモデルが複数のスレッドが立ち上げられるものであれば、各スレッドが実行する命令の実行タイミングによって生成される状態の数の組合せ的増大問題、いわゆる状態爆発問題をはらんでいる。さらに悪いことに、ゆるいメモリ一貫性モデルは従来は許容されなかった挙動を許容することがあるため、とり得る状態の数はさらに増大しがちである。今回、この状態爆発の悪化問題を緩和するための最適化とヒューリスティッ
クスを提案し、その実装と評価を行った。

起こり得る。そこまでゆるすぎないにしても、プログラマの直観に合わないゆるいメモリ一貫性モデルを採用したことにより、意図しない挙動を行うプロ
グラムをプログラマが書いてしまう恐れがある。こ2

ステージ最適化とリオーダリン
グ制御ヒューリスティックスのため、どのメモリ一貫性モデルを採用するかは処

ゆるいメモリ一貫性モデルの下では命令の持つ複理の効率とプログラムのわかりやすさとのトレード数の効果が互いにインタリーブして現れることがあ

オフの問題を抱えており、あるメモリ一貫性モデル

る。採用されるメモリ一貫性モデルを解析することがいつも優れているというものではない。ゆえに、で、可能な時は命令の持つ複数の効果を個別に扱うあるメモリ一貫性モデル一つだけを固定して考えてのでなく、それらを統合（ステージと呼ぶ）して扱おけばよいというものではない。

うことにより状態数の削減を行う、ステージ最適化

今回、著者らが開発しているメモリ一貫性モデルを提案しその実装と評価のための実験を行った [3]。

を入力にとれるモデル検査器 McSPIN のさらなるまた、モデル検査の反例発見法の一つであるヒュー改良を行った [1]。モデル検査とは、アルゴリズムやリスティク探索において、メモリ一貫性モデルを考プロトコルなどをモデリングして得られるプログラ

慮したモデル検査に特有のリオーダリング、特に、リム（モデルと呼ばれる）を、その挙動の網羅的探索

オーダリングされる命令の数の増減に着目したヒューを行うことで検査するというプログラム検証手法のリスティックスである、リオーダリング制御を提案一つであり、モデル検査器とはその実装をいう。既し、その実装と評価のための実験を行った [2]。

存の多くのモデル検査器はゆるいメモリ一貫性モデ

実験に使ったプログラムは並行コピーガーベジコレ

ルを考慮していないため、ゆるいメモリ一貫性モデ

クションアルゴリズム Chicken、Staccato、Stopless

ルを採用している現代の計算機上のプログラム検証

である。メモリ一貫性モデルは逐次一貫性（S）、トー

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 21

表 1: ステージ最適化の効果
MCM

algorithm
chickenTSO

staccato
stoplesschickenPSO

staccato
stopless

loop
1
2
1
2
1
2
1
2
1
2
1
2

✓
✓
✓
✓
✓
✓
✓
×
✓
✓
✓
×

without stages
statememory
(K)
(MB)
108
8,595
2,506
546,038141
14,9183,888
1,097,491
90
28,183564
564,635
308
25,208
1,136
264,857
143
15,1664,020
1,144,602
177
55,21045
45,416time(sec.)
132
8,637
236
16,022
378
7,705
430
4,248
243
16,975
833
630

✓
✓
✓
✓
✓
✓
✓
×
✓
✓
✓
×

with stages
statememory
(K)
(MB)
23
908
534
24,96026
1,032
733
43,432
14
1,404
87
18,88565
1,652
237
12,19026
1,032
751
44,920
30
2,520
8
2,148

time
(sec.)
8
433
11
735
19
585
28
227
11
768
46
41

mem.
ratio
(%)
10.6
4.6
6.9
4.0
5.0
3.3
6.6
4.6
6.8
3.9
4.6
4.7

net mem.
ratio
(%)
8.8
4.5
5.8
3.9
4.4
3.3
5.9
4.5
5.7
3.9
4.3
4.4

time
ratio
(%)
6.2
5.0
4.7
4.6
5.0
7.6
6.6
5.3
4.7
4.5
5.5
6.6

表 2: リオーダリング制御の効果
algorithm
chickenstaccato

stopless

M
C
M
S
T
P
R
S
T
P
R
S
T
P
R

the original
mem.
time
(MB)
(sec.)✓ 26522369.69✓ 48163709.24× 29328455.63———✓ 51050697.99✓ 825961198.35✓ 84669
1243.83———✓ 24361438.72✓ 37472658.28× 254342.64
× 15330
262.49the 1-local
mem.
time
(MB)
(sec.)✓ 26522768.32✓ 481631333.47
× 19240
300.38———✓ 51050
1443.83✓ 825962365.73✓ 84669
2432.60———✓ 24487
937.88✓ 374721349.94
× 193433.18
× 11309193.97the singular global
mem.
time
(MB)
(sec.)✓ 26522929.13✓ 481631372.40× 8586
175.73× 9792194.73✓ 51050
1756.46✓ 825962391.03✓ 84669
2463.15———✓ 24361
1407.18✓ 37472
1755.94
× 165338.54
× 166935.61

the increasing
mem.
time
(MB)
(sec.)✓ 26522367.57✓ 48163
707.45× 29328456.79———✓ 51050697.98✓ 825961197.39✓ 84669
1238.68———✓ 24361457.30✓ 37472715.72× 254339.40
× 15330264.55the decreasing
mem.
time
(MB)
(sec.)✓ 26522379.39✓ 48163712.32× 29662460.57———✓ 51050696.74✓ 825901202.74✓ 846631245.15———✓ 24361
386.11✓ 37472636.66× 240939.94
× 15129249.82

タルストアオーダリング（TSO または T）、パーシャ

終了しなかったケースはリラックストメモリオーダ

ルストアオーダリング（PSO または P）、リラック

リング下における Staccato の挙動の検査だけであっストメモリオーダリング（R）である。

た。このケースのみを京都大学スーパーコンピュー

表 1 はステージ最適化の効果を示すための実験のタシステム C により再実験を行い、このケースには結果であり、[3] の記載を引用・省略したものである。

反例が存在しない、つまり、リオーダリング制御のいずれのケースにおいてもメモリ使用量・実行時間

適用範囲外であることを確認した。

ともに改善している。通常、このような実験を行った場合の最適化による改善は、メモリ不足により実験が完了しないケースに最適化を適用すると実験が完了するというように質的に示されることが多いが、今回、1.5 TB という非常に大きな共有メモリを有する京都大学スーパーコンピュータシステム C を利用することによりどのくらい改善したかを量的に示すことに成功した。
表 2 はリオーダリング制御の効果を示すための実験の結果であり、[2] の記載を引用・省略したものである。表 1 とメモリ使用量や実行時間に違いが見られるのは、これら二つの実験の間にモデルがより詳細化されたことと McSPIN の性能が改善されたことによる。この実験では本稿では省略しているが多くのケースを扱ったため共有メモリを 128 GB に制限し、それを超えた場合は検査を打ち切った。メモリ不足によりリオーダリング制御を適用しても検査が22 − Academic Center for Computing and Media Studies, Kyoto University

参考文献[1] T. Abe and T. Maeda. A general model checking framework for various memory consistencymodels. In Proc. of HIPS, pages 332–341, 2014.
https://bitbucket.org/abet/mcspin/.
[2] T. Abe, T. Ugawa, and T. Maeda. Reordering
control approaches to state explosion in model
checking with memory consistency models. In
Proc. of VSTTE, 2017. To appear.
[3] T. Abe, T. Ugawa, T. Maeda, and K. Matsumoto. Reducing state explosion for softwaremodel checking with relaxed memory consistency models. In Proc. of SETTA, volume 9984
of LNCS, pages 118–135, 2016.

固液混相流の数値粘性測定
野村怜佳 ∗ ，高瀬慎介 † ，森口 周二 ‡ ，寺田賢二郎 ‡
∗

†

東北大学

八戸工業大学
‡

東北大学

はじめに

1

工学研究科

工学部

土木建築工学科

災害科学国際研究所

レオロジカルな挙動の捕捉が可能となる．本研究で未曾有の被害規模となった 2011 年東日本大震災に際しては，従来の津波予測シミュレーションで精確な評価がされなかった土砂移動による被害拡大，復旧・復興の遅延が報告され，土砂移動と津波の同時解析の必要性が強く認識された．このような状況を打開するために，土粒子レベルの微視的な運動機構が反映された新たな固液混相流モデルによる広域流体土砂連成解析の必要性に迫られている．本研究は流体運動を特徴づける代表的な物性値である粘性には，細管粘度計の測定原理である Hagen-Poiseuilleの法則を用い，管路内流れの数値シミュレーションを行うことによって等価な均質流体としての固液混
合流体の粘性を測定する方法を提案する．土砂粒子の混入した流体の管路内の流れの現象の数値シミュ
レーションを実施して，その後処理において管路内の流れのミクロ変数を空間および時間的に平均化し，その値を以下に示す Hagen-Poiseuille 式に代入する
ことで，粘性を得る．

着目し，固液混相流の等価な均質流体としての，見かけの粘性を評価するための数値粘性測定手法を提案しさらに，複数の数値粘性測定試験を行った．流体運動の記述には三次元 Navier-Stokes 方程式を用い，土砂を剛体として個別要素法（DEM）で表現す
る．また，水と土砂の相互連成解析には安定化有限被覆法 [1] を用いる．

2.1

πr4 Px8Q(1)

管路を用いた数値シミュレーションでは，入力条件として式 (1) の圧力勾配 Px を既知量として与えて付随するひずみ速度を算出することができる．したがって本研究では，圧力勾配 Px を任意に変化させて複数の解析を行い，提案する数値粘性測定法を適用して管路内での平均流量 Q を算出し，式 (1) から直接粘度を算定する．

粘性測定手法2µ=Hagen Poisuille 式による数値粘性
測定

一般的に粘性が一意に定まらない非ニュートン流体の場合，ひずみ速度−粘性関係を類推する必要があり，ひずみ速度を任意の値に変化させた粘性測定を複数行う必要があるが，そのような測定が行えるのは，種々の計器の中で回転粘度計のみである．数値解析による粘性測定の場合，どのような測定系であっても，解析条件を自由に設定でき，どのような測定原理を用いても，ひずみ速度を自由に設定することで粘性測定のための数値解析手法

3
3.1

安定化有限被覆法による流体固体連成
解析

本研究では流体土砂の相互連成解析として有限被覆法（Finite Cover Method，以降 FCM) を採用す
る．FCM は有限要素法 (FEM) と異なり，要素内に流体と構造の境界が存在することを許容するため，
その境界位置を正確に考慮した流れ場を求めることが可能である．本研究では，一次の四面体有限要素

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 23

の形状関数を重み関数として用いることにする．三2.0次元流れ場の支配方程式に対して SUPG/PSPG 法

近似曲線

3.2

µ /µ l[2] を適用し，安定化有限被覆法による解析を行った．

土砂の剛体球モデリング

N=243

1.5

N=116

本研究では，流体内における土砂粒子を流体力が

N=67

1.0

運動の駆動力として作用するものの，変形するほどの大きさを持たないものと仮定し，剛体として扱う

00.05
剛体粒子の体積分率ことにする．解析手法として，粒子間・粒子壁間接触の表現を簡易に行える個別要素法（Distinct Element

Method 以降，DEM）を採用する．土砂の形状に関しては，既往手法において用いられてきた土砂の代

0.1

図 1: 数値粘度測定の結果：体積分率−粘性関係
velocity[cm/s]
0.6

表粒径などを基に，全剛体粒子を完全な球形として0.4

モデル化する．

4

N=310

0.2
0.0

数値粘度測定結果-0.015(a) 管路内流れの様子
Velocity[cm/s]
0.0

0.028

-0.017

Velocity[cm/s]
0.0
0.015

土粒子の体積率と粘性の影響について明らかにするために，圧力勾配 Px を-10.0[Pa/m] に固定し，粒子数 N を 67，116，243，310 と変更させることで代表体積要素内における土粒子の混入数が異なるモデ(b) y2 方向流速

ルによる数値粘度測定を行った．数値実験より算定

図 2: 管路流れ解析の様子される平均流量 Q を Hagen-Poiseuille 式に代入する
ことで得られる粘性係数 µ についてまとめたグラフ
を図-1 に示す．縦軸は算出された粘性 µ を液相の粘性 µl (=0.01Pa·s) によって除した値であり，横軸は代表体積要素内に存在する土粒子体積 Vs を代表体

積要素の全領域体積 V で除した値であり，どちらも無次元量である．この図から，土粒子の占める体積率が上昇するにつれ粘性が増加する傾向が観察できる．さらに管路内領域での挙動を詳しく見るために，最も高い粘性値が算定された N =310 のケースの管路内 y1 -y2 断面における流速分布を図-2 に示す．同

図からは，平均流量に寄与しない接線方向成分 u2 ，u3 が，粒子の混入により乱れた分布を持っていることが観察できる．個々の粒子の存在が流れに対する抵抗となっているためと考えられ，これらのミクロな場における散逸が粘性として表出していることが示唆される．

(c) y3 方向流速

5終わりに
本研究では，固液混相流を均質な単一流体としの代表的な物性値であるマクロ粘性の算定が可能となる
数値材料試験手法を提案した．具体的には， Hagen-Poiseuille の法則を用いて管路内流れ解析による数
値粘性測定手法を提案し，粘性の算出を可能とした．
また土粒子の体積分率など異なる複数の解析条件を設定して数値粘度測定のケーススタディを行った．

参考文献[1] 高瀬慎介，森口周二，寺田賢二郎，小山直輝，金子賢治，車谷麻緒，加藤準治，京谷孝史：構造物の破壊を考慮した有限被覆法に基づく構造流体連成解析，
土木学会論文集 A2(応用力学), Vol. 71, No. 2（応用
力学論文集，Vol. 18），I 203–I 212，2015．
[2] Tezduyar, T. E. : Stabilized finite element formulations for incompressible flow computations, Adv.
Appl. Mech, 28, pp.1-44, 1991.

24 − Academic Center for Computing and Media Studies, Kyoto University

高プラントル数流体を用いた HZ 液柱内温度差マランゴニ対流の数値解析―特異な対流場における熱的相互作用の影響―
江口 俊和
東京理科大学理工学研究科機械工学専攻

1

緒言ばれる，特異な対流場の発現領域，及び発現メ
カニズムを明らかにすることを目的とした．ま自由界面を有する気液界面において，温度差を付与すると表面張力差が発生する．この対流た，JEREMI に向けた準備研究の一助も担っている．

はマランゴニ対流と呼ばれ，微小重力環境下やマイクロスケール下で顕在化する，自然対流の駆動力の一つとして着目されている．例えば，単結晶材料生成法のひとつである
floating-zone 法は，純度の高い単結晶材料を生成ことができる．一方で，材料生産の高効率化のために，微小重力環境下で利用するとマラン
ゴニ対流の不安定性によって生じる振動流が材料生成に悪影響を及ぼすことが知られている 1 )．故に，マランゴニ対流の不安定性に関する理解や制御の観点で研究が行われている．

Fig. 1 The schematic image of a liquid bridge

2

計算条件

私の所属する研究室は、国際宇宙ステーショ
ン日本実験棟「きぼう」にて行なわれていた，

液柱の高さを H，半径を R とし，アスペクト

MEIS( Marangoni Experiment in Space ) や比Г = H/R と定義する．今回の計算では，Г =していた。そして今後，JEREMI( Japaneseず，常にストレートな形状を保持しているもの

Dynamic Surf.と称される流体物理実験に参画
European Research Experiments on Marangoni

Instabilities )と称される，流体実験をおこなう予定である．図 1 は half-zone モデルと呼ばれる，2.0 で固定させた．液柱は表面の変形を考慮せとする．計算系全体は無重力状態であることを想定した．
流体は非圧縮性ニュートン流体であること

floating-zone 法を模擬した，最も典型的なモデを仮定し，Pr 数が 16 であるものとした．無次ルである．上下壁面間に液柱を形成し，上部を元化された流体の運動と熱輸送を記述する支加熱，下部を冷却することで液柱自由表面上に配方程式は，
Navier-Stokes 方程式，
連続の式，

温度勾配を付与し，マランゴニ対流を発生させエネルギー方程式とした．

る．JEREMI では，気体側に強制対流を付与し，

∂u*
Γ 2 *
+ u* ⋅∇ u* = −∇p* +
∇u
∂t *
Re

(2)

∇ ⋅ u* = 0(3)

自由界面での熱伝達制御をおこなうこととなっている．
本研究は，マランゴニ対流の不安定化メカニ
ズムの解明による理学的貢献，及び産業分野への波及効果を目指し，研究を実施している．本
研究は，MEIS にて観測された，X 状構造と呼

(

)

∂T *
Γ 2 *
*
+ u* ⋅∇ T
=
∇T
∂t *
Ma

(

)

(4)

また，液柱自由表面上での熱の授受を記述す京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 25

たな知見を示すことができた．また，先行研究

る境界条件は，以下で定義した．

∂T *=− Bi T * − T∞**∂r surf(*T∞=
TC* + ∆T * z *

)

(5)
(6)

Bi は自由界面での熱の授受を示す無次元数，ビオ数( Biot 数 )であり，hH/λで示される．

で見られた X 状構造は Bi が大きいところで，
かつ対流場が強い領域にて発現することが明らかとなった．一方で，Bi が小さいところでも

SW が発現し，X 状構造とは異なる新たな SWを初めて観測することができた．また，その内部構造もそれぞれ異なっている．図 5，6 はいそれぞれ，一周期分の X 状構造と新たな SW の

ただし，h : 周囲気体の熱伝達率，H : 液柱高内部構造で，赤と紫はそれぞれ時間平均温度かさ，λ : 液柱の熱伝導率，である．本研究ではらの高温，低温変動分を，水色は等温面を示し

周囲気体領域を設けず，液柱の自由表面の熱のている．X 状構造では，低温変動分のロッド上授受をビオ数でモデル化した．計算を行なう際部側が液柱自由表面に近づくと，二つに分かれは，ビオ数とマランゴニ数をパラメトリックにると，一つは液柱中心へと戻り，もう一つは液変化させて計算を行なった．特に，マランゴニ

柱自由表面に沿って下部ロッド側へと向かっ

数は，所定のマランゴニ数で十分発達させた後ていく．一方，新たな SW はそれぞれの温度変に，マランゴニ数を増大させていくという手順

動分が，入れ替わっているような挙動をしていを繰り返した．

ることが明らかとなった．SW は右回りの TW

3と左回りのTWの重ね合わせによって生じる振結果
図1 は横軸に液柱自由表面上での熱伝達を示

す無次元数，Bi を示し，縦軸にマランゴニ対流の強さを示す無次元数，マランゴニ数( Ma )を動流とされている 3)．それぞれの SW は，発現領域と内部構造の特性から，X 状構造は Bi が正側で発現する TW の重ね合わせ，新たな SW

は Bi が負側で発現する TW の重ね合わせであ示す．中抜きの点が定常流，中塗りの点が振動ることが推察される．今回の新たな SW は実現流を示し，×と実線はプロット点から予測され象でこの対流場を未だに観測していないため，る振動流遷移点，臨界マランゴニ数( Mac )であ

実験において再現されることが今後望まれる．

る．Bi が負のとき液柱全体は周囲気体から熱を受け取っており，正のときは熱が周囲気体に逃げていることを示す．Bi が負から正の方向に向
かってグラフをみてみると，Mac は減少から極小値を迎え，再度急上昇したのちに，一度 Macが落ち着く．その後，Mac は最大値を迎え，緩やかに減少していく．振動流の状態は脈動振動
流 ( standing waves ， SW ) と回転振動流( traveling waves，TW )の 2 つに分けられる．
本研究では，Mac の近傍で大きく分けて 4 つに分かれることが明らかになり，Bi の絶対値が比Fig. 2 The flow pattern map as a function of Bi.

較的小さいところでは TW，一方 Bi が比較的大きいところでは SW が安定的に発現することが明らかになった．本研究結果で見られた，
二種類の TW の類別は，先の宇宙実験でも同様の類別がされている 2)．従来の数値シミュレー
ションでは，安定的に SW を発現させることが難しく，十分発達した対流場では SW から TWへと遷移していくが，本研究は SW に関する新26 − Academic Center for Computing and Media Studies, Kyoto University

Fig. 3 The time series of X shape structure

を強くしたときに X 状構造が見られ，X 状構造の発現領域を明らかにした．Bi が負の領域では，
X 状構造とは異なる，新たな SW の振動流を発現し，内部構造を確認することができた．

5
Fig. 4 The time series of another SW.

参考文献

1) A. Cröll et al., J. Cryst. Growth, vol. 191, pp.
365–376, 1998.

2) T. Yano et al., Phys. Fluids, vol. 29, p. 44105,
2017.

3) K. Boronska and L. Tuckrman, J. Fluid Mech., vol.
559, pp. 279–298, 2006.

Fig. 5 The internal structure of X shape structure.

Fig. 6 The internal structure of another SW.

4

結言
Pr = 16 の流体で，無重力状態を想定したストレートな液柱において，自由表面上で液柱と自由表面の熱伝達を無次元数 Bi で制御したモデルにて，臨界マランゴニ数の遷移過程，X 状構造の発現領域，発現条件を調査した．Bi の大きさが比較的小さい領域では，TW の振動流が二種類見られ，その境界は，Bi が正のところで小さな値であるところとなっていた．これは，先の宇宙実験で見られた傾向に非常に類似していた．Bi の大きさが大きくなると，SW の振動流が二種類見られ，安定的に SW が発現することを数値シミュレーションで確認することができた．特に Bi が正の領域側では，対流場

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 27

カオス力学系におけるネットワーク解析
小林幹 ∗
∗

1

立正大学経済学部

はじめに
カオス力学系にはカオスアトラクタに埋め込まれた加算無限個の不安定周期軌道が存在することが知られている [1, 2]. そして, カオス軌道はこれら加算
無限個の不安定周期軌道間を不規則に遷移する現象
であると考えられている. この不安定周期軌道間の不規則な遷移現象にはどのような統計性質が潜んでいるのかを明らかにすることが本研究の目的である.
そのために, 我々はカオスの不安定周期軌道間の遷移現象をネットワークを使ってモデル化し, ネットワー
ク理論の観点からこの遷移現象を解析する. ネット
ワーク理論とは物 (ノード) と物 (ノード) との繋がを数値的に構成する. 時間 n (n = 1, 2, · · · , N ) におけるその軌道の位置を {xn } と表す.

(3) 時間 n におけるカオス軌道に沿った位置 xnからユークリッド距離の観点から最も近い不安定周期軌道を見つける. (例えば, その不安定周期軌道を

U P Oi とする.)(4) 時間 n + 1 におけるカオス軌道に沿った位置
xn+1 からユークリッド距離の観点から最も近い不安
定周期軌道を見つける. (例えば, その不安定周期軌

道を U P Oj とする.)(5) U P Oi から U P Oj へリンクを張る (もし i = j

であれば自己ループを構成する.).

り (リンク) を解析するためのもので, 自然, 生命や(6) (3) から (5) を繰り返す (n = 0, · · · , N − 1).

社会における様々な現象がネットワーク理論で解明

十分大きな数 N をもって以上の手続きを繰り返すされている [3].
ここでは, 不安定周期軌道をノードとみなし, カオ
ス軌道を用いてノード間のリンクを構築する. これにより, カオスの不安定周期軌道間の不規則遷移現象
を数理的にモデル化することが可能である. 本研究ことで, ノード数 M (これは不安定周期軌道の数に対応する) の有向ネットワークが構成される. 連続時間
カオス力学系に対してはポアンカレ写像を考えることで同様の扱いが可能である.

では, ネットワークの出次数に着目し, カオスのネッ
トワークを解析する [4].

2

ネットワークモデルの構築
時間離散カオス力学系を例にしてネットワークモ

デルの構築アルゴリズムを述べる.

(1) カオス力学系においてカオスアトラクタ (カオティックサドルでも良い) に埋め込まれた不安定
周期軌道をニュートン法などの方法を用いて数値的に見つける. 数値的に得られた不安定周期軌道を

U P Om (m = 1, 2 · · · , M ) と表すことにする.

(2) カオスアトラクタ (またはカオティックサドル)
上の適当な初期条件 x0 からカオス軌道に沿った位置

28 − Academic Center for Computing and Media Studies, Kyoto University

3

ローレンツモデルのネットワー
ク構造

ローレンツモデルにおけるカオスのネットワークを構成する: ẋ = −σx + σy, ẏ = −xz + rx − y, ż =xy − bz, ただし, ここでは σ = 10, b = 8/3 と r = 28の値を用いる. ネットワークを構成するために, ポア
ンカレ断面 z = r − 1, ż > 0 上で与えられるポアン

カレ写像を考える. ここで扱う不安定周期軌道の数

M は 1259 であり, これはポアンカレ写像周期で 14までのほぼ全ての不安定周期軌道を網羅している.
図 1 は, ポアンカレ写像周期 14 の不安定周期軌道のみに着目した, 周期と次数との関係を示している.

distribution

degree

1000

450
400
350
300
250
200
150
100
50
0

100

10

1
9.8

10

10.2

10.4

10.6

10.8

11

period110degree100

図 1: ローレンツモデルにおけるポアンカレ写像周

図 2: ローレンツモデルにおける周期 14 までの不安期 14 の不安定周期軌道の周期と各々が持つ次数との定周期軌道 1259 個を用いて構成したネットワークの関係. 周期と次数に明らかな相関は見られない.

次数分布. リンクの構成は N = 600000 の長さのカオス軌道を用いて行われている. 次数分布はべき的

図からは, 周期と次数に明らかな相関は見られない.
つまり, ネットワーク化して得られる次数という特徴な関数と見ることができるので, ネットワークがフ
ケールフリー性を有していると示唆される.

量は, 周期とは別の情報を含んでいると考えられる.
図 2 は, N = 600000 の長さのカオス軌道から構成されたネットワークの次数分布を表している. 次数分布は概ねべき的であると思われるので, ネット
ワークがスケールフリー性を持っていることが示唆される.
さらに, ローレンツモデルのネットワークの構造を明らかにするために, 平均パス長と, クラスター係数を示す. その際, ローレンツモデルのネットワークの定性的な特徴を明らかにするため, ローレンツモ
デルのネットワークと同じノード数 1259 と平均次数40.119 をもつランダムネットワークと比較する. 表
1 より, ローレンツモデルのクラスター係数はランダ
ムネットワークと比べて大きく, 平均パス長はロー
レンツモデルのネットワークの方が長いことが分かる. これにより, ローレンツモデルのネットワーク
がランダムネットワークよりは複雑な構造を持つが,いわゆる複雑ネットワークの構造と完全には一致しないことが示された.
Network
RandomLorenzProperties
Average path-length
Clustering coeﬃcient
Average path-length
Clustering coeﬃcientまとめ
本研究報告では, カオス軌道における不安定周期

軌道間の不規則遷移現象をネットワークを用いて数理モデル化する方法を示し, その方法によってネットワーク化されたローレンツモデルのネットワークの構造を示した. 結果として, ローレンツモデルのネットワークは, 次数分布がべき的であったり, クラスター係数がランダムネットワークと比べて大きいなど, 単縦なランダムネットワークとは異なる構造を持っていることが明らかになった. [4] ではさらに,次数とカオスの統計量との関連性も明らかにしており, カオスのネットワーク化がカオスを特徴付けるための方法として有益であることが示唆されている.

参考文献[1] P. Cvitanović, Physica D 51, 138 (1991).
[2] T. Kai and K. Tomita, Prog. Theor. Phys. 64,
1532 (1980).

values
1.98
0.03
2.58
0.18

表 1: ローレンツモデルのネットワークとランダム
ネットワークの構造の比較.

4[3] R. Albert and A.-L. Barabási, Rev. Mod. Phys.
74, 47 (2002); R. Albert, H. Jeong, and A.L. Barabási, Nature 406, 378 (2000); R. Albert,
H. Jeong, and A.-L. Barabási, Nature 401, 130(1999).
[4] M. U. Kobayashi, and Y. Saiki, Chaos 27, 081103(2017).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 29

Numerical investigation of the effect ofhydrogel surface shape on drag performance in turbulent flow
STOYANOVA Petya Valentinova
大阪大学大学院 基礎工学研究科

Direct numerical simulation of turbulent channel flow with varying slip velocity on thewalls was performed in order to investigate the drag reducing properties of hydrogel paints
with rough surfaces. It was found that, as long as the water trapping effect persist, surfaces
with higher mean peak-to-valley roughness can reduce drag better than smoother surfaces.

1 Introduction

flow rate in a turbulent channel flow with
infinitely-long parallel plates.
To

In maritime transportation more than halfsimulatetheeffectofhydrogelof the total resistance is due to frictional drag,roughness the velocity is considered at theresulting in great losses for the shipping

plane at the levels of the tips of the dimples ofindustry. Recently novel bioinspired hydrogel

wavy roughness such as the one shown in Fig.

paints have been developed to meet the present

1a. The following Navier slip boundary

requirements of non-toxicity. It was also foundconditions [2] were applied symmetrically tothat hydrogel coatings can reduce the frictional

the flat channel wall:

resistance while still possessing hydrodynamic
roughness. In order to take full advantage ofus  lxthat property the mechanism by whichhydrogel paints influence the turbulent dragu,
y wallws  l z

wy wall(1)where us and ws are the slip velocities in theshould be investigated. Shimada et al. [1] havestreamwise

suggested that hydrogels trap water in therespectively, and lx and lz are the slip lengths.

dimples of the rough surface. This work
presents a parametric study, where the waterandspanwisedirection,

The latter were varied via filtering functions:fx  2 
1 
x  ,
1  cos 2   x (2)channel flow.

fz  2 
1 z 
1  cos2   z (3)

2 Numerical method

l x , l z  f x f z l0

trapping effect is modeled as a local slip velocity
boundary condition on a flat wall in direct
numerical simulations (DNS) of turbulent(4)

to achieve the spatially changing slip velocity

The governing equations of the flow are(Fig. 1b) corresponding to the spatially varying

the incompressible Navier - Stokes equations

depth of the rough surface. l0 corresponds to theand the continuity equation, solved for constant

30 − Academic Center for Computing and Media Studies, Kyoto University

maximum depth of the valleys.

(a)3 Results
The drag reduction ratio (DR) was defined
as:

(b)

DR[%]  dp / dx   dp / dx 100 (5) dp / dx no slipslip

no slipwhere – dp/dx|no-slip and – dp/dx|slip are themean pressure gradients under no-slip and slip

Fig. 1 Schematic (a) the water trap effect byconditions, respectively.

rough hydrogel surface and (b) the local

In the simulations the wavelength was

slip model.
The governing equations were discretizedusing 2nd order finite difference method. For the
temporal integration of the convective andviscous terms 3rd order Runge-Kutta andCrank-Nicolson

are

used,

respectively.

A

computational domain of size 4πδ × 2δ × 2πδ
and 128 × 150 × 128 grid points was used in thestreamwise (x), wall-normal (y) and spanwise

chosen to be of the order of the turbulent
structures

in

order

to

investigate

theinteraction between the varying slip velocity
and turbulence. The effect of the different slip
conditions on the drag reduction ratio is shownin Fig. 2. For the smallest slip length l0+ = 0.9
the drag reduction is about 1%, it rises to ~ 2%for l0+ = 1.8 and varies from 3.3% up to 5.2% for

the largest slip length l0+ = 3.6. At the two small

directions (z), respectively. The simulations
were performed at Re  U c /  4200 orslip lengths the DR depends on the combinationRe  u  /  180 , where Uc is the centerline

for the same slip length the difference between

0velocity, δ is the channel half height, ν is the
kinematic viscosity and u 0 is the friction
velocity for no slip channel flow. In order toinvestigate how the size and height of thedimples affect the turbulent vortex structurenear the wall we varied the reference slip
length, and the streamwise and spanwise
wavelengths as shown in Table 1. Hereafter,
superscript “+” denotes quantities normalizedby u 0 and ν.
turbulent channel flow.

Reference slip length, l0+

two wavelength combinations can be as high as~ 1%. This correlation, however, is mostly lost
at l0+ = 3.6, where the DR increases with theincrease of slip length. A possible explanation
for this phenomenon might be the reduced
surface-vortex interaction at larger slip lengths,

then the no-slip areas act merely as obstacles
that disturb the near-wall flow. The larger
wavelength correspond to more gradual andless frequent disturbances leading to lower
friction.
Figure 3a shows the mean streamwise

Table 1. Parameters used in the DNS of

Parameterof streamwise and spanwise wavelengths andvelocity as function of the channel width for noslip and a slip case with l0+ = 3.6, λx+ = 550, λz+ =50. As a result of the surface slip pattern thevelocity is decreased everywhere in the channel.

Value
0.9 , 1.8, 3.6

Streamwise wavelength, λx+275, 368, 550, 735, 1100

Spanwise wavelength, λz+50, 183

It also leads to reduced root mean square (rms)velocity fluctuations (Fig 3b).
To investigate the effect on the vorticalstructures the isosurfaces of the second

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 31

invariant of the velocity gradient tensor:
1 u i u iQ2 x j x j

in all cases and it decreased with increasing(6)

slip length, reaching up to 5 % reduction. It was
found that for the small slip lengths the dragare plotted in Fig. 4. It can be seen from Fig. 4a

reduction depends strongly on the combination

and Fig. 4b that the slip conditions led toof wavelengths, but that the interaction

smaller number of vortices, but that the size of

between the different wavelengths is mostly

the structures was increased, i.e. under slip

lost for high slip lengths. Then the dragconditions the vortices are stabilized.

decreases with the increase in wavelength size.
Finally, it was found that the slip conditions led
to more stable vortical structures of lower
intensity. The results from this research imply
that, as long as the surface can trap water
effectively,

roughnesswithlargermean

peak-to-valley height is better at reducing drag.

No slip case

(a)

u’

Fig. 2 Dependence of the drag reduction on the

0.32

surface parameters.

(a)

y x

z

(b)

- 0.32
+

+

+

0

x

z

Slip case: l = 0.9, λ = 550 λ = 50

u’

0.32

(b)

z

y x

- 0.32

Fig. 4 Isosurfaces of the second invariant Q=
0.55 (in blue) for (a) no slip case and (b)Fig. 3 Comparison of the (a) velocity profiles

slip case l0+ = 0.9, λx+ = 550, λz+ = 50. In theand (b) root mean square (rms) velocity

plane: streamwise velocity fluctuations,fluctuations for no slip and slip case.

4 Conclusions
A

three-dimensionalu’.

References
directnumerical

simulation of turbulent channel flow withvarying local slip velocity condition was
performed to model the effect of trapped waterlayer of hydrogel coating. The drag was reduced

32 − Academic Center for Computing and Media Studies, Kyoto University

1.

N.

YamamoriandM.

Shimada,“Verification of the Reduction Effect of
New Antifouling Paint for Ship bottoms‘LF-Sea®’ on the Friction Resistance byUtilizing

ActualShips,”TECHNO-COSMOS, 22 (2009), pp. 28-33.
2.

T. Min and J. Kim, “Effects of Hydrophobic
Surface on Skin-Friction Drag”, Phys.

Fluids, 16 (2004), L55-L58.
3.

J. Kim, P. Moin and R. Moser, “Turbulence
Statistics in Fully Developed Channel
Flow at Low Reynolds Number”, J. Fluid

Mech., 177 (1987), pp. 133-166.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 33

MuSTAR MD に基づくタンパク質構造変化の自由エネルギー経路計算手法の確立と応用
山守 優 ∗∗

1

大阪大学基礎工学研究科

はじめに
タンパク質などの生体分子は、外部の摂動に伴う構造変化によって機能を発揮する。近年、計算機の発達に伴って、分子動力学 (Molecular Dynamics:MD)などの分子シミュレーション手法を用いて、生体高分子の機能メカニズム解明を行うことが強く動機づけられてきている。分子シミュレーションによる生体分子の機能解析のためには構造変化過程の再現 · メ

カニズム解明が重要である。しかしながら現在の計算資源をもってしても、通常のシミュレーション手法を用ちいただけでは、m 秒 ~ 秒のスケールにおよぶ大規模な構造変化過程を再現することは困難である。この困難に対処するために、レプリカ交換法 [2]、モデルの精度の高さを兼ね備えたサンプリング結果が得られるようになっている。筆者は、本手法を真空中のテスト系 (Alanine-dipeptide、Met-Enkephalin)に適用し、構造変化に関わる自由エネルギー変化を計算した。レプリカ交換法、レプリカ交換アンブレ
ラサンプリング法、Temperature Accelerated MD、などの既存の手法と比較してサンプリング効率と自由エネルギー計算精度の両面での優越性を確認した(図 1)。
本研究においては、1) 溶媒を陽に含む系にも適用可能なように MuSTAR MD の拡張を行うこと及び

2) CV 選択法として、機械学習の分野で発達した疎構造学習を用いた CV 選択法の導入を行った。

Temperature Accerelated Molecular Dynamics [3]、メタダイナミクス [4] などのさまざまな効率的サンプリング法が開発されてきた。
本研究においては、筆者の開発した Multi-Scale

Sampling using Temperature Accelerated and ReplicaExchange Molecular Dynamics(MuSTAR MD) [1]を拡張し、タンパク質などの生体高分子の機能メカ
ニズム解明に有用な方法として確立することを目的として研究を遂行した。MuSTAR MD は、構造変化を適切に記述可能なように少数の自由度を Collec-

tive Variable (CV) として選択し、CV 部分を独立に時間発展させる拡張系を導入している。CV 系を介して精密な全原子モデルと効率的な粗視的モデル
が相互作用するシミュレーションを行うことができる。更に、CV 系とその他のモデルとの結びつきの強さが異なる拡張系を複数個用意し、レプリカ交換によって結びつきの強さを交換することで、全ての拡
張系に渡って粗視的モデルの効率性の高さと全原子34 − Academic Center for Computing and Media Studies, Kyoto University

図 1: 左図は MuSTAR MD の概念図。右図は二面角空間に射影した真空中の Alanine-dipeptide の自由エネルギー地形。(a) MuSTAR MD、(b) レプリ
カ交換法による計算結果

2

結果および考察

列データへの異常箇所同定問題への解法として近年

2.1

MuSTAR MD の拡張

訓練データが与えられているとして、新しい標本や発達している方法である。異常箇所同定問題とは、

MuSTAR MD においては、以下のように全原子系、粗視系、CV 系の 3 種の系が同時進行で時間発展する拡張系を扱っている。

∂V α (xα , z) ∂W α (xα , z)αmα−
+ (thr, β α ) (1)i ẍi = −
∂xα
∂xα
i
i
mz z̈l = −N∑∂W α (xα , z)
α

∂xα
i

+ (thr, β z ) (2)αα
ただし、α = F G, CG1 , CG2 , . . .、mαi , xi , (thr, beta )はそれぞれの原子座標系の原子 i の質量、座標及び系の熱浴、mz , zl , (thrβ z ) は CV 系の仮想的な質量、座標、熱浴である。このときに全原子系と粗視系の計算を独立に行うように拡張した。このことによって、計算の大部分を占める全原子系の計算にオープ
ンソースで高速な GROMACS を採用でき、高速化が可能となる。現在、GROMACS を用ちいたコードを開発中である。

2.2

一点データが与えられた時に、各次元に対して元の訓練データとの異常度を計算せよという問題である。
疎構造学習においては、データの状態を変数間の関係に着目して記述する。このときに変数間の関係は、訓練データの精度行列 Λ として表現する。即ち、Λ
の (i, j) 成分が非零の値であるとき、i 番目の自由度と j 番目の自由度に相関があることを意味する。
このときに、精度行列 Λ が適切な疎行列で近似できるとき、多自由度のデータから重要な少数の自由度を抜き出したと言える。疎構造学習においては、精度行列が疎行列になりやすいような適切な事前分布の元で、事後確率最大原理に基づき、精度行列を求める。即ち、疎構造学習とは、直感的には理解不可
能な多自由度の時系列データに対して、重要な変化を引き起こした自由度を異常度つきで特定すること
で、意味のある解釈を与えるための方法であるといえる。Chignolin に対する予備計算を行った。10 残
基の de novo タンパク質 chignolin の Cα 間距離の時系列に対して、疎構造学習を適用した。その結果、

CV 選択法の導入

TYR2-TRP9 間距離と PRO4-GLY7 間距離に強い

MuSTAR MD の大規模系への適用にあたって重相関を持つことが検出できた。

要になるのは、系の変化を記述する適切な少数の自由度を CV として選択することである。現状では、対象系への化学的・生物学的な知見に基づいて、重要な自由度を予測し、反応座標とみなすなどの手段をとる。系に対する事前の知識を前提とすることな
く、重要な自由度を得ることを可能にする解析方法
として、主成分解析、がある。主成分解析は、原子座標のトラジェクトリに対して直交変換を行い、互いに相関のない座標の線形結合（主成分）を得る手法

3

結論
本研究では、テスト系において、既存の手法と比較してサンプリング効率と自由エネルギー計算精度において勝った効率的サンプリング手法 MuSTAR

MD を、1) 溶媒を含む系にも適用可能 2) 疎構造学
習を用いた CV 選択法の導入という二点において拡張するために研究を遂行した。

で、多くの研究で用いられている。この方法は、構
造ゆらぎの記述に適しているが、構造変化の引き金になる短時間の変化の抽出には適していない。他の方法も、適用可能な自由度数の限界や抽出したモードを物理的にどのように解釈するかが自明でないな
どの課題点を持つ。時系列データから少数の重要な相関を持つ変数の組み合わせを同定し、その組み合
わせの時系列に沿った変化をもとに、変化に関係する自由度を抽出する方法である疎構造学習を、MDに適用して、新しい解析方法を模索した。疎構造学
習とは、機械学習の分野において、多自由度の時系

参考文献[1] Y. Yamamori, A. Kitao, Jour. Chem. Phys., 139,
145105, (2013).
[2] Y. Sugita, Y. Okamoto, Chem. Phys. Lett., 314,141-151, (1999).
[3] L. Maragliano, E. Vanden-Eijnden, Chem. Phys.
Lett., 426, 168-175, (2006).
[4] A. Laio, M. Parrinello, PNAS, 99, 12562-12566,(2002).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 35

破壊力学に基づく損傷モデルによる鉄筋コンクリートの3 次元破壊シミュレーション
相馬 悠人
茨城大学 理工学研究科 都市システム工学専攻

1

はじめに
社会資本への投資削減により，土木構造物の設計・維持管理の合理性が要求されている．それに伴い，建設材料である鉄筋コンクリート(以下，RC)の力学挙動をより詳細に予測や把握することへの重要性が増している．
RC 部材はコンクリートの破壊特性および鉄筋とコンクリートの相互作用により，大小様々なひび割れが進展する複雑な破壊挙動を示す．数値解析による RC 部材の破壊挙動を再現する試みは盛んに行われているが，RC 部材の破壊挙動を 3 次
元で詳細に再現した例はほとんどない．
対して車谷らは，鉄筋とコンクリートの力学挙

図-1 領域分割法による並列有限要素解析

動を別々にモデル化し，鉄筋の幾何形状まで反映させた有限要素モデルを用いることで，RC 部材とコンクリートの力学挙動のモデル化が重要となの破壊挙動を詳細に再現できる数値解析手法を開る．本研究では，車谷らが開発した数値解析手法発した[1]．しかし，この方法は，鉄筋の幾何形状を用いる[1]．この手法では，鉄筋の材料モデルにを有限要素モデルで再現するため，細かなメッシ

非線形等方硬化則に基づくvon-Mises 塑性モデル，

ュが必要となる．
複雑な配筋を有する実規模のRC

コンクリートの材料モデルに破壊力学に基づく等部材への適用を考えた場合，内部構造を再現する

方性損傷モデルを適用している．

のに，大規模な有限要素モデルが必要となる．
そこで本研究では，RC 部材の数値解析手法に
領域分割法に基づくプロセス並列化を適用し，大規模シミュレーションを実行できる方法を開発する．そして，開発した数値解析手法の計算効率を検証し，実規模の RC 部材への適用性を示す．

2

RC 部材の数値解析手法

2.1 鉄筋とコンクリートの材料モデル
RC 部材の破壊挙動を再現するためには，鉄筋

36 − Academic Center for Computing and Media Studies, Kyoto University

2.2 並列化
RC 部材の数値解析手法では，連立一次方程式の解法に対角スケーリングのみを施した共役勾配法を用いている．本研究では，主に連立一次方程式の解法の部分に，領域分割法に基づくプロセス
並列化を適用している．
数値解析の手順を図-1 に示す．手順としては，
入力データの読み込み後，前処理として，対象とする有限要素モデルの分割数の設定，分割領域の割当，節点と要素番号の割当，境界条件の付け替

図-2 RC はりの有限要素モデル

図-3 計算時間とコア数の関係えを行う．次に，分割領域毎にコアを割り当て，
プロセス間で情報の送受信を行いながら並列有限要素解析を実行する．最後に，後処理として，分割領域毎の解析結果を集約し，有限要素モデル全体の解析結果を出力する．

3

大規模シミュレーションへの適用性

3.1 計算効率の検証
大規模シミュレーションへの適用性を示すため，開発した数値解析手法の計算効率を検証する．対象とする有限要素モデルは図-2 に示す RC はりである．対称性を考慮し，解析範囲は 1/4 とした．

図-4 3 次元でのひび割れ分布

一次の四面体要素を用い，節点数と要素数は135473 と 750042 とした．強制変位 12 mm を 600
step で与えた．この有限要素モデルにおいて，使用するコア数を 8 (2×2×2)，16 (4×2×2)，32 (8×2×2)，64 (16×2×2)とし，計算時間を比較した．
カッコ内は，x, y, z 軸方向の分割数を示している．
図-3 に計算時間とコア数の関係を示す．使用するコア数の増加とともに，計算時間が短縮されており，並列化による計算効率の向上が確認できた．
以上より，開発した数値解析手法は，実規模の RC
部材における大規模な有限要素モデルにおいても高速で解析が可能である．

4

おわりに実規模の RC 部材の詳細な破壊挙動を再現する

ため，RC 部材の数値解析手法に領域分割法に基づくプロセス並列化を適用した．RC はりの有限要素モデルを対象とし，コア数と計算時間を比較した結果，計算効率の向上が確認でき，実規模のRC 部材の数値解析において必要となる大規模な有限要素モデルにおいても，高速に数値解析を実行可能であることを示した．

参考文献

3.2 解析結果
3 次元でのひび割れ分布を図-4 に示す．この可[1] 車谷麻緒，根本優輝，相馬悠人，寺田賢二郎：

視化結果は，コンクリートの材料モデルに適用しコンクリートの破壊力学を考慮した鉄筋コン

た損傷モデルにおける損傷変数の結果を 3 次元でクリートの 3 次元破壊シミュレーションとそ

可視化したものである．3 次元で複雑に進展するの性能評価，日本計算工学会論文集，Vol.2016,

ひび割れを詳細に再現できていることが分かる．

pp.20160004, 2016.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 37

分子性導電・磁性材料の基底状態および励起状態に関する理論的解析
中野 義明
京都大学大学院理学研究科化学専攻

1 緒言

制御する一方で、TPC のような有機分子がイオンラ
ジカル塩の結晶構造を制御する役割を果たす。TPC

金属フタロシアニン（MPc）は、特異な光学・導電・分子は、有機イオンを収容するのに適した六角形の磁気物性を示すことが期待され、光学、電子、光電空隙を形成し、多成分錯体において層状構造を提子素子のセンサーや材料として利用されている。例

供する。具体的には、TiOPc、VOPc のラジカルアニ

えば[M L2Pc] （M = Co, Fe; L = CN, Cl, Br）から作

オンを含む同形錯体 1、2 の結晶は以下のように得III−製された部分酸化状態の塩は、半導体的、または、
金属的導電性を示すようになる。また、た。Me4PBr が過剰に存在するジクロロベンゼン中で

MIVOPc をナトリウムフルオレノンケチルで還元した。

[FeIII(CN)2Pc]−から成る塩は巨大磁気抵抗を示す。

さらにベンゾニトリル、TPC を加えた後、ヘキサンを一方、金属フタロシアニンは還元されることによって反応液にゆっくりと拡散させることにより目的の単結も、金属状態や超伝導状態になることが理論的に予晶を得た。

測されている。しかしながら、還元状態の金属フタロ

X 線構造解析の結果、錯体 1、2 は同形であり、シアニンは大気中で不安定であるため、結晶構造の{[MIVO(Pc•3−)]•−}2 ダイマー層と Me4P+-TPC-C6H4Cl2

報告が少なく、構造－物性の相関が明らかになって層が交互に積層した層状構造となっていることが分いない。

かった。また、結晶構造を基に計算した錯体 1、2 に本研究では、多成分法によって作製した 2 種類の同 形 錯 体 (Me4P+)[MIVO(Pc•3−)]•−(TPC)0.5·C6H4Cl2

おける{[MIVO(Pc•3−)]•−}2 ダイマー間の重なり積分は、

ダイマー内の重なり積分の 10 分の 1 以下であり、ダ（M = Ti (1), V (2)）の結晶構造、および、光学的・磁イマー間の相互作用は極めて小さく、各ダイマーは気的特性を調べた。ここで、TPC はトリプチセンであ孤立していることが分かった。

る（Figure 1）。

錯体 1、2 の UV-vis-NIR スペクトルを測定したところ、Pc•3−に特徴的な吸収帯が観測される一方で、0.85 eV 付近にも吸収帯が観測された。また、磁化率を測定したところ、[(TiIVOPc•3−)•−]2 において J/kB =−123 K の反強磁性相互作用、また、[(VIVOPc•3−)•−]2において J1/kB = −105.0 K、J2/kB = −15.2 K の反強磁性的相互作用が見積もられた。
Figure 1. Molecular structures of MIVOPc (M = Ti, V)(left) and TPC (right).

2 結果と考察
多成分法は金属フタロシアニンやフラーレンのラそこで錯体 1、2 の電子状態を明らかにすべく、

CAM-B3LYP/cc-pVTZ/cc-pVDZ レベルで、C4v 対称性の[MOPc]0、C2v 対称性の[MOPc]•−、[(MOPc)•−]2（M = Ti, V）について理論的解析を行った。[TiOPc]0、[TiOPc]•−、[(TiOPc)•−]2 のフロンティア軌道、電荷・ス

ピン密度を比較することにより、[(TiOPc)•−]2 の半占

ジカルアニオン塩を作製するために開発された合成（SO）軌道は、[TiOPc]0 の Pc 部位に広がった 2 重縮手法である。この方法では、有機カチオンが電荷を退した最低空（LU）軌道に由来することが分かった。

38 − Academic Center for Computing and Media Studies, Kyoto Universityすなわち、[(TiIVOPc•3−)•−]2 という電子状態を支持する（Figure 2）。

分子内の V−Pc 間には弱い反強磁性的相互作用が働いていることが示唆された。これにより、実験的に求められた J1、J2 がそれぞれ分子間、および分子内反強磁性的相互作用に由来することが分かった。

Figure 2. (a) Calculated [(TiOPc)•−]2 dimer. Monomers

are highlighted in red and blue. (b) Calculated spin
density. (c) Energy diagram for the frontier Kohn–
Sham orbitals of the 1A state. The lines in red and blue

are ascribed to the orbital levels mainly contributed bymonomers in red and blue, respectively. HO, SO, and

LU in red and blue indicate the highest occupied, singly

occupied, and the lowest unoccupied orbitals ofmonomers shown in red and blue, respectively. The
α-orbitals are shown on the right side.

次に、0.85 eV 付近の吸収帯の起源を明らかにする た め に 同 レ ベ ル で [TiOPc] 、 [TiOPc] 、
0•−[(TiOPc)•−]2 の TD-DFT 計算を行い、実測スペクトルとの比較を行った。その結果、Q 帯の分裂等の分子内遷移を 再現で き 、 0.85 eV 付近の 吸収帯が[TiOPc] 分子間の SO → SO/LU 遷移による電荷移

Figure 3. (a) The observed UV–vis-NIR spectrum of 1

in KBr. Calculated spectra of (b) the 1A state in[(TiOPc)•−]2 dimer, (c) the 2B1 state in [TiOPc]•−monomer, and (d) the 1A1 state in [TiOPc]0 monomer.

The SO → SO, SO → LU, and HO → LU chargetransfer transitions are indicated as red lines.

3 結言

•−

動吸収帯であることが分かった（Figure 3）。[VOPc]0、(Me4P+)[MIVO(Pc•3−)]•−(TPC)0.5·C6H4Cl2 の 単結晶[VOPc]•−、[(VOPc)•−]2 についても同様に電荷移動吸作製に成功し、結晶構造を明らかにした。また、光収帯と帰属できた。一方、低スピン状態と高スピン状学・磁気的特性を実験・理論的に調べ、錯体中では態の解析から、[(TiOPc) ]2 について、J = −196 K のPc•3−（S = 1/2）、TiIV（S = 0）、VIV（S = 1/2）の電子状態•−

分子間反強磁性的相互作用、[VOPc]•−について J =をとっていること、Pc−Pc 間、Pc−M 間の相互作用に−23 K の分子内反強磁性相互作用が見積もられた。

より物性が発現していることを明らかにした。今後は、

すなわち、2 量体[(VOPc) ]2 において、分子間にはより広範な金属フタロシアニン化合物について検討•−

Pc の積層により強い反強磁性的相互作用が働き、し、一般性を導出したいと考えている。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 39

ポスト「京」100 ナノ電子状態計算むけの
超並列構造保存波束ダイナミクスソルバー開発
井町 宏人
鳥取大学

1 はじめにに対するエネルギー保存連続段 Runge-Kutta 法を取り上げる。これは自然な並列性を持つという有機物質の量子伝導シミュレーションのある定長所があり、その高並列実装を検討する。

式化では時間依存非線形 Schrödinger 方程式𝑖𝑖𝜓𝜓̇(𝑡𝑡) = 𝐻𝐻 (𝑡𝑡)𝜓𝜓(𝑡𝑡)

並列アルゴリズム設計の際には一般的に、計算対の数値解が必要になる。ここで、𝐻𝐻はハミルトニなど）にどう割り当てるかが重要な選択になる。

アンと呼ばれる行列である。次元𝑛𝑛は離散化の手

例えば、非常に大まかな見積もりであるが、密行法によるが、原子に局在する電子軌道の線形結合

列アルゴリズムでは 1 ノードあたり最低 1000×で解を表す LCAO 近似では、原子数に比例したオ

1000 次元程度の部分行列を一度に扱わなければーダーになる。有機高分子集合体においては構造

並列化のオーバーヘッドが並列化による速度向上のもつ乱雑さが伝導度の決定において本質的であを上回ってしまう。この割り当ての選択を変更するため、必然的に原子数は大きくなり、現実の電ると実装も大きな影響を受けるため、各並列化方

子デバイスのスケールの系を扱おうとした場合

向のサイズを事前にある程度想定しておくことが、

108 のオーダーに達する。このような系を扱うた

高い実効効率の並列プログラムを素早く実装する

where 𝐻𝐻 (𝑡𝑡) ≔ 𝐻𝐻0 − 𝛾𝛾 diag(|𝜓𝜓1 , … , |𝜓𝜓𝑛𝑛
𝐻𝐻0 ∈ 𝑅𝑅

𝑛𝑛×𝑛𝑛

|2

|2 )

,

, 𝐻𝐻0𝑇𝑇 = 𝐻𝐻0 , 𝜓𝜓(𝑡𝑡) ∈ 𝐶𝐶 𝑛𝑛 , 𝛾𝛾 ∈ 𝑅𝑅.

象に存在する並列化方向（行列の行方向 / 列方向など）のそれぞれを、計算機アーキテクチャに存在する並列化方向（ノード内並列 / ノード間並列

めには、その特性を生かした並列性の高い数値解ために重要である。本稿では京都大学スーパーコ法が必要になる。

ンピュータシステムを含む多くのスーパーコンピ

時間依存非線形 Schrödinger 方程式は、ハミルトュータで利用可能な、ノード間を MPI による分ン系という著しい特性を持つ時間発展方程式のク

散並列化、ノード内を OpenMP でスレッド並列

ラスに属する。その特性の一つはエネルギー保存

化するハイブリッド並列パラダイムを前提にする。

性であり、この場合次で定義されるエネルギー
1𝛾𝛾
𝐸𝐸 (𝑡𝑡) ≔ 𝜓𝜓(𝑡𝑡)† 𝐻𝐻0 𝜓𝜓(𝑡𝑡) − ∑|𝜓𝜓𝑖𝑖 (𝑡𝑡)|42
42 ハミルトン系に対するエネルギー保存連続段 Runge-Kutta 法

が時間非依存、すなわち

𝑖𝑖

文献[2]では、連続段 Runge-Kutta 法と呼ばれ𝑑𝑑𝑑𝑑 (𝑡𝑡)=0𝑑𝑑𝑑𝑑

る常微分方程式の数値解法のクラスについて、そであることが示せる。エネルギー保存性やシンプれが離散版のエネルギー保存性を満たすためのレクティック性といった、系が本来持つ構造を離(実質的な)必要十分条件、次数条件、および後述

散化後も再現させるような数値スキーム（の導出の並列化のための条件を明らかにし、これらを満法）は構造保存数値解法と総称され、定性的によたす 6 次精度までの数値スキームを実際に構成しい数値解を生成することが知られている[1]。本稿ている。導出の詳細は省くが、𝑁𝑁次元ハミルトン

では特に文献[2]で提案されている、ハミルトン系

系に対する𝑠𝑠段スキームの計算手順は、ある𝑠𝑠𝑠𝑠次40 − Academic Center for Computing and Media Studies, Kyoto University

の非線形方程式に帰着され、さらにその解法に簡点である。全体のコア数が105 を超える近年の最易ニュートン法を採用すると、最終的に𝑠𝑠個の独大規模のスーパーコンピュータにおいても、行列

立な𝑁𝑁次元連立一次方程式

次元𝑁𝑁が使用する計算機システム全体のスレッド

(𝑘𝑘)
(𝑘𝑘)
𝐴𝐴𝑖𝑖 𝑥𝑥𝑖𝑖 = 𝑏𝑏𝑖𝑖 (𝑖𝑖 = 1, … , 𝑠𝑠, 𝐴𝐴𝑖𝑖 ∈ 𝑅𝑅𝑁𝑁×𝑁𝑁 ),数𝑝𝑝よりも十分大きいと仮定できるので、各𝑚𝑚𝑖𝑖 の𝐴𝐴𝑖𝑖 = 𝐼𝐼𝑁𝑁 − ℎ𝜆𝜆𝑖𝑖 𝐽𝐽0

計算を 1 スレッド内で完結させても十分に並列性

い実数、𝐽𝐽0 は現在時刻におけるヤコビ行列である。

コストがほとんどかからない。前処理の適用がある。これらの独立な連立一次方程式は並列に解よる疎行列ベクトル積である点も、実装の容易さが現れる。ここで、𝑘𝑘は簡易ニュートン法のステを利用できることになる。これは計算結果の分配

ップ数、ℎは時間刻み幅、𝜆𝜆𝑖𝑖 は方程式に依存しなといったわずかな後処理部分を除いて並列実装の非線形 Schrödinger 方程式においては𝑁𝑁 = 2𝑛𝑛でGMRES 法本体で使用するのと同次元の行列にくことができる。本研究ではこの部分の並列化をという点で有利である。

MPI のグループ機能により実装した。行列𝐴𝐴𝑖𝑖 の疎第二の利点は次回の連立一次方程式を解く際に再性はハミルトニアン𝐻𝐻のものを引き継いでいる。

利用が可能な点である。簡易ニュートン法の反復

多くの場合、原子間相互作用の空間的な局所性にでは係数行列は一定であるから反復間で前処理行

より、これは疎になる。次節ではこれらの連立一列𝑀𝑀を再利用できるのに加え、時間方向に連続的次方程式の効率的な並列解法を検討する。

な現象の時間発展を扱っているので、次の時間ス

3 近似逆行列前処理付き GMRES 法

テップ𝑡𝑡𝑖𝑖+1 における係数行列𝐴𝐴は今の時間ステッ

前処理行列𝑀𝑀は𝑡𝑡𝑖𝑖+1 でも有効と考えられる。こうした再利用による総計算コストの低減は多項式近非対称連立一次方程式
𝐴𝐴𝐴𝐴 = 𝑏𝑏, 𝐴𝐴 ∈ 𝑅𝑅

プ𝑡𝑡𝑖𝑖 におけるそれと近く、結果として𝑡𝑡𝑖𝑖 で計算した

𝑁𝑁×𝑁𝑁

𝑇𝑇

𝐴𝐴 ≠ 𝐴𝐴似などの一部の前処理にはない優れた性質である。

の数値解法のうち、ここでは残差最小化の原理に基づく疎行列反復解法である GMRES 法を取り

4 おわりに上げる。疎行列反復解法で最も基本的かつ実行時間を占める操作は疎行列ベクトル積である。ここ

本稿では時間依存非線形 Schrodinger 方程式にでは詳細は省くが、GMRES 法、疎行列ベクトル

適用可能な数値手法である、ハミルトン系に対す積ともに高効率な並列化の研究およびパッケージ

るエネルギー保存連続段 Runge-Kutta 法を取り化された実装が多数存在する。

上げ、その並列実装を与え、その際に考慮すべき反復解法では一般に前処理点を論じた。これらの論点は時間依存非線形
Schrodinger 方程式以外にも広く通じるものであ

𝑀𝑀𝑀𝑀𝑀𝑀 = 𝑀𝑀𝑀𝑀の選択が非常に重要であるとされる。ここで𝑀𝑀はる。

前処理後の係数行列𝑀𝑀𝑀𝑀の条件数が小さくなり、<謝辞> 本共同研究制度を活用させて頂きました

かつ𝑀𝑀自身と行列ベクトル積 𝑦𝑦 ↦ 𝑀𝑀𝑀𝑀 の計算コことに感謝致します。

ストも小さいような行列である。𝑀𝑀は陽に計算さ[1] E. Hairer, C. Lubich and G. Wanner,

れるとは限らない。本稿では前処理として、𝑀𝑀の

Geometric

各列を、連立一次方程式

Structure-Preserving Algorithms for Ordinary
Differential

𝐴𝐴𝐴𝐴𝑖𝑖 = 𝑒𝑒𝑖𝑖 (𝑖𝑖 = 1, … , 𝑁𝑁)

Numerical
Equations,Integration:
2ndの荒い近似解𝑚𝑚𝑖𝑖 とする近似逆行列前処理を採用

Springer-Verlag, Heidelberg, 2006.

ed.,する。ここで𝑒𝑒𝑖𝑖 は𝑖𝑖番目の要素のみが 1 で他が 0 で[2] Y. Miyatake and J. C. Butcher, A

ある単位ベクトルである。荒い近似解の生成にはCharacterization

疎行列疎ベクトル積を用いた GMRES 法を数反

Methods and the Construction of Parallel

復で停止させるという手法を取った。

Integrators for Hamiltonian Systems, SIAM

近似逆行列前処理は本稿の問題設定において次のJournal on Numerical Analysis Vol. 54, Issue 3,

二点で有利である。第一は容易に並列化が可能な

pp. 1993-2013, 2016.

ofEnergy-Preserving

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 41

固相光化学反応の計算化学による解明
鈴木聡
京都大学 福井謙一記念研究センター

1 はじめに2 理論と計算方法

本研究では固相光化学反応の機構を調べるため

光化学反応では、熱平衡状態の分子が吸光し、の方法論の整備を行った。具体的には、励起状態

分子構造はそのままで電子状態のみが励起されたでの最小エネルギー経路(minimum energy path,

フランクコンドン状態（以下、FC 状態）が生成

以下 MEP)を計算するため、string 法による構造

する。FC 状態から構造変化が起こり、蛍光や非最適化を実装した。また、ONIOM 法を用いた固断熱遷移が起こる。蛍光は励起状態の極小、非断相中での円錐交差構造最適化プログラムを実装し熱遷移は状態間交差の周辺で起こりやすい。そこた。本稿では前者について中心に述べる。

で、励起状態の PES 上で、FC 状態、極小点、TS、

熱化学反応の機構を考える際、反応物と生成物

状態間交差を結ぶ MEP を計算し、熱化学反応とをつなぐ、経路に沿った極大点すなわち遷移状態

同様の議論を行いたい。問題は、FC 状態や状態(transition state, 以下 TS)を議論するのが通常の間交差は極小値ではないため熱化学反応で行われ手順である。典型的な手続きでは、まず、反応物

るような IRC 計算を行ってもそれらは得られなと生成物の間の中間的な分子構造で Hessian をい、ということである。TS から IRC 計算を行う

計算し、最寄りの極大値を得る。TS は一次の鞍以外の方法として Nudged elastic band (NEB)法点である必要があるので Hessian はちょうど一や string 法による MEP 最適化が存在する。これ

つ、負の固有値を持つ必要がある。TS で のらは反応経路そのものを初期推測から最適化してHessian の固有ベクトルの方角に沿って少しずつゆく方法である。これらの方法であれば端点を固

ポテンシャルエネルギー面（以下 PES）を降りて定した条件下での MEP を計算することができる。

行く IRC 計算を行うことで極小値に到達し、反応

String 法では以下のような最適化サイクルに物と生成物を結ぶ経路を得ることができる。

より MEP が得られる。

一方、光化学反応では、熱平衡状態の分子が吸① 反応経路の端点を結ぶ経路上の N 点に対し各光し、分子構造はそのままで電子状態のみが励起点でエネルギー勾配と経路の接ベクトルを計されたフランクコンドン状態（以下、FC 状態）

算するが生成する。FC 状態から構造変化が起こり、蛍② 勾配を接ベクトルに直交する成分に射影し、光や非断熱遷移が起こる。非断熱遷移は円錐交差

その射影に沿って各点のエネルギーを最小化の周辺でおこりやすい。しかしながら、FC 状態

するように各点を更新するや円錐交差は非平衡点であるから、熱化学反応と同様の方法では MEP の計算ができない。そこで
反応経路全体を最適化するstring法を用いた実装③ 更新された点をスプライン補間することで経路全体を補間する④ 経路長を等分するように新たに N 点を取り直により、励起状態での MEP 計算を行えるようにし、接ベクトルも更新するした。

このような更新を繰り返すことで、勾配の接ベ
クトルに直交する成分が 0 に収束し、MEP が得られる。ここで、この iteration の間で Hessianは必ずしも計算する必要がないのがポイントであ42 − Academic Center for Computing and Media Studies, Kyoto Universityる。円錐交差付近では多参照電子状態理論を用いは蛍光量子収率が低く、凝集により蛍光量子収率る必要があるが、多参照電子状態理論の Hessian

が 向 上 す る こ と が 知 ら れ て い る [J. Org.

を実装していないプログラムが多く、Hessian がChem., 2017, 82 ,pp 6865]。得られた円錐交差はアル不要であることは大きなメリットとなる。本研究キルアミノ基の回転とピレン環の面外変角により

では②のステップで各点のエネルギーの和を大きな構造変化を伴うので、凝集により嵩高いア

L-BFGS 法で更新する方法[Sheppard, JCP, 128,

ルキルアミノ基の回転が抑制され、結果的に内部

134106]を用いた。典型的な場合には端点も最適転換が阻害されるものと考えられる。

化し極小点を探すが、本研究では必要に応じて端点は固定することにした。

3 結果と考察
テスト計算として、ベンゼンの S1 状態の極小構
造と S1/S0 状態間の最小エネルギー円錐交差の間に存在する TS を求めた。まず S1 状態の極小構造と円錐交差をそれぞれ計算し、その間を 13 点で内挿した。エネルギーと勾配は 2 状態平均(6,6)CASSCF/6-31G で計算した。String 法による最適化は 10 回でほぼ収束し、極小値に比べ105.0kJ/mol 高い近似 TS を得た。この近似 TSから Hessian を用いた Netwon-Raphson 法を行ったところわずか 3 回の iteration で TS を得ることができた。string 法で TS の十分そばまで収束しているためである。
これを踏まえて、4,5-bis ジアルキルアミノピレンの S1 状態の極小構造と S1/S0 状態間の最小エネルギー円錐交差の間に存在する TS を求めた。円錐交差は事前に求めておき、円錐交差から基底状
態の極小までの MEP を初期推測とし、S1 状態の
PES 上で string 法による最適化を行なった。ここで、円錐交差から基底状態の極小までの MEPを初期推測とする理由は、簡便に得られる初期推測であり化学的に自然な（不自然に結合長が短いような構造を取らない）経路だからである。このような経路を初期推測とすれば比較的簡単に励起状態の TS を求めることができる。この結果、この分子では、S1 面の極小点からみてまずアミノ基の回転が起こり 60ｋJ/mol の高さの TS に到達し、一度構造が緩和したのち、ピレン環の面外変革により 120ｋJ/mol 程度の円錐交差に到達することがわかった。これらのエネルギーは FC 状態より低いため、光励起された 4,5-bis ジアルキルアミ
ノピレンは円錐交差に到達し内部転換により失活
していくものと考えられる。この分子は溶液中で京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 43

漂流物の衝突力に関するばらつきの基礎的検討
小谷 拓磨 ∗ ，森口 周二 † ，高瀬 慎介 ‡ ，寺田 賢二郎 †∗東北大学 工学研究科 土木工学専攻†‡

東北大学 災害科学国際研究所

八戸工業大学 工学部 土木建築工学科

1 はじめに津波による構造物の被害を考える上では，遡上した津波の波圧によって構造物が損壊する一次被害に加え，津波漂流物の衝突による二次被害も考慮するべきである．近年では，漂流物の衝突を考慮した防波施設や津波避難ビル，原子力発電所などの沿岸域の重要構造物の耐力設計に関するガイドライン [1] や，漂流物の衝突力想定に含まれる不確実性を考慮し確率論的にリスク評価を行おうとする試みがある [2]．
これらに用いられる予測式では，漂流については解析的なアプローチ [3] が，衝突力については実験的なアプローチ [4] がなされている傾向にある．しかの漂流物を構造物の衝撃面から 400mm の位置に 1つ配置した Case-2, Case-3, Case-4 によって，漂流物の質量の変化に伴う流体力と衝突力への影響を検証する．加えて，構造物の衝撃面から 100mm の位置に500kg の漂流物を 1 つ配置した Case-5，Case-5 の漂流物の後方 100mm に同質量の漂流物を加え計 2 つを配置した Case-6，Case-5 の漂流物の側方各 50mmに同質量の漂流物を加え計 3 つを配置した Case-7 を検討することで，漂流物の量と配置の違いによる影響を検証した．なお本解析には，剛体と流体の相互作用を表現し得る手法として，有限被覆法に基づいた
流体構造連成手法である中村ら [5] の手法を用いた．

しながら，現状の津波避難ビル等の沿岸構造物の耐力設計には，3 次元的な漂流物同士の相互作用による漂流挙動によるばらつきや，漂流物の量による衝
突力のばらつきを考慮した検討は十分ではない．そこで本稿では，複数の漂流物による複雑な漂流挙動
を３次元的に考慮した衝突力に関する不確実性のパ
ラメトリックスタディによる定量化のための基礎研究として，衝突力と各不確実性との関係を定性的に評価し，検討すべき不確実性の推定を行う．本稿ではその一例として，漂流物の質量・数・配置位置に関する衝突力の不確実性を仮定する．

3

解析結果

3.1

質量の変化による作用力の変化

Case-1～4 における構造物への流体力および漂流物の衝突力を図-2 に示す．Case-1 では漂流物が存在しないため流体力のみプロットしている．まず Case-1～Case-3 を比較すると，漂流物の質量が増加するにしたがって流体力は減少し，衝突力が上昇している傾向にある．これは質量が増加したために衝突力が増加したこと加え，流体力が漂流物を加速に必要な力

2 解析モデルへ変換され流体力が減少したものと考えられる．一方 Case-4 では，Case-2 および Case-3 に比べて流体力

本研究における解析モデルを図-1 に示す．固定さが卓越し，衝突力は小さくなった．すなわち，Case2れた構造物の前に剛体を仮定した漂流物を配置し，ダおよび Case3 とは逆の傾向となった．これは図-3 にムブレーク問題を解くことで，構造物の衝撃面におけ示すように，段波先端で漂流物が加速した Case-2,

る流体力と衝撃力を得る．解析モデルは，まず漂流物

Case-3 に比べ，Case-4 では漂流物の加速に時間を要を配置しない Case-1 と， 500kg， 1000kg， 1500kg

した為に段波のみが流れ方向に先行し，衝突時には

44 − Academic Center for Computing and Media Studies, Kyoto University

的な衝突力が発生しており，Case-5 と Case-6 との間

Pressure-sensitive surface(Top view)

でピーク値の差は見られなかった．一方，Case-5 とCase-7 を比較すると，全漂流物が同時に衝突するこ

Debris blockxとにより，衝突力のピーク値がおおよそ 3 倍程度に

400

100

1000

y

上昇していることがわかる．これらの結果から，瓦

100

1220

1170

(Side view)

よる衝突までの時間差についても十分に考慮しなけ

160 670

ればならないことが分かった．

Unit: mm

550

160

4

z
x
dtFEM = 0.005(sec)

160

100

450

礫の数は衝突力に支配的であるが，それらの配置に

dtDEM =

まとめ
本研究では，漂流物の衝突力に関する不確実性のdt FEM(sec)
3000

定量化に向け，漂流物の質量・数・位置が異なる解

810,857 tet. elements

析モデルを用いて定性的な評価を行った．結果として以下の 3 点を得た．(i) 漂流物の質量の変化による衝突力の傾向は３次元的な挙動を考慮した定量化がCase-1: Without debris

400

必要であり，多くの計算ケースを検討すべきである．

100(ii) 衝突力に対して支配的なのは瓦礫の量であるが，傾向が単純であり，定量化には多くの計算ケースを

Case-2: 500kg/m� x1

400

要さないと考えられる．(iii) 衝突力に対する瓦礫のCase-5: 500kg/m� x1

100

位置に関しては，衝突時の時間差を考慮すべきである．今後，更なる漂流物に関する不確実性の検討を

100
Case-3: 1000kg/m� x1

400

行い，考慮すべき不確実性の選定を進めていく．

Case-6: 500kg/m� x2

100
50
50

Case-4: 1500kg/m� x1

参考文献

Case-7: 500kg/m� x3

図 1: 解析モデル
先に構造物に到達した波先による鉛直上向きの力と，先行段波による衝突時の緩衝効果が作用したためと考えられる．このことから，質量の変化に対する衝撃力の傾向は単純な傾向ではないことが分かった．3
次元的な挙動も考慮した定量化には，より多くの解析ケースを実施する事が必要である．

3.2 数・配置の違いによる作用力の変化
Case-5～7 の解析結果を図-4 に示す．Case-5 と Case6 を比較すると，漂流物の増加に伴って流体力が減少していることが見て取れる．これは，流体力を受ける漂流物が増加したことによる明示的な結果と言える．しかし衝突力に関しては，後方漂流物による断続[1] Federal Emergency Management Agency: Guidelines
for design of structures for vertical evacuation fromtsunamis, 2nd edition (FEMA P-646), prepared by Applied Technology Council for Federal Emergency Management Agency, 2012.
[2] 吉田郁政, 武内大記, 大森政則, 藤井直樹, 柳沢賢：津
波漂流物による衝突力ハザードの算定方法について,
土木学会地震工学論文集, Vol. 29, pp. 912–917, 2007.
[3] 木原直人, 松山昌史, 藤井直樹：漂流物挙動解析による津波漂流物衝突に関する確率論的評価手法, 土木学会論文集 B2（海岸工学）, Vol. 69, No. 2, pp. 341–345,
2013.
[4] Ikeno, M., Takabatake, D., Kihara, N., Kaida, H., Miyagawa, Y. and Shibayama, A.: Improvement of collision force formula for woody debris by airborne andhydraulic experiments, Coast. Eng. J. Vol. 58, No. 4,
1640022, 2016.
[5] 中村正人，高瀬慎介，樫山和男，寺田賢二郎，車谷麻緒：有限被覆法に基づく自由表面を有する流体–構造
連成解析手法の構築, 土木学会論文集 A2（応用力学）,
Vol. 67, No. 2（応用力学論文集 Vol. 14）, I 199-I 208,2011.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 45஺ణ੔ഈ峑ၝ૴

Fluid force้఍&DVH 岴峲峝 &DVH峕岴岻峵้఍ৎ峘଱৿Collision force

150๗ઉ্਱峢峘ৡ岶੿৷

100
50
0
1.0

&DVH峕岴岻峵้఍ৎ峘଱৿

図 3: Case-2～Case-4 における衝突時の挙動の違い

50

500

0

0
0.0

0.5
time [sec]1.0

Fluid force [kN]

200

1500

150
1000

100
50

500

0

0
0.0

0.5
time [sec]1.0Collision Force [kN]

Case-2: 500kg/m� x1

200

1500

150
1000

100
50

500

0

0
0.0

0.5
time [sec]1.0Collision Force [kN]

Case-3: 1000kg/m� x1

Case-4: 1500kg/m� x1

図 2: 解析結果（Case-1～Case-4）

Collision force

200

1500

150
1000

100

500

50
0

0
0.0

0.5
time [sec]1.0Case-5: 500kg/m� x1

200

1500

150
1000

100

500

50
0

0
0.0

0.5
time [sec]1.0Case-6: 500kg/m� x2

200

1500

150
1000

100

500

50
0

0
0.0

0.5
time [sec]1.0Case-7: 500kg/m� x3

図 4: 解析結果（Case-5～Case-7）

46 − Academic Center for Computing and Media Studies, Kyoto University

Collision Force [kN]100

Fluid force

Collision Force [kN]

1000

Collision Force [kN]

1500

150

Collision Force [kN]

Fluid force [kN]200

Fluid force [kN]

0.5time [sec]Fluid force [kN]0.0

Case-1: Without debris

Fluid force [kN]

 ੔ষ஺ణ峕峲峵೹้஍ટ஺ణ峘੔ষ

Fluid force [kN]

Fluid force [kN]200

気象条件に応じた都市の境界層乱流の組織構造に関する研究
吉田
∗

1

敏哉 ∗

京都大学大学院理学研究科地球惑星科学科専攻

はじめに

2

観測との比較

人口構造物により被覆された都市地表面は気象環

本研究では計算領域として京都市内の南北 11 km・境を大きく変えるインパクトをもつため、ヒートア

東西 2 km・高さ約 1 km の領域を再現し使用した。

イランドや大気汚染といった都市特有の環境問題を

この計算領域は大気乱流観測が行われた京都大学防引き起こす。これらの問題解決には、運動量・スカ

災研究所宇治川オープンラボラトリーを含むようにラー輸送を担う乱流に対して都市構造物が及ぼす影設定しており、観測データと LES の結果を比較する

響を明らかにすることが重要である。

ことが可能である。観測は鉄塔に取り付けられた超大気境界層内の乱流は空間的に組織だった構造（乱流組織構造）をしており、運動量やスカラーの大部

音波風速計とその鉄塔付近に設置したドップラーラ
イダーによって行われた。

分はこの構造に伴い輸送される。近年、観測や数値

Fig. 1a に観測で得られた４事例（D1-D4）と LES

シミュレーションによって、都市構造物が立ち並んの平均主流風速の比較を示す。図は両対数軸で示しだような粗面上にはストリーク構造やヘアピン渦とた。LES と各観測事例のプロファイルはおおよそいった構造が存在することが示されてきた [1][2][3]。

1/4 の傾きに沿っていることが分かる。これは平均

一方、実在する都市の構造物は非一様性が大きく、風速が高度のべき乗に比例していることを示してお特に地面に近い高度において流れ場に与えるその影
響は強まる。例えば Park et al. [4] は高層ビルの後り、そのべき指数が 1/4 程度であることを意味する。
Counihan[5] によると、都市および郊外地域におけ

流域に形成された局所的な渦構造の存在を示しておるべき指数は 0.21 から 0.28 の範囲に入るため、観り、都市構造物の幾何的な特徴と乱流組織構造の関

測・LES 共に観測地点におけるべき指数を正しく示係が示唆される。

している。

そこで本研究では、建物近傍の詳細な流れ場を再Fig. 1b に高度 25 m における主流風速成分のパ現可能な Large-eddy simulation (LES) を用いて、実ワースペクトルを示す。ここでは LES と超音波風

在都市における乱流組織構造の解析を行った。特に速計のプロファイルの比較を行っている。LES のプ

実在都市の構造物の特徴の１つである建物高さのばロファイルは観測データと共に乱流スペクトルに現らつきの影響を評価した。また、地面近傍では機械れる慣性小領域での-5/3 乗則（ここでは規格化によ的な作用による乱流生成が卓越するため、本研究でり-2/3 乗）が確認できるため、LES は風速の変動成

は中立成層の気象条件に応じた組織構造を対象とす

分をよく再現できているといえる。

る。以下では、LES の精度を検証するため行った実以上の結果より、本研究で使用した LES は中立成大気の観測結果との比較と、建物高さのばらつきの層下における都市構造物上の乱流をよく再現できる

影響を評価するために行った感度実験の結果を示す。

ということが示された。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 47

フラックスに対応する第 2 象限の ejection 及び第４象限の sweep は組織構造に伴って発生することが知

(a)
200

D1 (ane)D2 (ane)D3 (ane)D4 (ane)D1 (dop)D2 (dop)D3 (dop)D4 (dop)
LES (CTL)

150

z [m]

100

50
1/4

られており、これらの大きさを調べることにより乱流構造の特徴を把握することが可能である。ここでは、sweep と ejection の差 ∆S0 を用いて、ejectionと sweep のどちらが卓越した流れ場であるかを調べた。Fig. 2b は空間平均された ∆S0 の鉛直プロファイルである。CTL と UNI の ∆S0 は高度とともに減少している。これは、建物が多く存在する地面近く

(a)
10

2

20

3 4 5 67
-1

10

-2/3

2.5

–u [m s ]100

z/ hallfE(f)/σu2

(b)

1

-1

10

D1 (ane)D2 (ane)D3 (ane)D4 (ane)KaimalLES (CTL)

-2

10

-3

10 -3
10

-2

10

-1

10

fz/ u–

1
0.5

0

10

10

1

Fig 1: (a) Vertical proﬁle of mean streamwise velocity in the observations and the LES. (b) Power

(b)

CTL
UNI
0.1
0
0.001 0.002 0.003——
2-<u'w'>all / U∞2010

建物高さのばらつきの影響

2.5

3

本研究では京都市の構造物を再現した計算 (CTL)
及び、その構造物全てを計算領域内の平均建物高さ(hall ) に揃えた場合の計算 (UNI) を実行した。Fig.
2a は計算領域内で空間平均されたレイノルズ応力の鉛直プロファイルである。図から hall 以下では CTLと UNI の差はほとんどないことが分かる。一方、hall
以上では CTL の値が UNI より大きい傾向にあり、ピークの高度は CTL のほうが高いことが示された。

z/ hall

spectra of streamwise component obtained from thesonic anemometer and the LES at 25 m.

1
0.5
CTL
UNI
0.1
-0.4 -0.2 0

0.2 0.4<∆S0>all

続いて、乱流組織構造への影響を調べるため運動量フラックスに関する４象限解析を行った。4 象限解
析は、水平風速の変動成分 u′ と鉛直風速の変動成分

w′ の符号に基づいて鉛直運動量フラックス u′ w′ を
４象限に分類し、各象限の値が運動量輸送に寄与す
る割合を調べる手法である。ここで、下向き運動量48 − Academic Center for Computing and Media Studies, Kyoto University

Fig 2: Vertical proﬁles of (a) Reynolds stress normalized by free stream velocity U∞ and (b) ∆S0
averaged spatially over the analysis region.

の高度では sweep が卓越するが、高度が上がるにつ向が見られる。この特徴は Fig. 3b のレイノルズ応れ ejection の寄与が増加し、約 2.5hall を超えたあた

力と一致する。よって、λp が 0.3 以上でみられたレ

りで ejection が卓越するということを示している。

イノルズ応力の増加は sweep により寄与されるとい

また、ほとんどの高度で CTL の値が UNI より大きうことが示唆された。

いため、建物高さのばらつきは sweep の寄与を増加させる傾向があるということが分かる。
次に建物高さのばらつきの影響を定量的に示すため、計算領域を 1 km 四方の領域に区切り、領域毎に乱流統計量と粗度パラメータ λp (単位面積に含まれる建物の割合) を計算した。Fig. 3 はレイノルズ
応力及び ∆S0 と λp の関係を示している。ここではFig. 2 の結果に基づき、0.5hall と 2.5hall の高度のプロファイルを調べた。Fig. 3a から 0.5hall の高度
では CTL と UNI 間でのレイノルズ応力の差はほとんどないことが分かる。一方、Fig. 3b より 2.5hall
では λp が約 0.3 を超えたところで CTL と UNI の差が現れる。これは建物高さのばらつきの効果は λp が
約 0.3 を超えたときに顕在化することを意味している。続いて、Fig. 3c と 3d に ∆S0 と λp の関係を示す。図から 0.5hall の高度では CTL と UNI における

∆S0 はほとんど λp に依存していないことが分かる。
一方、2.5hall の高度では、CTL の ∆S0 は λp が 0.3を超えたところで λp の増加とともに大きくなる傾4

まとめ
本研究では都市構造物の高さのばらつきと乱流組織構造の関係を調べるため、京都市の建物を再現した LES を実行した。計算領域は乱流観測が行われた施設を含むように設定し LES と観測の比較を行った。これにより、本研究で使用した LES は観測結果をよく再現できることを確認した。運動量フラック
スに関して４象限解析を行うことにより、組織構造の特徴を示すパラメータ (∆S0 ) を算出し、このパラメータに対して建物高さのばらつきがどのように影響するかを調べた。その結果、組織構造に伴う下向き運動量輸送である ejection と sweep のうち、建物
高さのばらつきは sweep を増加させる傾向があることを示した。さらに、この傾向は建物の密度に対応する粗度パラメータ λp が大きい領域で顕著になることを示した。

参考文献
(a)

(b)
0.004

CTL
UNI

——
2
-<u'w'>1 km2 /U∞——-<u'w'>1 km2 /U∞2

0.003

0.003

0.002

0.002

0.001

0.001

0

0
0

0.1 0.2 0.3 0.4 0.5 0.60λp

0.1 0.2 0.3 0.4 0.5 0.6

λp

(c)

(d)
CTL
UNI

0.4
0.2
0
-0.2
-0.4

CTL
UNI

0.4

<∆S0> 1 km2<∆S0> 1 km2[1] Kanda, M, ”Large-eddy simulations on the eﬀects
of surface geometry of building arrays on turbulent organized structures,” Boundary-Layer Meteorol 118 (2006), pp. 151-168

0.004
CTL
UNI

[3] Yagi, A et al, ”Nature of streaky structures observed with a Doppler Lidar,” Boundary-Layer
Meteorol 163 (2017), pp. 19-40[4] Park, Seung-Bu ”Large-eddy simulation of turbulent ﬂow in a densely built-up urban area.” Environ Fluid Mech 15 (2015), pp. 235-250.

0.2
0
-0.2
-0.4

0

[2] Coceal, O et al, ”Structure of turbulent ﬂow overregular arrays of cubical roughness,” J Fluid Mech
589 (207), pp. 375-409

0.1 0.2 0.3 0.4 0.5 0.6

λp0

0.1 0.2 0.3 0.4 0.5 0.6

λp[5] Counihan, J, ”Adiabatic atmospheric boundary
layers: a review and analysis of data from theperiod 18801972,” Atmos Environ 9 (1975), pp.
871-905

Fig 3: Variations of Reynolds stress normalized by
U∞ at (a) 0.5hall and (b) 2.5hall , and ∆S0 at (c)0.5hall and (d) 2.5hall with λp .

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 49

海色衛星データを用いた沿岸および陸棚海域における河川由来物質の解析手法の開発
小林志保1 中田聡史 2
1 京都大学フィールド科学教育研究センター

2 神戸大学海事科学研究科

1 はじめにShimadzu）を用いて溶存有機炭素（DOC）濃度を測定するとともに，分光光度計（U-2900,

本研究課題では，近年公開された高頻度・高解Hitachi）を用いて有色溶存有機物の吸光係数

像度の衛星データ(COMS/GOCI)を用いて沿岸お（aCDOM）を測定した．現場で得られた DOC

よび陸棚海域における河川由来物質の挙動を解析

濃度と aCDOM の関係，現場の aCDOM と海色する手法を開発した．申請者らは平成 27 年 8 月

衛星から得られる aCDOM の関係，および塩分とと 11 月において，大阪湾における河川由来の溶aCDOM の関係を調べ，河川由来の溶存有機物濃

存有機物濃度と衛星から得られる海色データとの度を推定するモデルを作成した．溶存有機物濃度

間に一定の関係がみられることを見出しているの動態を調べるため，メディアセンターの可視化（Kobayashi et al.,2017）．平成 28 年には，対象

アプリケーション等を利用して，定点データの抽領域を瀬戸内海沿岸から陸棚海域に広げ，数年分

出・可視化，作図等を行なった．

の衛星データを用いて河川由来物質の挙動を解析
淀
川

大阪湾

しモデル化する．本研究で使用する衛星データは500 m の水平解像度，１時間に１回の時間解像度

播磨灘であり，100 km 四方の広さがある対象領域のデ

大和川ータを数年分解析するためには汎用 PC のメモリ
ー容量・計算能力では到底不可能である．そこで吉野川

今回スーパーコンピュータを用いて対象領域のデータを解析した．

2 方法

紀の川

紀伊水道

那
賀
川

太平洋（紀伊水道外域）

平成 28 年 10 月から平成 29 年 3 月（11 月中旬
から 2 月中旬までの繁忙期を除く）の期間中に，黒潮

瀬戸内海沿岸から陸棚海域の毎時の衛星データ

図 1 現地観測データの位置図．背景色は水深を(COMS/GOCI)をダウンロードし，スーパーコン
ピュータにアップロードした．一方，平成 28 年 5

表す．

月，8 月，11 月に船舶により大阪湾から紀伊水道3 結果と考察の 69 地点（図 1）において現地観測を行ない，表層水を採取するとともにその塩分を測定した．採取した水を濾過し，全炭素分析計（ TOC-L,

50 − Academic Center for Computing and Media Studies, Kyoto University

現場採水によって得られた aCDOM と DOC
濃度との間の相関係数は高く（r = 0.95），aCDOM

の分布から溶存有機物の分布がわかることが確認から陸棚海域の毎日・毎時の河川由来物質の分布された．また塩分と aCDOM との間の相関係数もが明らかになった．河川由来の溶存有機物につい高く（r = −0.80），塩分の低い河口域を由来とすては，近年内湾における水質環境管理を困難にす

る溶存有機物の分布は，aCDOM の分布とほぼ一る項目として，沿岸海域から陸棚海域におけるそ

致することが示唆された．また海色衛星から得らの挙動を解明する必要性が高まっている．本研究れる aCDOM と現場の aCDOM との間の相関関の成果は，河川由来の溶存有機物が陸棚海域の水係は有意であった（r = 0.55 ~ 0.78）．

質に及ぼす影響の解明につながり，内湾の水質環

以上のことから，衛星データを用いて河川，河境管理に寄与する．また，河川由来の溶存無機物

口域から陸棚海域に広がる溶存有機物の挙動を調（栄養塩）については，沿岸で養殖されている海

べた．
図2 に，
2014 年7 月25 日9 時の衛星𝑎𝑎𝐶𝐶𝐷𝐷𝑂𝑂𝑀𝑀

藻の生産に大きな影響を及ぼすことから，その挙分布図を示した．このとき紀伊水道外域を通過す

動を解明することには水産上の意義がある．本研る黒潮は，紀伊半島の潮岬から近い距離にあり，究の成果は，海藻養殖の適地および最適なタイミ

紀伊水道の東側に aCDOM の低い外洋系水が流ングを提案するために用いることができる．

入していた．一方，aCDOM，すなわち溶存有機物濃度は流量の大きい淀川が流入する大阪湾奥で

謝辞

最も高く，
大阪湾の東岸に沿って高かった．
また，
紀伊水道北部から西部にかけてもその濃度は高く，本研究における衛星データの解析は，平成 28

紀伊水道西岸に沿って比較的高濃度の溶存有機物

年度スーパーコンピュータ共同研究制度を利用しが紀伊水道外域に流出していることが伺われた．

て行いました．また，現地データの収集は，徳島

このような分布は，河川の流量が平水流量に近く，

県立農林水産総合技術支援センター 水産研究所，

黒潮が接岸している際の典型的な例として頻繁に大阪府立環境農林水産総合研究所，和歌山県水産

見られたが，一方で台風や大雨等の気象擾乱の後試験場，京都大学農学研究科海洋生物環境学分野には，１日ごとに aCDOM の分布が大きく変化すの皆様のご協力の下に行なわれました．心からのることが衛星データを用いて示された．

感謝を申し上げます．

引用文献
Kobayashi

S,

NakadaS,

Nakajima

M,

Yamamoto K, Akiyama S, Fuchi M, Hayashi M,
Ishizaka J; Visualization of the Distribution of
Dissolved Organic Matter in Osaka Bay Usinga Satellite Ocean Color Sensor (COMS/GOCI),
Journal

of

Water

and

Environment

Technology,2, 55–64, 2017.

図 2 海色衛星から得られた aCDOM の分布

4 まとめ
本研究により，衛星データの公開された平成 23
年から現在までの晴天日における，瀬戸内海沿岸

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 51

硫酸ランタン La2(SO4)3 結晶中における水拡散チャネルの発見
豊浦 和明
京都大学大学院工学研究科材料工学専攻

本研究では，100 ~ 250 ºC の低温域において高速脱水・水和反応が報告されている硫酸ランタン La2(SO4)3
中における水の存在状態および拡散メカニズムについて，第一原理計算に基づく理論解析を行った [1]．その結果，この結晶格子中の H2O は水分子の形状を保ったまま存在し，その移動過程においても分子構造を維持していることが明らかとなった．また，MD シミュレーションで確認された水拡散の素過程についてエネルギ
ープロファイルを評価した結果，b 軸方向に沿った 1 次元拡散チャネルが存在することを明らかにした．一方，
この拡散チャネルに沿った水の拡散係数は水和反応が進行する 250 ºC で 10-10 cm2/s 程度であり，拡散距離は
高々1 μm であった．このことから，本系における高速脱水・水和挙動は，ここで見出された 1 次元拡散チャネルの存在だけでなく，反応時の結晶粒微細化や高速粒界拡散にも起因していることが示唆された．

1 緒言

面波打ち切りエネルギーは 400 eV とした．水の存在状態計算に対してはユニットセルを 1×2×2

100 ~ 250 ºC の工場排熱の有効利用において，倍に拡張したスーパーセル (272 原子)，拡散挙動

化学反応により蓄熱・放熱が可能な化学蓄熱技術

解析にはプリミティブセルを 2×2×1 倍に拡張しの開発が期待されている．近年，我々の研究グル

たセル(272 原子)を用いた．また，k-point サンプ

ープでは，この温度域における脱水・水和反応に

リングは点のみとした．

より吸熱・発熱を繰り返すことのできる硫酸ラン
タン La2(SO4)3 を見出した [2]．この化合物は，他3 結果と考察の化学蓄熱材と比較して，脱水・水和反応が可逆
的かつ速やかに進行するという特長を有している．

まず，H2O の存在状態について，網羅的な構造

本研究では，高速脱水・水和反応の起源を探るこ

最適化計算によって見出された H2O の安定構造とを念頭に，第一原理計算を用いて La2(SO4)3 結を図 1 に示す．H2O の OH 結合は SO4 四面体ユ晶中における水の存在状態と拡散メカニズムを明

ニットを構成する近接 O イオン (ONN)を向いてらかにすることを目的とした．

おり (ONNH 距離：~ 2 Å)，水素結合の形成によりさらに安定化することも明らかとなった．

2 計算手法まず，La2(SO4)3 結晶中における水の存在状態について，水分子の様々な位置・配向に対する網
羅的な構造最適化計算を行った．次に、水の拡散挙動を定性的に把握するために，分子動力学(MD) 法に基づく拡散シミュレーションを行い，確認された水拡散の各素過程 (経路) に沿ったエ
ネルギープロファイルを算出した．これらの計算
は，すべて平面波基底 PAW 法 (VASP code) に基づいて行った．交換相関項には PBE _GGA，平52 − Academic Center for Computing and Media Studies, Kyoto University

図 1 La2(SO4)3 結晶中に見出された水の最安定
サイトの局所構造．

次に，図 1 に示す水の最安定構造を始状態として，1500 K で 100 ps の MD シミュレーションを行った．その結果，結晶内の等価な最安定サイト(図 2a 参照) を行き来する様子が見られたが，その移動過程においても，H2O は水分子の構造を維持していることが確認された．
図 2b および 2c に，MD シミュレーション中に確認された隣接サイト
間を繋ぐ移動経路を模式的に示す．最も近接した
サイト間を繋ぐ path 1 は a 軸方向の短距離経路

であり，path 2 および path 3 は，それぞれ b 軸
方向にジグザグに繋がる経路である．第一原理計
算を用いてこれら 3 つの経路に対するポテンシャ
ル障壁を評価した結果，それぞれ，0.27, 0.81, お
よび 1.50 eV であった．したがって，この結晶中においては，path 1 および path 2 で構成される 1

図 3 path 1 および path 2 で構成される b 軸方次元チャネル (図 3) に沿った水拡散が支配的で向に沿った水の 1 次元チャネル．

あることが明らかとなった．
この結果を踏まえて，実験的に水和反応が開始
する 250 ºC における水の拡散係数 DH2O を見積もったところ，10-10 cm2/s 程度であることがわかった．この値から見積もられる水の拡散距離は高々1 μm 程度 (t4 結言
本研究では，La2(SO4)3 の高速脱水・水和反応~ 102 s) であることを考えると，の起源を明らかにすることを念頭に，第一原理計

La2(SO4)3 の高速脱水・水和反応は，本研究で明算に基づく理論解析により結晶中の水拡散を定量らかとなった水の 1 次元チャネルの存在だけでな評価した．その結果，結晶中の H2O は水分子のく，反応時の結晶粒微細化や高速粒界拡散に起因形状を保ったまま存在し，その移動過程においてすることが示唆される．今後は，実験グループとも分子構造を維持していることが明らかとなった．

連携して，本結晶の微細組織を明らかにすること

また，MD シミュレーションで見出された移動経を計画している．

路のエネルギープロファイルを評価した結果，b軸方向に沿った 1 次元水拡散チャネルの存在を見出した．ただ，その移動に対するポテンシャル障
壁は 0.81 eV と比較的高く，水和開始温度における拡散距離は高々1 μm 程度であった．このことから，本系における高速脱水・水和挙動は，本研究で明らかとなった 1 次元拡散チャネルの存在に加え，反応時の結晶粒微細化や高速粒界拡散にも起因している可能性が示唆された．

図 2 (a) La2(SO4)3 結晶中に存在する等価なすべ[1] K. Toyoura, H. Tai, N. Hatada, K. Shizume,ての水の最安定サイト．(b)(c) MD シミュレーシ

T. Uda, J. mater. Chem. A, submitted.

ョン中に見出された移動経路（paths 1 ~ 3）．

[2] N. Hatada, K. Shizume, T. Uda, Adv. Mater.
29, 1606569 (2017).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 53

粒子との接触を伴う液体挙動の直接数値解析
自由界面―固体物体間相互作用を考慮した多相直接数値解析井上 幹允
東京理科大学大学院理工学研究科機械工学専攻

1 緒言
生活の中でもよく見られる流体現象である「濡れ」に注目する．濡れ現象とは液体が固体面に接触している現象である濡れは身近に見られる現象のみならず，工業的にも広く見られる現象である．
近年，この現象を理解し応用・制御する動きが活発化している．例えば，テフロン加工が施されたフライパンや車のフロントガラスは撥水を促進する
もの，一方，化粧品やコンタクトレンズなどは保水性を高めるため親水性を高めたものが作られている．また，製鉄工場において金属の精錬時に発生す
る煙霧の浄化効率向上のため，濡れたフィルターを用いている例(1)もある．その上で重要なものが，固体との接触を伴った液体の挙動である．この液体と固体の接触は非常に単純な現象と思えるが，液体表面の変形，固液気接触界線の移動，接触で生じる液体‐固体間の相互作用など様々な物理現象が複合したものであるため，極めて複雑なものとなる．界面の挙動に着目した研究は従来広く行われている．
Glasheen & McMahon(2)は円形ディスクを液体に押し当てた時の界面の挙動に着目した実験を行った．また，Deen et al. (3)は液面への固体粒子の自由落下に関する数値解析を行った．界面の変形および粒子挙動について解析はされているが，しかし，粒子周りの濡れ性を考慮したものはいまだ少ない．そこで本研
究では，液面への固体粒子の自由落下における，固液接触角（濡れ性）が粒子周りの流体および粒子そのものの挙動に及ぼす影響に注目し，volume of fluid(VOF)法(4)と immersed boundary (IB)法(5)を組み合わせた直接数値解析（direct numerical simulation, DNS）を行った．

54 − Academic Center for Computing and Media Studies, Kyoto University

2 計算方法
2.1 計算対象
今回解析対象として，等温環境下の容器内に静置した液面に対する，固体粒子の自由落下の系を用いた．図 1 に解析系を示す．3 次元の矩形容器内に液体を静置し，液面近傍から粒子を自由落下させる．
容器壁と液体との接触角は 90º とした．

2.2 計算手法
𝝏𝝏𝝏𝝏𝒊𝒊
𝝏𝝏𝝏𝝏𝒊𝒊

𝝏𝝏𝝏𝝏𝒊𝒊

= 𝟎𝟎

(1)
𝝏𝝏𝝏𝝏

𝝏𝝏𝝏𝝏

𝟏𝟏

𝝏𝝏𝟐𝟐 𝒖𝒖+ 𝒖𝒖𝒋𝒋 𝝏𝝏𝝏𝝏 𝒊𝒊 = − 𝝏𝝏𝝏𝝏𝒊𝒊 + 𝐑𝐑𝐑𝐑 𝝏𝝏𝝏𝝏 𝝏𝝏𝝏𝝏𝒊𝒊
𝝏𝝏𝝏𝝏
𝒋𝒋

𝒊𝒊

𝟏𝟏

𝒋𝒋

𝒋𝒋

𝜿𝜿𝜿𝜿

+ 𝐅𝐅𝐅𝐅 𝟐𝟐 𝜹𝜹𝟐𝟐𝟐𝟐 + 𝐖𝐖𝐖𝐖 + 𝒇𝒇𝐩𝐩𝒊𝒊 (2)式(1)は連続の式，式(2)はナビエストークス方程式である．代表長さ L = Dp/2，代表速度 U，代表時間 L/U，
代表圧力U2 として無次元化を行った．ui, xi, pi は各軸方向の速度成分，距離成分，圧力成分を示し，t は
時間, は界面の曲率である．式(2)中の Re はレイノ
ルズ数，Fr はフルード数，We はウェーバー数であり ， そ れ ぞ れ Re = 𝜌𝜌𝜌𝜌𝜌𝜌/𝜇𝜇, Fr = 𝑈𝑈/√𝑔𝑔𝑔𝑔, We =𝜌𝜌𝑈𝑈 2 𝐿𝐿/𝜎𝜎と表される．ここでは表面張力係数を表す．本研究では外力項として重力項 (1/Fr2)の他に二つ加えた．一つは continuum surface force (CSF)法(6)により表面張力(F = n) を考慮した表面張力項(n/We)，もう一つは固体を流体場として解くことで生じる流体‐固体間の相互作用力を考慮した相互作用力項 (fpi = up-uft)である．ここで n は界面に対する法線ベクトル，は格子内に含まれる粒子

の割合，up は粒子速度，uf は流体速度を表す．表 1

ラフを示す．この結果より，粒子上部と下部においに計算条件および解析に用いた初期条件を示す．

て著しい圧力の変化が生じていることがわかった．
この粒子上部と下部での圧力の違いによって，粒子速度に変化が出たと考えられる．これより，粒子周囲流体の圧力が固体-液体間の相互作用に対して大きく影響を与えていることが推察される．また，t =2.4 において顕著な圧力変化が見られるが，この瞬間の液体形状に着目すると，粒子上部において左右から濡れ上がってきた液体が互いに接触する瞬間ということがわかった．このことから，液体の接触の瞬間の液体形状は粒子挙動に最も大きく寄与することが考えられる．

Fig. 1 Computational domain of the impact of a sphericalparticle upon a liquid in a pool.
Table 1 Computational conditions.
Fig. 2 Temporal variation of the particle velocity at eachcontact angle.

3 結果と考察本研究では界面形状，特に固体粒子と液体形状のFig. 3 Pressure distribution around the particle at each time.

動的変化とそれに伴う圧力変化に注目した．まず，粒子と液体との接触角をパラメータとして行った

4 結言

粒子の液面への自由落下を考えた直接数値解析で
は，液面接触前後の粒子の落下速度の変化をみた．

本研究では直接数値解析を用いて，静置した液体

その結果を図 2 に示す．この結果より，液面形状がに対して粒子を落下させるという数値解析を行っ大きく変化する t = 1.5 以降において大きく違いが現た．主に，粒子落下挙動および液体形状に対して，れた．この結果から，粒子に対する接触角が粒子挙

粒子の液体に対する濡れ性，粒子周辺圧力がどのよ

動に対して大きく影響を及ぼすことがわかった．

うに影響を与えるかを検証した．まず，接触角の影

次に粒子周辺の圧力分布に着目した．図 3 に各時響を見たところ，粒子の液体に対する濡れ性つまり刻における粒子周辺圧力の分布をプロットしたグ

接触角によって粒子落下挙動が変化することを確

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 55

認した．次に，粒子周辺圧力に着目したところ，粒子上部と下部において圧力差が生じていることを確認した．そこでこの圧力差によって大きく粒子挙動が左右されることを推察した．
したがって今後としては，粒子挙動と圧力差の因
果関係についてさらに追究していこうと考えている．また，濡れ性によって粒子周辺圧力がどのように変化するかの検証も平行して行っていこうと考えている．

5 参考文献[1] Agranovski, I. E. and Whitcombe, J. M., Chemical
Engineering & Technology, 24-5 (2001), 513.
[2] Glasheen, J. W. and McMahon, T. A., Physics of
Fluids, 8.8 (1996), 2078.
[3] Deen, N. G., van Sint Annaland, M. and Kuipers, J. A.
M., Fifth International Conference on CFD in the Process
Industries CSIRO, Melbourne, Australia, 13-15 December(2006).

[4] Brackbill, J. U., Kothe, D. B. and Zemach, C., Journalof Computational Physics, 100-2 (1992), 335.

[5] 酒井幹夫編著, 粉体の数値シミュレーション,(2012), 153.

[6] Hirt, C. W. and Nichols, B. D., Journal ofComputational Physics, 39-1 (1981), 201.

謝辞：本研究を進めるにあたり，京都大学学術情報
メディアセンターの方々には多大なご協力を賜りました．ここに記し，感謝の意を表します．

56 − Academic Center for Computing and Media Studies, Kyoto University高精度視覚質感記憶の心理学的基盤と神経機構の解明

藤道 宗人
京都大学大学院人間・環境学研究科共生人間学専攻

1 問題・目的歳（SD=2.8）であった。全ての参加者は正常な視力を有していた。実験は京都大学こころの未来研

ヒトは様々な物体に囲まれて生活している。したがってヒトが安全に生活していくためには，こ究センターの MRI スキャナー及び関連実験設備を用いて行った。

れらの物体に対して適切に働きかける必要がある。
物体の正確な認知はどのようなメカニズムに支えられているのだろうか。物体の形状や色の認知メ
カニズムについては多くの知見が存在するものの，物体表面の性質，すなわち質感の認知メカニズムについては研究が進み始めたところである。これまでヒトの質感認知メカニズムについて，知覚研
究から，ヒトが物体の素材を瞬時に特定できることや腹側高次視覚野で質感を知覚処理していることが明らかにされてきた[1]。しかしながら，視覚記憶研究はほとんど行われていない。知覚は経験

2.2 刺激Ward（1992）の反射モデルのパラメータ[2]をFerwerda（2001）の変換式[3]を用いて変化させることで，粗さと光沢感がそれぞれ知覚的に等間隔に連続的に変化する球体画像セットを用意した。
球体画像にはキッチンあるいはビル街の照明環境が照射されていた。

2.3 手続きと学習無しに成立し得ないことを考えると，ヒト

実験は粗さ記憶課題と光沢感記憶課題から成っの質感認知メカニズムを解明するためには，視覚ていた。したがって，参加者は継時呈示される 2

記憶についても検討する必要がある。また，知覚枚のサンプル球体画像を見て，どちらの課題を行から視覚記憶に至る過程において，質感の視覚表っているかを識別する必要があった。例えば 2 枚

象にどのような差異と共通性が存在するのかを解のサンプルが粗さの次元で変化していれば，粗さ明する必要があるが，そのためには脳の情報表現記憶課題であった。2 枚のサンプル球体画像の呈を探ることが最も根本的である。したがって，本示後に，保持しておくべきサンプルを指示した。

研究では質感の視覚記憶メカニズムを担う脳領域

実験参加者は，試行の最後に呈示されるテスト球を解明することを目的とした。

体画像と保持した画像を比較して，粗さの程度が小さい，または光沢感の強い方の画像を選択する

2 方法

ことが求められた。照明環境に依存しない質感の視覚記憶を検討するために，2 枚のサンプルとテ

2.1 参加者

スト球体画像の照明環境は常に異なっていた。

24 名の京都大学生および大学院生（男性 15 名，
女性 9 名）が実験に参加した。平均年齢は 22.9

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 57

3 解析
3.1 マルチボクセルパターン解析（MVPA）
 質感の視覚記憶情報を表現する脳領域を同定するために，本研究ではマルチボクセルパターン解析（MVPA）を適用した。質感の見えは照明環境
や形状，輝度ヒストグラムの歪度など，多次元の情報に影響を受けるため[4][5]，脳活動パターンの規則性を抽出する MVPA が有効な手法であると考えたのである。
本研究では，全脳に対して MVPA を行うサーチライト MVPA を実施し，質感視覚情報を保持する脳領域を探索した。その際，特に腹側高次視覚野と初期視覚野，そして頭頂間溝に着目した。

関連の脳活動が生じるまでのラグを考えると，課
題開始 6 秒前後は記憶課題が粗さと光沢感のどちらであるかの弁別に対応していると考えられる。
一方でその後の時間からチャンスレベル以上で課題の判別をすることはできなかった。
初期視覚野において，課題開始 6～12 秒後の全ての時間で課題を判別することはできなかった。
また，頭頂間溝において，課題開始 6～12 秒後の全ての時間で課題を判別することに成功した。これらの結果から，質感の視覚記憶課題の情報が初期視覚野ではなく，腹側高次視覚野と頭頂間溝で表現されていることや，これらの領域の質感視覚
記憶への関与は機能的に異なる可能性が示唆された。

これは，従来の質感知覚研究から質感知覚に腹側
高次視覚野が関連すること[1]，そして視覚性短期記憶研究から初期視覚野や頭頂間溝が記憶保持に関与すること[6]が示されているからである。
さらに本研究では課題のタイムコースに沿ってMVPA を実施し，着目した各脳領域の MVPA の正答率の変化を検討した。課題の中で参加者はまずどちらの記憶課題なのかを弁別して，その後は指示されたサンプルを保持しなければならない。
したがって，弁別と保持という異なる課題要求によって関連領域に違いがあるか検討することができる。
図 MVPA の正答率の推移。網掛けになっているのは考3.2 解析手順
 遅延期間中の脳活動に MVPA を適用した。実験の結果得られたある脳活動パターンが，粗さの察対象の時間（6～12 秒）。下図はサンプルと手がかり，
そしてテストが呈示されたタイミングを表している。

5 参考文献記憶課題実行中の脳活動パターンなのか，光沢感の記憶課題実行中のパターンなのかを判別するデ[1] Wada, A., et al. (2014). Neuroimage, 98, 243–

コーダを作成した。
そして，作成したデコーダに，

257.

デコーダ作成時には使用していない脳活動パター[2] Ward, G., J. (1992). Comput Graph, 26(2),ンをテストさせて，それがどちらの質感記憶課題

265–272.

の活動であるかを判別させた。

[3] Ferwerda, J., A., et al. (2001). Proceedings of

4 結果・考察

SPIE, 4299, 291-301.
[4] Fleming, R., W., et al. (2003). J. Vision, 3,
347–368.

腹側高次視覚野と初期視覚野，頭頂間溝における MVPA の正答率の推移を図に示す。腹側高次
視覚野において，課題開始 6 秒前後の応答からチャンスレベル以上で課題の判別に成功した。課題

58 − Academic Center for Computing and Media Studies, Kyoto University[5] Motoyoshi, I., et al. (2007). Nature, 447, 206–
209.
[6]  Katherine, C., B., & Yaoda, X. (2016).
Nature Neurosci, 19, 150-157.

臨界レイノルズ数における球の流体力に対する周囲の乱流場の影響
武藤昌也
京都大学大学院工学研究科機械理工学専攻

臨界レイノルズ数条件下における一様気流中の球に付加する流体力が，周囲気流の乱流場の乱れの強さによってどのような影響を受けるかを調べた．計算の結果，乱れのない一様流中の球と比較して，球表面の圧力係数の変動や球後流の変動流が抑制されることが分かった．

1 緒言
乱流中での粒子や気泡のような分散相の運動を理解するには，それに掛かる様々な力について考慮することが重要である．しかし，従来，粒子の抗力係数は乱れの小さい一様流中で計測された実験値に基づいている．実際の工業プロセス内の分散固気二相流では，微粒子が分散する流れ場が乱れの小さい一様流であることはほとんど無い．こ図 1 球周り流れの解析領域（Muto et al.1)2)3)）の場合，境界層の発達やはく離の位置はもはや一様流中とは同様ではなく，非対称，非一様な流体

ズム要素，その他は四面体要素で空間的に分割さ力の作用により，その運動が一様流中とは大きくれている．解析手法としては，非構造有限体積法

異なることが容易に推測される．昨今，工業プロに基づく Large-eddy simulation を用い，解析領域のセスにおける燃焼，
粉体搬送技術の高度化に伴い，流入境界からは，一様流に乱れを重ね合わせた流シミュレーションには解析精度の向上が要求され速を流入させた．乱れの生成手法として，Reynoldsているが，多くの場合，微粒子の運動に対する周応力のコレスキー分解に基づく方法 4）を用いた．

囲乱流場の影響が考慮されることはほぼ無い．

解析条件としては，平均流に基づくレイノルズ数

そこで本研究は，臨界レイノルズ数条件下におが臨界レイノルズ数に近い 2.0×105 になるようにいて，一様気流中の球に付加する流体力が，周囲主流流速を与えた．数値解析コードには，非構造

気流の乱流場の乱れの強さによって受ける影響を格子有限体積法による非定常乱流解析コード

詳細に調べることを目的とする．

FrontFlowRed/Comb を用いた．

2 計算手法

3 計算結果と考察

本解析では，臨界レイノルズ数条件下での球周りの流れ場解析

1)2)3)

図 2，3 に，球周りに一様流のみを与えた場合，に用いられた，解析領域およおよび一様流に乱れを重ね合わせた流速を与えた

び解析手法を用いる．解析領域の概略図を図 1 に場合の一様流方向無次元速度の瞬時分布および球示す．三次元矩形流路の中心に固定された直径
44.7 mm の滑球周りの流れ場が解析される．解析
格子は，球表面近傍および壁面近傍においてプリ

表面の圧力係数 Cp の分布を示す．図 2 に対して図3 では，球の上流澱み点から 90 度以降の下流表面において，Cp の変動や後流流れ場の流速変動が抑京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 59

図 2 球周りの瞬時流速場（変動風なし，Muto et al.2)）

図 4 球表面の圧力係数分布（乱れ無しのデータは武藤ら 3）による）

4 結言
本研究では，臨界レイノルズ数条件下における一様気流中の球に付加する流体力が，周囲気流の乱流場の乱れの強さによってどのような影響を受けるかを調べた．計算の結果，乱れのない一様流中の球と比較して，球表面の圧力係数の変動や球
図 3 球周りの瞬時流速場（変動風あり）

後流の変動流が抑制されることが分かった．
今後，
流体抗力を含む詳細な解析の実施が望まれる．

制されていることが分かる．これは，一様流に乱れを付加することにより，球表面の乱流遷移が促

参考文献

進され，一様流のみを与えた場合に見られた空間
的に大きな流速変動がむしろ抑制されたことによると推測される．
次に，図 4 に図上流澱み点を 0 度とした y = 0
断面内の圧力係数分布をプロットする．ただし，一様流のみを与えた場合は時間平均分布であり，一様流に乱れを重ね合わせた流速を与えた場合は瞬時分布である．一様流のみを与えた場合には，上流澱み点から 92 度における層流剥離に伴う圧力上昇と 98 度における再付着に伴う圧力低下の

1) M. Muto, M. Tsubokura, N. Oshima, Numericalvisualization of boundary layer transition whennegative Magnus effect occurs, Journal of Visualization,
Vol. 15(3), 261-268, 20122) M. Muto, M. Tsubokura, N. Oshima, Negative

Magnus lift on a rotating sphere at around the critical

Reynolds number, Physics of Fluids, Vol. 24, 014102
20123) 武藤昌也，坪倉誠，大島まり，大島伸行，周囲に特徴が見られる．しかしながら，一様流に乱れを周期変動が与えられた球に負荷する流体力，日本

重ね合わせた流速を与えた場合は，一様流のみを機械学会論文集 B 編，Vol. 76，563-569，2010

与えた場合と比較して，90 度近傍での圧力低下が大きく，より層流状態に近づくことが分かる．このことも，球周り流れの乱流遷移が促進され，境界層の剥離および再付着現象が下流に移動することを示唆している．今後，流体抗力を含むより詳細な解析の実施が望まれる．

60 − Academic Center for Computing and Media Studies, Kyoto University

4) Z.T. Xie and I.P. Castro, Efficient Generation of

Inflow Conditions for Large Eddy Simulation ofStreet-Scale Flows, Flow, Turbulence and Combustion,

Vol. 81(3), 449-470, 2008

3 次元ナビエストークス流の過渡的乱れ
佐々木 英一同志社大学理工学部

系のエネルギーを一定に保つ線形外力によって駆動される 3 次元周期箱領域の Navier-Stokes
流れの過渡的乱れについて調査した．この系における自明解である 1 次元せん断流が任意のReynolds 数で線形安定であることを示した．
乱流状態を初期条件に数値シミュレーションを行うと，
動粘性係数が増加すると乱れは過渡的になり十分時間がたつと定常 Beltrami 流に漸近する．

1 はじめに
3 次元周期箱流れは乱流の普遍的な統計則が発現する等方乱流の典型的なモデルである．近年，エネルギー注入率を一定に保った外力の下での 3次元周期箱流れについて，乱流状態が過渡的となり，最終的に速度ベクトルと渦度ベクトルが平行

𝒇𝒇̂(𝒌𝒌, 𝑡𝑡) ={𝜖𝜖(𝑡𝑡)
̂(𝒌𝒌, 𝑡𝑡)𝒖𝒖
2𝐸𝐸𝑓𝑓 (𝑡𝑡)0

for 0 < |𝒌𝒌| < 𝑘𝑘𝑓𝑓otherwisê(𝒌𝒌, 𝑡𝑡)はそれぞれ外力ととした．ここで𝒇𝒇̂(𝒌𝒌, 𝑡𝑡)と𝒖𝒖
速度場の波数𝒌𝒌の Fourier 係数であり，𝐸𝐸𝑓𝑓 は波数ベクトルの大きさ𝑘𝑘𝑓𝑓 = 2.5で Low-path フィルタな流れである定常 Beltrami 流れとなることが報ーした速度場のエネルギー，𝜖𝜖はエネルギー散逸告された(Linkmann & Morozov 2015, McComb率である．この外力により系の運動エネルギーはet al. 2015)．壁面を有する流れでも過渡的乱れが一定に保たれる．

存在することが報告されており(例えば Avila et

数値計算は古典的ルンゲクッタ法と Fourier 展

al. 2011)， 壁の有無にかかわらず過渡的乱れの発開によるスペクトル法を用いた．エイリアジング

現は遷移過程における共通した性質であることが除去のため，一辺の格子点数𝑁𝑁と切断波数𝑇𝑇は𝑁𝑁 >示唆される．
我々は系のエネルギーを一定に保つ線形外力によって駆動される 3 次元周期箱 Navier-Stokes 流

3𝑇𝑇 + 1を満たす．格子点数は𝑁𝑁 = 48とし，時間差分間隔は10−3 とした．

3 結果と今後の課題

れの過渡的乱れについて調査した．単純なせん断流が任意の Reynolds 数で線形安定であることを一つの方向にだけ速度を持つ１次元せん断流が証明した．乱流状態を初期条件に数値シミュレー

任意の Reynolds 数で線形安定であることを証明

ションを行うと，動粘性係数が増加すると過渡的した．このことは，乱流状態と層流状態を相空間

乱れとなり定常な Beltrami 流に近づくことを見で分割するエッジ状態が存在することを意味する．

つけた．

乱流状態を初期条件に数値シミュレーションを行

2 問題設定うと，動粘性係数を増加させると乱れが過渡的になり，定常な Beltrami 流に近づくことが分かっ我々が扱う方程式は３次元周期境界条件におけた．このことは 1 次元せん断流以外にも安定解がる非圧縮 Navier-Stokes 方程式である：

存在することを意味するが，安定解の吸引領域が∂𝑡𝑡 𝒖𝒖 + 𝒖𝒖 ∙ 𝛁𝛁𝒖𝒖 = −𝛁𝛁𝑝𝑝 + 𝜈𝜈∆𝒖𝒖 + 𝒇𝒇,

相空間をどのように分割しているかわかっていな𝛁𝛁 ∙ 𝒖𝒖 = 0.

ここで，𝒖𝒖は速度，𝒇𝒇は外力，𝜈𝜈は動粘性係数，𝑝𝑝はい．乱流状態と安定解のエッジ状態を数値的に見圧力である．周期箱寸法は一辺の長さ𝐿𝐿 = 2𝜋𝜋の立つけ，不変解によって乱流状態からの離脱の機構

方体領域である．外力はを特徴づけることが今後の課題である．

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 61

プログラム高度化支援事業研究報告

動的／静的水～土骨格連成有限変形解析コードの高度化～固有値解析に関わる高速化、並列化の検討～
野田利弘

豊田智大

名古屋大学 減災連携研究センター

名古屋大学大学院 工学研究科

1 平成 28 年度事業の概要と目的

2 Foss の方法による固有値解析の概要

筆者らが開発した地盤の動的／静的変形解析コー

当該解析コードでは、材料非線形性を考慮するたド（GEOASIA® [1], [2]）に対して、平成 24 年度（下めに増分型の弾塑性体構成式を用い、幾何学的非線半期）[3], 25 年度[4], 26 年度[5]に続き、平成 28 年

形性を考慮するために updated Lagrange 法を用

度も京都大学学術情報メディアセンターのプログラいる。このような方法で動的境界値問題を扱うため

ム高度化支援事業に採択していただいた。この概要には，速度型運動方程式を解いて運動を求める必要を報告する。

がある。そこで、速度型運動方程式の弱形式を有限本解析コードは、粘土や砂そして細粒分を多く含

要素離散化した上で、土骨格の構成式として有効応

む中間土までを一貫して記述する弾塑性構成式（All

力を用いて記述した弾塑性構成式を適用する。
また、

Soils）を搭載し、地盤の変形から破壊までの挙動

水～土骨格連成式については、飽和土の連続式およ（All States）を、動的と静的とを問わずあらゆる

び平均的な間隙水の流速式に基づいて、各要素の中外力条件（All Round）の下で解析することを目指

心に間隙水圧 u を代表させる田村流[8]の物理モデ

した解析コードである。
本事業における平成 24～26

ルを拡張して用いる。さらに、有限要素法の速度場

年度の目的は、大規模・高速化計算が実現できるよに線形制約条件が課せられる場合を考慮すると、最う当該解析コードの高度化・高性能化を図ることで終的に解くべき連立常微分方程式は式(1)のようにあり、OpenMP と MPI とのハイブリッド並列化、表される[1], [2]。

反復法ソルバである PETSc の組み込みや、領域分割法による並列処理による計算手法の導入などを実施した。
平成 28 年度の目的は、固有値解析に関する部分

 



 Mv C* v  Kv LT u CT    f
 Lv Lv  Hu Gu  f u  Cv  0


 

(1)

ここに、M は質量マトリクス、K は接線剛性マトの高度化で、具体的には、一般化固有値解析に対応

リクス、L は土骨格の変位速度を土骨格の体積変化

した数値演算ライブラリの調査に加え、当初用いてに変換するマトリクス、v は節点の変位速度ベクト

いた Intel MKL（LAPACK）の GGEV[6]をオープ

ル、u は要素中心における土骨格から見た間隙水圧

ンソース SLEPc（Scalable Library for Eigenvaue
Problem Computations）[7]に変更し、性能・計算結果の比較を実施した。

速度、L は L を変形して作られるマトリクス、H

は透水係数マトリクス、 G は間隙率と水の圧縮率で与えられるマトリクスである。また、 C*  は粘性

本稿では、当該解析手法における固有値解析の概境界を導入したことによって生じる減衰マトリクス

要を説明した後、高度化を行った際の動作検証と性であり、非比例減衰である。C は節点の運動に課さ

能評価について紹介する。

れる線形制約条件を表すマトリクス、   はLagrange の未定乗数であり、 CT   は制約条件を実現するのに必要な節点に作用する束縛力としての62 − Academic Center for Computing and Media Studies, Kyoto University

式(1)で示す M と K が疎行列で値が非対称である

物理的意味を持つ。
ここで、Foss の方法[9]に倣って、 v  w とおくことを考慮して、行列形式を圧縮行格納方式（CSR）と、式(1)の斉次方程式は、式(2)のように書きかえるに変更した上で、並列化とあわせて非対称行列に対ことができる[10]。

応可能なオープンソースSLEPc に変更した。
なお、
この固有値ライブラリは、平成 26 年度までの高度  M w  C* v K v L T u  C T    0 M w M v  0 Lv  L v H u  G u  0  Cv  0
(2)

示す 2 次元有限要素メッシュである次の 2 種類を使  M  L K 
OO   M B   
 L  O  C O OO OG  OO O
T

O
O
O

用した。データ①：行列規模 1600×1600（非零設
定箇所：24856（全要素の 0.97%））, データ②：行
列規模 3660×3660（非零設定箇所：101152（全要 v
O  CT O O  , x  w 
H  O  u
 O O 

素の 0.76%））である。なお、データ①は一相系の(3),

x0   v0 T w0 T u0 T  0 T  etT

場合、データ②は二相系の場合について固有値解析を実施している。また、データ②では、フルマトリ
クスでメモリを約 102MB／行列を使用している。

ここで、式(3)の解 x0  が、x  x0 et

ラリ群である。
今回の高度化による動作検証を行うため、図 1 にで次式のように表される。

Ax  Bx ,Toolkit for Scientific Computation）[11]を利用して開発されている固有値解析に関する数値演算ライブ

式(2)の連立常微分方程式はさらに、マトリクス表現

 C*A   M L  O

化の際に搭載した PETSc（Portable, Extensible(4)(a) データ①：1600×1600（左端節点を拘束）と表されると仮定すると、次式で表される一般固有値問題が得られる。
 Ax 0   Bx 0 (5)複素固有値  の虚部によって減衰を考慮した固有振動数が，複素固有ベクトル x0  の一部分である

v0  の実部と虚部の線形結合によって固有振動モードが与えられる。なお、この Foss の方法を用いる(b) データ②：3660×3660 （下部節点を拘束）と、解くべき元の数は、変形解析時に比べて約２倍

図１ 動作確認に用いたメッシュ図になる。

3 高度化（SLEPc へ変更）の概要および動作検証と性能評価a)性能計測結果
性能計測結果として、Intel MKL（LAPACK）のGGEV との比較を表１に示す。なお、性能計測にあたっては、SLEPc のデフォルトの Krylovschur 法

前章で述べた解析手法を用いたパソコンによる解や Jacobi-Davidson 法など固有値ソルバによる違い

析事例を文献[10]に示すが、この事例では、フルサや求解固有値数などの影響を調べている。
ここでは、

イズ（密）行列を入力として Intel MKL の GGEV

シリアル実行時で最も効果が得られた場合について、を適用しているが、このライブラリには(1)密行列を固有値求解解析に要した部分のみ時間を示している。

対象としているため、規模が大きくなるに従い、行列の値の保管に要する配列のメモリ量が非常に大き

表１ 高速化（速度向上比）の効果

くなる、(2)並列化に対応していない、などの問題が

データ

GGEV

PETSc

速度向上比

あるため、特にパソコンレベルでは十分な計算がで①

67.43 s1.67 s40.4

きなかった。

②

812.55 s

1.70 s480.0

そこで、一般固有値解析の高度化を図るのに際し、京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 63

データ規模が大きくなるに従い、GGEV に比べ速度値の計算結果の際は生じていないため、固有値解析

向上が見られる。なお、SLEPc において並列化適用の入力に使われる材料の判定条件（弾性 or 弾塑性）（プロセス並列、スレッド並列）時も調べている。

がコンパイラ最適化オプションなどの設定条件に敏

このデータ①に、②については、シリアル実行時と感に影響していた。

比べて計算結果に差異はほとんど見られなかったが、データ②において行と列をそれぞれ 10 倍にしたデータなどデータ規模や条件設定によっては並列化による速度向上が得られる場合もあった。
b)計算結果の検証
固有値、固有ベクトル（固有モード）などについて行った計算結果の比較検証結果について述べる。
(a) 1 次モード

図２と図３はそれぞれ、データ①と②について、固有値が最も小さい固有モードから順に 3 つ目までを示す。
データ①については設定条件による計算結果の差
異はほぼ見られなかった。一方で、データ②については、図３に示すように、設定条件により固有値や固有モードにやや差異が生じた。なお、固有値ライ
ブラリの置き換え（GGEV or SLEPc）に伴い固有(b) ２次モード(c) ３次モード
図３ 計算結果（固有モードの比較）（データ②）(a) 1 次モード

4 おわりに
平成 28 年度では、当該解析コードの一般固有値
解析に関する部分において、高度化（SLEPc への変更）を図り、シンプルな形状に対する固有値解析を実施した。

(b) ２次モード

・SLEPc 適用の際、行列の格納形式を密行列から疎行列（CSR 形式）へ変更し、これにより固有値解析の行列（剛性マトリクス、質量マトリクス）に対応する配列のメモリ量を大幅に減らすことができた。

(c) ３次モード

・SLEPc での各種条件設定を行うことで、固有値求

図２ 計算結果（固有モード）（データ①）

解箇所の性能比較では、GGEV と比較して、データ

64 − Academic Center for Computing and Media Studies, Kyoto University①（節点数 805、要素数 640）で約 40 倍、データ②

Evaluation, and Application, Hong Kong, China, 11-27.

（節点数 1891、要素数 1800）で約 480 倍の速度向[2] Noda, T., Asaoka, A. and Nakano, M. (2008):

上できた。

Soil-water coupled finite deformation analysis based on a・並列化（プロセス並列、スレッド並列）では、シ

rate-type equation of motion incorporating the SYS

リアル実行の最速ケースでの更なる性能向上は今回

Cam-clay model, Soils and Foundations, 48(.6),の使用データでは見られなかった。一方で、データ

771-790.

規模を大きくしたケースでは、使用するデータや設[3] 野田利弘 (2013): 京都大学学術情報メディアセ

定条件によって、並列化による速度向上が見込めるンター全国共同研究利用版広報，Vol.12, No.2, 33-35.

ことを確認できた。

[4] 野田利弘 (2014):京都大学学術情報メディアセ

なお、現時点では性能計測までは実施していない

ンター全国共同研究利用版広報，Vol.13, No.2, 44-47.

が、パソコンレベルでは時間的にも容量的にも計算[5] 野田利弘 (2015): 京都大学学術情報メディアセできなかった３次元有限要素メッシュ（節点数

ンター全国共同研究利用版広報，Vol.14, No.2, 31-34.

13971、要素数 12000）に対して今回の解析コード[6] Intel MKL（LAPACK）, GGE: https://software.intel.

を適用した事例を図４に示しておく。ここでは直方

com/en-us/mkl-developer-reference-c-ggev.

体と形状は単純であるが、今後は複雑な３次元要素[7] SLEPc: http://slepc.upv.es/.

についても計算を実施してゆきたい。

[8]赤井浩一, 田村武(1978): 弾塑性構成式による多次元圧密の数値解析, 土木学会論文集, 第 269 号,

謝辞 京都大学学術情報メディアセンターの中島

95-104.

浩教授および牛島省教授ならびにクレイ・ジャパ[9] Foss, K, A. (1958): Coordinates which uncouple theン・インクの方々には、上述以外の点においてもごequations of motion of damped linear dynamic systems,

尽力いただきました。ここに深く謝意を表します。

Journal of Applied Mechanics, ASME, 32(3), 361-364.
[10] Shimizu, R. and Yamada S. (2013): Effect of seismic
waves with different dominant frequencies on the delayed

参考文献

failure behavior of a soil structure-ground system, Proc. of[1] Asaoka, A. and Noda, T. (2007): All soils all states all

the 18th International Conference on Soil Mechanics andround geo-analysis integration, International Workshop on

Geotechnical Engineering, Paris 2013, 1599-1602.

Constitutive Modelling - Development, Implementation,[11] PETSc/Tao: https://www.mcs.anl.gov/petsc/.

(a) １次モード(b) ２次モード(c) ３次モード

図４ ３次元の固有値解析の事例（直方体、二相系、行列規模：51690×51690、底面節点を拘束）

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 65

津波による瓦礫浮遊物を考慮した構造流体連成解析
高瀬
∗

1

慎介 ∗

八戸工業大学 工学部 土木建築工学科

はじめにし，この数学被覆と物理領域の交わりを物理被覆として定義する．物理被覆には，物理変数を近似する

災害大国である我が国において，防災・減災を目ための被覆関数と呼ばれる未知パラメータを含む関的とした対策事業は幅広く展開されているが，それ数を導入する．被覆関数には，PU 重み関数との積にでもなお毎年のように災害の被害は繰り返し起こっより構成される基底関数が一次独立性を満たす限りている．これらの災害被害を効率的に制御するには，任意の関数形を利用できるが，本研究では，定数項想定する災害に対して信頼性の高い被害予測を行う

のみを未知パラメータとして用いることにする．

ことが極めて重要な課題である．

そこで本研究では，これらの自然災害の中でも津波漂流物および流体力による構造物の破壊に着目し，流体力および漂流物の衝突による構造物の変形・破壊過程が引き起こす災害の解析手法の開発を行う．解析手法として，
流体力の計算には安定化有限被覆法，構造物の移動・衝突には個別要素法を用いて，相互の連成解析を行う．昨年度の本プログラム高度化により MPI による並列が未対応であった個別要素法の解析部分が並列化され，MPI を用いた並列計算が可能となった．し

本研究では，FCM の物理要素の基本形状としては

1 次の四面体要素を用いることにする．そして，被覆関数は定数項のみとするので，本研究における FCMの物理要素は，必然的に FEM における 1 次の四面体

有限要素と等価なものとなる．ただし，FCM では流
体の数学領域内に構造物が存在してもよく，流体の
物理要素の境界は数学要素と独立に動くことができる．この特徴を利用して，構造物の境界位置を正確に考慮したうえで流れと構造物の剛体運動の相互作用を評価可能な数値解析手法を構築する．なお，流かし，個別要素法では，逐次領域内の演算量が変化

体と構造物の境界位置は，構造物表面をゼロとする

するため，広範囲での浮遊する構造物の解析では，そ

距離関数である Levelset 関数作成し，流体領域の定の部分がボトルネックになっている．本手法の開発義を行っている．

により，複数の構造物をモデル化し，実地形等を考慮した，実際の条件に近い津波被害予測解析が可能となり，工学的にも有効な被害予測手法になる．

2.2

流体に関する定式化

本研究で用いる支配方程式は，非圧縮粘性流体に

2

流体の数値解析手法

2.1

有限被覆法

有限被覆法 [1]（FCM) は，近似関数が数学領域 Ω Mと支配方程式が満たされるべき物理領域 ΩP を独立しおける Navier-Stokes の運動方程式と連続式で以下のように表される．
)( ∂u+ u · ∇u − f − ∇ · σ (u, p) = 0
ρ
∂t

∇·u=0

(1)
(2)

て定義する．そして，FCM では有限要素法（FEM）

ここで，ρ は密度，u は流速ベクトル， p は圧力，fにおける形状関数の代わりに，解析対象となる物理

は物体力ベクトル，σ は応力テンソルである．また，領域とは独立な PU 条件（Partition of Unity) を満た

Newton 流体を仮定し，構成則には次式を用いる．

す近似関数（以下，PU 重み関数）を数学被覆に導入

66 − Academic Center for Computing and Media Studies, Kyoto University

σ = −pI + 2µε(u)(3)

ここで，I は 2 階の単位テンソル，µ は粘性係数であり，ε(u) は次式で定義される変形速度テンソルで
ある．

2.3

)
1(
∇u + (∇u)Tε(u) =
2(4)的手法の１つである Phase-ﬁeld 法を採用することにする．Phase-ﬁeld 法では，次式に示す保存形式に修正された Allen-Cahn 型移流方程式 [5] を解くことで自由表面位置を決定する．

∂ϕ∂t

安定化有限被覆法

3 次元流れ場の支配方程式 (1) に対して有限被覆法を用いた SUPG/PSPG 法 [2],[3] を適用すると，次式のような安定化有限被覆法による離散化方程式が得
られる．
( h
)
∫
∂u
h
h
h
+ u · ∇u − f dΩ
ρ
w ·ρ
P
∂t
∫Ω
∫
h
h h
+
ε(w ) : σ(u , p ) dΩ +qh ∇ · uh dΩ
ΩP

e=1

∫

h

また，ϕ は Phase-ﬁeld 変数を表し，気体であれば 0.0，

液体であれば 1.0，自由表面上であれば 0.5 をとるものとする．そして，各要素における流体の密度 ρ と粘性係数 µ は，液体（水）と気体（空気）の密度 ρl , ρgと粘性係数 µl , µg ，Phase-ﬁeld 関数 ϕ を用いて次式のように求められる．

ρ = ρl ϕ + ρg (1 − ϕ)hΩePp̄(uh − ūh ) · wh = 0

(5)

(6)

ここで，ϵ は移動度，δ は界面領域の代表幅である．

ΩP}
1
τsupg u · ∇w + τpspg ∇q+Pρ
e=1 Ωe
{ ( h
)
}
∂u
h
h
h h
+ u · ∇u − f − ∇ · σ(u , p ) dΩ
·
ρ
∂t
nel ∫
∑
+
τcont ∇ · wh ρ∇ · uh dΩ
{

nel ∫
∑

ϵ
+ u · ∇ϕ = ∇ · (δ(∇ϕ) − Fa ) ,δ∇ϕ
Fa = ϕ(1 − ϕ)|∇ϕ|µ = µl ϕ + µg (1 − ϕ)

(7)
(8)

式 (6) に対して，SUPG 法に基づく安定化有限要素法を適用し，解析を行っている．詳しい定式化については，文献 [6] を参照されたい．

3

構造の数値解析手法

ここで，ΩP は Navier-Stokes 方程式の物理領域，uh ，

3.1

個別要素法を用いた剛体要素の近似式，p̄ はペナルティー定数，ūh はペナルティー

結することで複雑な剛体モデルを表現することが可+Γgph は，それぞれ速度と圧力の有限要素近似式，wh ，qhは，それぞれ運動方程式と連続式に対する重み関数

境界面での速度である．安定化項である第 4,5 項は要素ごとに定義される不連続量であるため，要素ご

球要素を用いる個別要素法 [7] では，球要素を剛能である．このとき，剛体を構成する球要素がそれぞれ個別に接触判定を行っており，球要素ごとに計との領域 ΩeP での積分の総和となる．ここで，e は算された接触点と接触力の情報を剛体の重心に作用流の卓越に対して安定化を施す SUPG 項，および圧

すなわち，剛体要素 g の重心に関する次式の運動方

要素番号，nel は要素数である．また，第 4 項は移する力に換算することで剛体の運動を表現している．

力振動を回避するための PSPG 項であり，第５項は程式を解く．

自由表面の数値不安定性を回避するための衝撃捕捉（Shock-Capturing) 項 [4]，ペナルティ項は構造壁面
上での Dirichet 型境界条件を考慮するために付加したものである．また，τsupg ，τpspg ，τcont は，すべて安定化パラメータである．

d(mg vg )
= Fg
(9)
dt
d(Ig ωg )
= Mg
(10)
dt
ここで，mg は剛体の質量，vg は剛体の速度ベクトル，Ig は慣性モーメントテンソル，ωg は剛体の角速度ベクトルである．重心に作用する力 Fg 及びモーメ

ント Mg は剛体を構成する球要素 i に対して接触し2.4 Phase-ﬁeld 法を用いた界面捕捉ているすべての要素に関する接触力の総和として以本研究では，砕波等の複雑な自由表面形状を表現する必要性があるため，固定メッシュを用いた Euler

下のように書き下せる．

Fg =

i
∑

(Fn + F s )(11)京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 67

図 1: Cohesive model による DEM 剛体モデルの連結

Mg =

i
∑

(N × F s )(12)

このとき，xi を球要素 i の重心座標，xg を剛体の重心座標として N = (xi + ri n) − xg である．なお，ri は球要素 i の半径，n は球要素 i の重心から接触している球要素の重心へと向かう法線方向の単位ベクトルである．また，本研究では, 球要素は剛体モデルの表面にのみ配置し，計算負荷を軽減している．剛体の姿勢管理には，クォータニオン（四元数）を導入している．
式 (9) から，時間ステップ n と時間ステップ n + 1

における剛体の速度ベクトルを vng ，vn+1
とすると，
g

vn+1
= vng + (Fg /mg )∆tg(13)となる．ここで，∆t は微小時間増分量である．また，
式 (10) から，時間ステップ n と時間ステップ n + 1 にとすると，おける剛体の角速度ベクトルを ωng ，ωn+1
g
−1 n
n
ωn+1
= ωng + (I−1gg Mg − Ig ωg × Ig ωg )∆t

図 2: Cohesive model のイメージ(14)となる．また，剛体重心の位置座標 xg の更新は，剛( )[( ) ( )2 ]


27
δ

 4 σc δ 1 − 2 δδ + δδ
,

f
f
f
σcoh = 



 0,

δ ≤ δf ,(16)

δ > δf ,ここで，σcoh は結合力，σc は結合力の上限値，δ は開口量，δ f は完全な破壊に至る開口量を意味する．

このモデルのイメージ図を 図-2 に示す．なお，本研

究では，負荷時による破壊に着目しているおり，簡便のため，除荷時も式 (16) を用いるモデル化を行っ
ている．

4

数値解析例
数値解析例として，複数ブロックの構造流体連成

解析を示す．昨年のプログラム高度化において，非並列であった DEM の計算に対して，DEM の計算で用いるバックセルを領域分割することで MPI による並列計算が可能となった．図-3(a) に示す赤いボック
スが各プロセッサに割り当てられる領域を表してい体の速度ベクトルを用いて次式により更新する．

る．この方法では，図-3(b) に示すように計算が進む

xn+1= xng + vng ∆tgにつれ，剛体要素が偏ることにより，各プロセッサ(15)のロードバランスが悪くなり，並列度が低下していた．今年度のプログラム高度化では，図-3(c)(d) に示
すように，剛体要素を包括するようにバックセルを

3.2 Cohesive model による破壊の表現
本研究では，図-1 に示すように，剛体モデルの表面

動的に配置し，領域分割することで，並列密度を高くした．そのため，剛体要素の集積状態に合わせて計算領域を再設定しているので，並列度の低下も防を潜在的な不連続面として，この面を Cohesive model

げるようになっている．この方法により，図-4 に示を用いて連結することで破壊挙動を表現する．本研

すように，並列化性能の向上が見られた．

究で用いた Cohesive model は，多項式を用いたものであり，次式で表現される．

68 − Academic Center for Computing and Media Studies, Kyoto University図 4: 並列化性能

参考文献[1] Terada, K., Asai, M. and Yamagishi, M.: Finite Covermethod for linear and nonlinear analyses of heterogeneous solids, International Journal for Numerical Methods in Engineering, Vol.58, Issue 9, pp.1321-1346,
2003.
[2] Brooks, A.N., Hughes, T.J.R.:
streamlineupwind/Petrov-Galerkin formulations for convection
dominated ﬂows with particular emphasis on theincompressible Navier-Stokes equations, Computer
Methods in Applied Mechanics and Engineering, 32,
pp.199-259, 1982.
[3] Tezduyar, T.E.: Stabilized ﬁnite element formulations
for incompressible ﬂow computations : Advanced in
Applied Mechanics, 28, pp.1-44, 1991.
[4] Aliabadi,S., and Tezduyar, T.E.: Stabilized-ﬁniteelement/interface-caputuring technique for parallel
computation of unsteady ﬂows with interfaces, Compute Methods in Applied Mechanics and Engineering,
190, pp.243-261, 2000.

図 3: バックセルと解析結果

5[5] N. Tanaka, J. Matsumoto and S. Matsumoto: Phaseﬁeld model-based simulation of motions of a two-phaseﬂuid on solid surface, Journal of Computational Scienceand Technology, Vol.7 No.2, pp.322-337, 2013.

おわりに本研究では，構造物の破壊過程も考慮可能な構造

流体連成解析手法を構築し，その並列化性能向上のため，DEM を用いた剛体要素の並列化について検討を行った．今年度導入した並列化方法により，昨年度に比べ，約 1.5 倍の並列化性能の向上が見られた．
これにより，より大規模な実問題への適用が可能となる．

[6] 高瀬慎介，森口周二，寺田賢二郎，小山直輝，金子賢治，車谷麻緒，加藤準治，京谷孝史：構造物の破壊を考慮した有限被覆法に基づく構造流体連成解析，
土木学会論文集 A2(応用力学），Vol.71,No.2（応用力学論文集 Vol.18), I 203-I 212, 2015.
[7] Cundall, P. A. : A Computer model for simulating progressive, Large Scale movement in blocky rock system,
Proceedings of ISRM Symposium, pp.11-18, 1971.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 69

サービスの記録・報告システム A 運転状況 ( 2016 年 10 月 ～ 2017 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時

サービス再開日時

システムダウン障害発生状況
保守時間[h]障害発生日時

サービス再開日時

ダウン時間[h]2016/10/01

0:00

2016/10/03

9:00

57.00

2016/11/22

15:12

2016/11/22

18:46

3.57

2016/10/10

6:30

2016/10/12

9:45

51.25

2017/01/21

0:25

2017/01/22

21:15

44.83

2016/11/07

9:00

2016/11/08

9:30

24.50

2017/01/23

15:30

2017/01/23

16:28

0.97

2016/12/06

9:00

2016/12/06

18:00

9.00

2017/02/06

8:53

2017/02/06

10:00

1.12

2016/12/19

9:00

2016/12/27

10:00

193.002017/02/07

12:30

2017/02/08

12:00

23.50

2017/03/27

9:00

2017/04/01

0:00

111:002017/03/07

20:32

2017/03/08

9:40

13.13

2017/03/15

2:36

2017/03/15

15:00

12.40

２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

692.75

16,48530,992

12,416,000

8,738,750

1499.5

21

%

11 月

692.5015,19258,473

21,149,300

14,375,800

1577.2

30

%

12 月

542.0028,102

68,27917,676,300

12,430,700

1630.437%

1月

698.2043,034

183,775

40,817,900

24,715,500

1770.453%

2月

647.3837,557

105,909

43,118,300

29,644,300

1783.662%

3月

605.83

89,115

155,250

44,271,800

27,890,000

1733.055%

計3878.66

229,485

602,677

179,449,600

117,795,050

1,665.7

43

%

1月

2月

ノード
利用率

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

70 − Academic Center for Computing and Media Studies, Kyoto University

システム B 運転状況 ( 2016 年 10 月 ～ 2017 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時

サービス再開日時

システムダウン障害発生状況
保守時間[h]障害発生日時

サービス再開日時

ダウン時間[h]2016/10/01

0:00

2016/10/03

9:00

57.00

2016/11/22

15:12

2016/11/22

18:46

3.57

2016/10/10

6:30

2016/10/12

9:45

51.25

2017/01/08

14:00

2017/01/08

16:25

2.42

2016/11/07

9:00

2016/11/08

9:30

24.50

2017/01/21

0:25

2017/01/22

21:15

44.83

2016/12/19

9:00

2016/12/27

10:00

193.002017/01/23

15:30

2017/01/23

16:28

0.97

2017/03/27

9:00

2017/04/01

0:00

111:002017/01/31

20:40

2017/02/01

5:30

8.83

2017/02/20

21:00

2017/02/21

9:15

12.25

2017/03/13

0:00

2017/03/13

22:00

22.00

2017/03/15

2:36

2017/03/15

15:00

12.40

２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

635.7536,467280,1132,514,810

2,014,580

585.8

49

%

11 月

695.5042,320491,993

3,841,410

3,295,640

588.8

58

%

12 月

551.0061,283

443,578

3,314,070

2,725,630

646.0

54

%

1月

692.4589,628481,289

10,063,000

8,072,530

809.2

45

%

2月

654.25

174,571584,524

9,578,440

8,187,060

839.3

46

%

3月

596.00

194,598

473,067

9,177,120

7,862,150

807.1

35

%

計

3824.95

598,867

2,754,564

38,488,850

32,157,590

712.7

48

%

ノード
利用率

※ 10 ‐ 12 月は旧システム、1‐3 月は新システムの稼働状況

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 71

システム C 運転状況 ( 2016 年 10 月 ～ 2017 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]障害発生日時

サービス再開日時

ダウン時間[h]2016/10/01

0:00

2016/10/03

9:00

57.00

2016/11/22

15:12

2016/11/22

18:46

3.57

2016/10/10

6:30

2016/10/12

9:45

51.25

2017/01/21

0:25

2017/01/22

21:15

44.83

2016/11/07

9:00

2016/11/08

9:30

24.50

2017/01/23

15:30

2017/01/23

16:28

0.97

2016/12/19

9:00

2016/12/27

10:00

193.002017/01/31

20:40

2017/02/01

5:30

8.83

2017/03/27

9:00

2017/04/01

0:00

111:002017/03/15

2:36

2017/03/15

15:00

12.40

２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

635.75300

3,201

22,97919,848

12.2

22

%

11 月

695.50612

11,13864,14659,756

12.8

34

%

12 月

551.00

1,406

5,152

31,921

25,625

12.7

22

%

1月

694.87

1,623

12,882

88,60447,019

14.8

15

%

2月

666.50925

58,977

121,068

151,266

16.0

9

%

3月

618.00601

3,471

44,92830,486

15.3

5

%

計

3861.62

5,467

94,821

373,647

333,999

14.0

18

%

ノード
利用率

※ 10 ‐ 12 月は旧システム、1‐3 月は新システムの稼働状況

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

72 − Academic Center for Computing and Media Studies, Kyoto University

システム D 運転状況 ( 2016 年 10 月 ～ 2017 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]2016/10/01

0:00

2016/10/03

9:00

57.00

2016/10/10

6:30

2016/10/12

9:45

51.25

2016/11/07

9:00

2016/11/08

9:30

24.50

障害発生日時

サービス再開日時‐なし

ダウン時間[h]‐２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

635.755,744

35,4862,910,210

191,477

415.6

39

%

11 月

695.505,959

45,9294,534,940

117,787

415.9

55

%

12 月

465.0010,43933,1252,937,660

199,382

416.0

58

%

1月

0.00

0

0

0

0

0

0

%

2月

0.00

0

0

0

0

0

0

%

3月

0.00

0

0

0

0

0

0

%

計

1796.25

22,142

114,541

10,382,810

508,646

207.9

25

%

ノード
利用率

※ システム更新に伴い、システム D は 12 月をもってサービス休止

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 73

システム E 運転状況 ( 2016 年 10 月 ～ 2017 年 3 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]2016/10/01

0:00

2016/10/03

9:00

57.00

2016/10/10

6:30

2016/10/12

9:45

51.25

2016/11/07

9:00

2016/11/08

9:30

24.50

2016/12/19

9:00

2016/12/22

10:00

73.00

2017/03/30

9:00

2017/04/01

0:00

39.00

障害発生日時2017/01/08

サービス再開日時
2:15

2017/01/08

ノード
利用率

15:40

ダウン時間[h]13.42２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

635.75627

3,759

1,420,300

22,752

299.8

71

%

11 月

695.502,690

11,305

1,648,370

37,664

290.4

81

%

12 月

671.002,604

26,646

1,818,630

61,915

356.1

66

%

1月

685.75

2,635

54,161

2,593,210

91,561

481.2

74

%

2月

672.001,006

21,164

1,088,230

24,537

481.2

33

%

3月

705.003,923

67,1952,342,000

82,688

481.1

65

%

計

4065.00

13,485

184,230

10,910,740

321,118

398.3

65

%

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

74 − Academic Center for Computing and Media Studies, Kyoto University

汎用コンピュータシステムのサービス状況
1.

ホスティング・ホームページサービス利用状況

VMホスティングサービス
650
600

ストリーミングサービス
607

602

599

595

591

589

WEBホスティングサービス

550
500
450
400
350
300

280

279

277

277

281

277

250
200
150
100
50
0

2

2

0

平成28年10月 平成28年11月 平成28年12月

0

0

0

平成29年1月

平成29年2月

平成29年3月(平成28年10月から平成29年3月)

大型計算機システム利用承認件数について

平成29年3月末現在、大型計算機システムの利用件数は、3,015件となっています。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 75

資料

別表１ スーパーコンピュータシステム
提供サービス
利用負担額
コース

タイプ

セット

エントリ

タイプA
タイプB
タイプC
タイプE

基本
基本
基本
基本
基本
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位

パーソナル

タイプA1
タイプA2
タイプA3
タイプB1
タイプB2
グループ

タイプB3
タイプC1
タイプC2
タイプE1
タイプE2
タイプE3
タイプA
タイプB

大規模ジョブタイプCタイプE
専用クラスタ-ストレージ容量追加
ライセンスサービス

システム バッチ

12,600 円/年
100,000 円/年
100,000 円/年
100,000 円/年
100,000 円/年
230,000 円/年
115,000 円/年
276,000 円/年
69,000 円/年
690,000 円/年
345,000 円/年
240,000 円/年
120,000 円/年
288,000 円/年
72,000 円/年
720,000 円/年
360,000 円/年
150,000 円/年
150,000 円/年
180,000 円/年
90,000 円/年
280,000 円/年
140,000 円/年
336,000 円/年
84,000 円/年
840,000 円/年
420,000 円/年
23,000 円/週(7日）
11,500 円/週(7日）
24,000 円/週(7日）
12,000 円/週(7日）
15,000 円/週(7日）
7,500 円/週(7日）
28,000 円/週(7日）
14,000 円/週(7日）
720,000 円/年
360,000 円/年
10,000 円/年
20,000 円/年

システム資源

経過時間 ストレージ
無料
（時間）
（TB)
利用者数

最大1ノード相当((36コア、128GBメモリ)×1)
共有
1
最大4ノード相当((68コア、16+96GBメモリ)×4)
共有
最大4ノード相当((36コア、128GBメモリ)×4)
共有
168
最大1ノード相当((72コア、3072GBメモリ)×1)
共有
共有 最大2ノード相当((10コア、32GBメモリ+ 1MIC)×2)
4ノード((68コア、16+96GBメモリ)×4)
優先
2ノード((68コア、16+96GBメモリ)×2)
8ノード((68コア、16+96GBメモリ)×8)A
準優先
2ノード((68コア、16+96GBメモリ)×2)
8ノード((68コア、16+96GBメモリ)×8)
占有
4ノード((68コア、16+96GBメモリ)×4)
4ノード((36コア、128GBメモリ)×4)
優先
2ノード((36コア、128GBメモリ)×2)
8ノード((36コア、128GBメモリ)×8)
B
準優先
2ノード((36コア、128GBメモリ)×2)
8ノード((36コア、128GBメモリ)×8)
占有
336
4ノード((36コア、128GBメモリ)×4)
1ノード((72コア、3072GBメモリ)×1)
優先
1ノード((72コア、3072GBメモリ)×1)
C
2ノード((72コア、3072GBメモリ)×2)
準優先
1ノード((72コア、3072GBメモリ)×1)
4ノード((10コア、32GBメモリ+ 1MIC)×4)
優先
2ノード((10コア、32GBメモリ+ 1MIC)×2)
8ノード((10コア、32GBメモリ+ 1MIC)×8)E
準優先
2ノード((10コア、32GBメモリ+ 1MIC)×2)
8ノード((10コア、32GBメモリ+ 1MIC)×8)
占有
4ノード((10コア、32GBメモリ+ 1MIC)×4)
8ノード((68コア、16+96GBメモリ)×8)A
占有
4ノード((68コア、16+96GBメモリ)×4)
8ノード((36コア、128GBメモリ)×8)
B
占有
4ノード((36コア、128GBメモリ)×4)
168
2ノード((72コア、3072GBメモリ)×2)
C
占有
1ノード((72コア、3072GBメモリ)×1)
8ノード((10コア、32GBメモリ+ 1MIC)×8)
E
占有
4ノード((10コア、32GBメモリ+ 1MIC)×4)
8ノード((36コア、128GBメモリ)×8)
B
4ノード((36コア、128GBメモリ)×4)
ストレージ容量10TBの追加につき可視化ソフト（AVS,ENVI/IDL)およびプリポストウェアの1ライセンスにつき
B
A
B
C
E

0.2
2.0
2.0
2.0
2.0
16.0
8.0
19.2
4.8
32.0
16.0
16.0
8.0
19.2
4.8
32.0
16.0
16.0
16.0
19.2
9.6
16.0
8.0
19.2
4.8
32.0
16.0

-

-

-

32.0
16.0

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコースを、
年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除した後、利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． 大型計算機システムの全ての利用者は、上記表のサービスの他、次のサービスを受けることができる。
1） 大判プリンタサービス
2） その他、大型計算機システムが提供するサービス、機器の利用
3． 上記表の大規模ジョブコース、ストレージ容量追加、ライセンスサービスの申請には、スーパーコンピュータシステムの利用者であることが必要である。
4． 「共有」： 当該カテゴリのユーザ間で一定の計算資源を共有するベストエフォートのスケジューリングを行う。
「準優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/4の計算資源が確保されることを保証する。
「優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/2の計算資源が確保されることを保証する。
「占有」： 稼働状況によらず記載値の計算資源が確保されることを保証する。
5． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。
6． グループコース及び専用クラスタコースの利用者番号は利用者あたり年額5,000円を負担することで追加できる。
7 ． 機関・部局定額制度
他機関又は学内における部局（『国立大学法人京都大学の組織に関する規程』第3章第2節から第11節で定める組織をいう。）の組織が、その組織単位でグループコースサービスを利用申請する場合の利用負担額は、別表１に規定する額の1.5倍の額とする。
なお、利用負担額が年額150万円未満の場合は100人、年額150万円を超える場合は、150万円毎に100人までの利用を認める。
ストレージは、1.5倍の容量とする。
8 ． スパコン連携サービス
学術情報メディアセンターのスーパーコンピュータシステムと密な連携により、学内における部局の組織が計算サーバ等を設置する
場合、下記の負担額を支払うものとする。
冷却方式
水冷
空冷

利用負担額
11,300 円/月
13,200 円/月

利用負担額算定単位
水冷冷却方式の計算サーバ等の定格電力 1kWにつき
空冷冷却方式の計算サーバ等の定格電力 1kWにつき76 − Academic Center for Computing and Media Studies, Kyoto University

8
4
16
4
16
8
8
4
16
4
16
8
8
8
16
8
8
4
16
4
16
8

16
8

別表２ 汎用コンピュータシステム
区分

利用負担額

仮想サーバホスティングサービス

36,000円／年

単位
1仮想サーバにつき備考
1． 利用負担額は、総額表示である。
2． 上記表の仮想サーバホスティングサービスを利用するためには、スーパーコンピュータシステムの利用者であること。
3． 1仮想サーバに割当てるシステム資源は、CPU：2コア、メモリ：4GB、ディスク：100GBである。
4． 仮想サーバホスティングサービスにおいて、下記の負担額を支払うことによりCPU、メモリ、ディスクを増量することができる。
区分

利用負担額

単位

CPU増量

3,000円／年

2コアにつき(最大8コアまで)

メモリ増量

3,000円／年

4GBにつき(最大64GBまで)

ディスク増量6,000円／年100GBにつき(最大1,000GBまで)

5． 利用負担額は、当該年度（4月から翌年3月まで）の利用に対して年額として算定するが、年度途中から利用を開始する場合には月数に応じて減額する。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 77

※民間機関利用
別表３ スーパーコンピュータシステム
システム

A

B

C

E

経過時間
（時間）

ストレージ
（TB)

無料
利用者数

8ノード(68コア、16+96GBメモリ)×8)

336

19.2

16

1,104,000 円/年

12ノード(68コア、16+96GBメモリ)×12)

336

28.8

24

1,656,000 円/年

16ノード(68コア、16+96GBメモリ)×16)

336

38.4

32

2,208,000 円/年

8ノード(36コア、128GBメモリ)×8)

336

19.2

16

1,152,000 円/年

12ノード(36コア、128GBメモリ)×12)

336

28.8

24

1,728,000 円/年

16ノード(36コア、128GBメモリ)×16)

336

38.4

32

2,304,000 円/年

2ノード(72コア、3072GBメモリ)×2)

336

19.2

16

720,000 円/年

3ノード(72コア、3072GBメモリ)×3)

336

28.8

24

1,080,000 円/年

4ノード(72コア、3072GBメモリ)×4)

336

38.4

32

1,440,000 円/年

8ノード(10コア、32GBメモリ+ 1MIC)×8)

336

19.2

16

1,344,000 円/年

12ノード(10コア、32GBメモリ+ 1MIC)×12)

336

28.8

24

2,016,000 円/年

16ノード(10コア、32GBメモリ+ 1MIC)×16)

336

38.4

32

2,688,000 円/年

システム資源利用負担額

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコースを、年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除した後、利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。

78 − Academic Center for Computing and Media Studies, Kyoto University

全国共同利用版広報・Vol.15(2016)総目次［巻頭言］
Vol.15，No.1 号の発刊にあたって ................................................................................................................1-1
Vol.15，No.2 号の発刊にあたって ................................................................................................................2-1［新スーパーコンピュータサービス開始］
新スーパーコンピュータ利用ガイド .............................................................................................................1-3
新スーパーコンピュータ Camphor2 のベンチマーク評価報告 ...............................................................1-15［研究会開催報告］
The 35th JSST Annual Conference（JSST2016）開催報告..................................................................1-20［スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告］
液膜内振動マランゴ二対流のメカニズム解明 ..............................................................................................2-3
乱流相分離制御による自己組織化構造の能動的選択 ..................................................................................2-6
有機系太陽電池への応用を指向した光機能性有機分子の構造と電子構造の相関解明 ............................2-9
バイオフィルム成長に関する三次元モデルの構築 ....................................................................................2-11
間欠性ダイナミクスの不安定周期軌道解析................................................................................................2-14
運動量輸送と熱輸送の非相似性の最大化....................................................................................................2-16高分子特性を考慮した界面活性剤添加溶液の抵抗低減効果.....................................................................2-18
周囲気体を考慮した HZ 液柱内温度差マランゴニ対流場の数値解析 .....................................................2-21
高プラントル数を用いた HZ 液柱内温度差マランゴニ対流の数値解析―カオス・乱流化過程の解明に向けて―....................................................................................................2-23
パターン学習を用いた経験的セル・オートマトン構成法とそれによる偏微分方程式モデルの模倣 ...2-25
共有結合性有機骨格構造体(COF)の電荷移動シミュレーション..............................................................2-29
理想化モデルを用いた，竜巻を模した渦の力学に関する研究.................................................................2-31
オフロード路面における車輪の走行性能解析―惑星ローバ用車輪におけるラグ周辺部分のミク口な砂の挙動解析―..................................................2-34［プログラム高度化共同研究報告］
非線形有限要素法による鉄筋コンクリ一卜の大規模破壊シミュレーション..........................................2-36
津波流体力による構造物の破壊過程を考慮した構造流体連成解析 .........................................................2-41
産業集積の空間的同期と都市規模・空間パターンに関する研究 .............................................................2-45「メッシュ状透過性平行平板間乱流の直接数値シミュレーション」に関するプログラム高度化 .........2-48
乱流輸送コード GKNET への運動論的電子の実装...................................................................................2-52［スーパーコンピュータ共同研究制度（大規模計算支援枠）研究報告］
地震前兆における山岳地表面を伝搬する電磁プラズマ波の特性解明 .....................................................2-56「メッシュ状透過性平行平板間乱流の直接数値シミュレーション」に関する大規模計算.....................2-60

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 79

［会議参加報告］
SC16 参加報告 ...............................................................................................................................................2-64［サービスの記録・報告］
スーパーコンピュータシステムの稼働状況................................................................................................1-22
汎用コンピュータシステムのサービス状況................................................................................................1-27
スーパーコンピュータシステムの稼働状況とサービスの利用状況 .........................................................2-65
センター利用による研究成果（平成 27 年度）..........................................................................................2-71［資料］
大型計算機システム利用負担金 別表...............................................................................................1-28,2-74
全国共同利用版広報・Vol.14(2015)総目次 .................................................................................................1-31
サービス利用のための資料一覧 ..........................................................................................................1-33,2-77［編集後記］
編集後記、奥付 .....................................................................................................................................1-34,2-78

80 − Academic Center for Computing and Media Studies, Kyoto University

―1.

サービス利用のための資料一覧―

スーパーコンピュータシステム・ホスト一覧
システム A：camphor.kudpc.kyoto-u.ac.jp



システム B・C：laurel. kudpc.kyoto-u.ac.jpシステム B（SAS 利用時）
：sas.kudpc.kyoto-u.ac.jp

システム E：camellia.kudpc.kyoto-u.ac.jp※ ホストへの接続は SSH(Secure SHell) 鍵認証のみ、パスワード認証は不可

2.

問い合わせ先 & リンク集


情報環境機構のホームページ
http://www.iimc.kyoto-u.ac.jp/



学術情報メディアセンターのホームページ
http://www.media.kyoto-u.ac.jp/



スーパーコンピュータシステムに関する問い合わせ先

利用申請などに関する問い合わせ先【情報環境支援センター】
E-mail : zenkoku-kyo@media.kyoto-u.ac.jp / Tel : 075-753-7424
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/

システムの利用など技術的な問い合わせ先【スーパーコンピューティング掛】
E-mail : consult@kudpc.kyoto-u.ac.jp / Tel : 075-753-7426
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/contact.html

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.16 No.1（ Jan. 2018）− 81

編集後記

編 集 後 記No.1 の発刊が遅くなってしまいましたが，年度内に No.2 の発刊も予定しております．No.2 では，2017
年 10 月よりスパコンのサブシステム B/C で提供を開始したバーストバッファ（高速ストレージ）に関する特集を予定しております．最新の機能であり，やや使い所が難しいところもありますが，うまく利用できればファイルアクセス性能を大きく改善できる可能性もありますので，ご参考にしていただ
ければと思います．
(副部会長)

京都大学学術情報メディアセンター全国共同利用版広報 Vol．16，No．1
2018
年 １月
1 月３０日
xx発行
日発行
２
０１８年
広報編集部会
編集者 京都大学学術情報メディアセンター

深沢 圭一郎（部会長）

全国共同利用版広報編集部会

平石 拓（副部会長）

発行者 〒606-8501 京都市左京区吉田本町

水谷 幸弘

南雲 円

京都大学学術情報メディアセンター

熊谷 真由美尾形 幸亮

Academic Center for Computing and Media Studies

高見 好男

Kyoto University
Tel. 075-753-7400
http://www.media.kyoto-u.ac.jp/
印刷所 〒616-8102 京都市右京区太秦森ヶ東町 21-10
株式会社エヌジーピー

82 − Academic Center for Computing and Media Studies, Kyoto University

表紙デザイン：谷 卓司（ティアンドティ･デザインラボ）