Vol.14, No.2 2 15

センター共同研究報告（平成26年度）【巻頭言】「Vol.14，No.2 号の発刊にあたって」牛島 省 【スーパーコンピュータ共同研究制度（若手研究者奨励枠）】山本 卓也◎茂木 孝介◎Suddhasattwa Das, Yoshitaka Saiki, Evelyn Sander, James A Yorke◎鬼頭（西岡)宏任，安藤 耕司◎中野 直人◎高木 洋平◎水口 尚◎中本真義◎梅山 有和◎原田 浩平，谷口 貴志◎中村 翔一 【プログラム高度化支援事業】中畑 和之◎野田 利弘◎車谷 麻緒◎森 知也 【スーパーコンピュータ共同研究制度（大規模計算支援枠）】山本 義暢

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

牛島 省

本号では、学術情報メディアセンターにおいて実施された平成 26 年度の共同研究報告について特集いたします．平成 26 年度は，若手研究者の支援事業に関して 11 件，プログ
ラム高度化支援事業について 4 件，大規模計算の支援事業について 1 件の共同研究が実施されました．計算科学の広範な分野にわたる研究報告が掲載されていますので，関連する分野の研究者や学生の方の参考になると思います．
なお，これら３種類の共同研究制度のうち，「若手研究者奨励枠」は，平成 27 年度に実施された追加募集から「若手・女性研究者奨励枠」と変更され，(1)40 歳未満の若手研究者(学生を含む, 性別は問わない)，(2) 女性研究者 (年齢は問わない)，の方を対象とする奨励制度となり，支援対象者の範囲を広げております．また，平成 28 年度より，この「若手・女性研究者奨励枠」は，文部科学大臣の認定を受けた JHPCN (学際大規模情報基盤共同
利用・共同研究拠点) の活動の一環となり，スーパーコンピュータを利用することで学術的にインパクトがある成果を創出できると期待される課題に対して，計算機利用負担金の全部または一部を本センターが負担しますので，ぜひ積極的にご利用ください．今回報告のあった 11 件には各種の流体計算，電荷移動問題，時系列データ解析，有機分子材料計算など多様な研究課題が含まれており，また，東北大学，東京理科大学，一橋大学，名古屋大学，京都大学，大阪大学，琉球大学と全国の若手研究者との共同研究となりました．
プログラム高度化支援事業では，スーパーコンピュータをグループコースまたは専用クラスタコースでご利用の研究グループを対象に，プログラムの制御構造・データ構造の改良による性能チューニングや並列化手法の改良，問題分割・負荷分散方式などの改良による大規模計算プログラムの高度化・高性能化を支援しています．平成 26 年度は，有限積分法 (FIT) による空洞を含む固体中の波動解析，地盤の有限変形解析，鉄筋コンクリー
ト中のひび割れ進展解析，都市の人口規模・空間分布における秩序形成解析を行うプログ
ラムの高度化を支援させていただきました．メニーコア (MIC) 環境でのプログラム高度化や，連立１次方程式の並列解法，経済学分野の問題を解析するためのハイブリッド解法など，様々な高度化・高性能化が行われており，数値計算研究者にとって最適化の良い指標となると思います．
また，大規模計算の支援事業に関しては，スーパーコンピュータをグループコースまたは専用クラスタコースでご利用の研究グループを対象に，例えば最大で 128 ノード × 2
週間の大規模ジョブコース利用を無償で実行できる支援を行っています．平成 26 年度は，
高レイノルズ数チャンネル内の乱流解析の大規模計算の支援が行われ，その研究成果が報告されています．
今後も皆様の研究，教育にご活用いただけるようにセンター教職員も尽力していきます
ので，ご利用・ご支援のほど、よろしくお願い申し上げます．

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

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

三次元非定常マランゴニ対流の解明
Ä wXħĩñð÷ôáęoµ»1ºģÀõ§ Navier-Stokes §ġþ§

緒言

Energy §1çÓÛêįT!HKP7AÚÓÛêįTÜ¦tģ)Ý¦tģ æêÅ]Q¿àÅ]TU>395!įÅ]Q.ĉį©ptģ æÚ/ íæ.

Ü¦1ėTU>395 Ü¦ +Ú/!HKP7AÚ(ĕü![UHKP7AÚ1²Ě¾ä(ć¾[1-2]¾[3-4])1zß)"/v

ĭ¼u·Ö!øĦ 2 Ðû¦ÿªĊĦÖÂĦ` aċEO:9ù ¯çį(ĤČ.

 2 Ðû¦ Backward Ö1çĮÍĠÊF,v` ċĠ!Ú` VĝÁ°+R9 ·Ö.'Laplacian ıĮïYĊÑ

-½eďÆ £¥ ð÷ç,/.

1Ĉ®Č .[8]Z !ĪîïYĊÑ1ĈĜ¤½eďÆ!CLR(4REP); Over-relaxed approach[5]1çđúÊ Ìë1 Fig. 1 (c)òVùĦī

R9 Ùç {"/-Ú`ďÆĮÍĠÊFR9ÃĨ`öÖĥí/

HM<DO=5ÍĠÊ1çđúÊďg

OpenFOAM[5-6] ¥%-.½eďÆ!

¦!hĈð÷ ďg¦1zĂ(Fig. 1 (c)W

V&V(Verification and Validation)Ĕc Į¢ĤČ

½ Nr = 35, Nz = 50, Nφ = 35)Z ç. ç6R@!ì-CLR;R9

OpenFOAM

BR8JP! 2.1.1 .!Ąm V&V Ĕc1Ĉ0.1«Ä!CLR;R9.

3

結果と考察

OpenFOAM 1çđú ďÆû¦ Ô',/.Ú`V° ĳËē1ĈË Ú/ Î1 Fig. 2 òĄèĉįÕēĳTĒ HKP7AÚ ÛÉ

ĴÜ>395 ,_Ü>395f~ HK

.V° đú1bç

P7AÚ íæ Ú/!ÛÉWĢ1Ğ
k#ĴÜ>395f¶. Ú/!i

2

数値解析手法îÛÉj HKP7AÚ â­.
 9?NR?ÛÉ | ÛÉWĴ¿į

½eďÆĲ1 Fig. 1 (a), (b)òđúĘğ¦F5?M1 Fig. 3 ò×½ 2

IR@Ĭ^ IMA-2 %',/Ëēíæ-/!¾äe[7]Są.

ĳ[7]1bç.TUqàQlx>3951

 ĪîĮïYĊÑ1Ĉ|(corrected,ģāÛÉ1± ÛÉ ªã!9?N

Fig. 3 (a) and (b))Ĉ|(uncorrected,

R?(a)ªĵ¹ĐĎ 60 ¦Ķ(b) |1ĘFig. 3 (c) and (d))Ú/ Î ì.ĮïYĊÛÉ 29G5?Ò(d/R)!1.2.Ñ1Ĉ|!íæ$ ! ě%ÛÉ ĉį!ÓÛêį.Û`Ģč,/ĮïYĊÑ+.½eĖ

&1đúÛ`Prandtl ½

0.01



ě1

¨ ĚĂ,/.

Ú`1±ÞĤpÅ]ÓÛêį sîª

 ªÛÉ | Ü¦tģ m¡1 Fig. 4 ò1\

Īî laplacian ı ĊÑ1.|(corrected)

2 − Academic Center for Computing and Media Studies, Kyoto UniversityĊÑ1|(uncorrected)vmĊÑ.|Methods for Fluid Dynamics (Springer, Berlin, 2002).

(limited 0.5) ýÇ1//ò.ĊÑ1
%!vmĊÑ1.|! Ferziger and
Peric [8].+¸sď(ĀÏÎ Ü¦tģm¡) íæ9?NR?ÛÉ ýÇ|0Ă´.ïYĊÑ1Ĉ ĮïY°1.®Č .*hĈð÷ù[7]!
%-®Č!/À~ đúÊďg

Fig. 1 Schematics of (a) straight and (b) deformed¦1.®Č .À~Ê½1qliquid bridges, and (c) the framework of the numericalđú1ĈĮïYĊÑÃ-đú1Ĉ

grid.

hĈð÷}Î IR@ íæ¸sďù(č,/ 4

結言Ä!CLR;R9. OpenFOAM

1çHKP7AÚ

V°.

V&V Ĕc1ĈOpenFOAM !ĮÍĠÊF
R9 ÃĨ`öÖ 6R@.'đúÊ
ĮïY°+-laplacian ıĮïYĊÑ1çFig. 2 Velocity distribution in a cross section in theperpendicular to the liquid bridge in (a) straight and (b)deformed liquid bridge.

.®Č .  ĮïYĊÑ ĮåîÚ/1¨ Ě| . m Ĵû¦đú ®Č.|đúÊ Įï
Y°1 .Ø³.®Č .ÄđúĘ!hĈð÷!%-®Č
/À~ Êďg¦

OpenFOAM

!®Č. m 

引用文献[1] L. E. Scriven and C. V. Sternling, Nature, 187,
186-188 (1960).
[2] H. Kuhlmann, Thermocapillary convection inmodels of crystal growth (Springer, Berlin, 1999).
[3] N. Imaishi, Int. J. Microgravity Sci., 31, S5-S12

Fig. 3 Velocity vectors in a cross section of horizontal
plane in the middle height: (a) corrected at 10 s, (b)corrected at 30 s, (c) uncorrected at 10 s, (d)uncorrected at 30 s.

(2014).
[4] S. Matsumoto et al., Int. J. Microgravity Sci., 31,
S13-S35 (2014).
[5] H. Jasak, Ph. D. Thesis of the University of London,
Imperial College of Science, Technology and Medicine(1996).

Fig. 4 Temperature gradient along the free surface: (a)uncorrected, (b) limited 0.5 and (c) corrected.

[6] OpenFOAM (see http://www.openfoam.com/)

謝辞[7] V. M. Shevtsova, J. Cryst. Growth, 280, 632-651

 Äð÷n¦ĵĆ·ð÷ărÈĶ1Ùç(2005).

İ1 1d-y ¬óéT

[8] J. H. Ferziger and M. Peric, Computational.

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

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

高プラントル数流体を用いたHZ液柱内温度差マランゴニ対流の数値解析
日欧共同宇宙実験に向けて
茂木孝介
東京理科大学大学院理工学専攻機械工学研究科

1. 緒言
気液界面上に温度分布が存在する場合、表面張力の温度依存性によりマランゴニ対流が生じる。マランゴニ対流は微小重力環境下やマイク
ロナノスケール環境下で顕著化し、本質的な問題となりうる。
微小重力環境下は浮力対流が抑制されるため高品質な材料生成を行う上で理想的な環境である一方で、材料生成過程においてマランゴニ対流の不安定性により生じる振動流が悪影響を及ぼす事が知られている。そのため、マランゴニ
対流の不安定性の理解、制御の観点から研究が進められてきた。
私の所属する研究室は、国際宇宙ステーショ
ン日本実験棟「きぼう」における流体物理宇宙実験に共同研究グループとして参画している。
今回、2008~2013年にかけて実施された宇宙実験「マランゴニ対流におけるカオス・乱流とその遷移過程」の結果の理解を目的として、数値シミュレーションを行った。微小重力実験では背の高い液柱の長時間安定した実験が実施できる事が魅力の一つである。そのため微小重力実
験では、地上では実現不可能な長液柱を作成し実験を実施した。
本研究は微小重力環境下におけるマランゴニ
対流の不安定性の理解と対流場の遷移過程を明らかにする事を目的としている。、微小重力実
験の結果と比較を行うため、長液柱を対象とし、液柱のアスペクト比を変化させ、対流場が不安定となる臨界マランゴニ数への依存性を調べた。

2. 計算モデル
高品質な材料生成法の一つであるF l o a t i n g zone法の代表的なモデルであるHalf-zoneモデル(HZモデル、図１）を研究対象とする。HZモデ
ルはﬂoating-zone法の下半分の領域を模擬した
簡易モデルであり、これまでのマランゴニ対流の不安定性の研究で慣例的に使われてきた。HZ
モデルにおいて、Maを大きくしていくと軸対称
定常流が回転波に対して不安定となり、非軸対称周期流へと遷移する。
液柱の高さをH、半径をRとし、アスペクト比をΓ=H/Rと定義する。上下ロッド間に液柱を形成し、上部ロッドを加熱、下部ロッドを冷却する事で液柱表面に温度分布を付加し、マランゴ
ニ対流を発生させる。液柱は表面変形がない事を仮定し、液柱の周囲気体の影響を無視する。
また、無重力であるとする。
支配法的式はN a v i e r- S t o k e s方程式、連続の式、エネルギ方程式である。液柱の物性値としてプラントル数（Pr）を 28として計算を行った。

マランゴニ数（Ma）はマランゴニ対流の強さを表す無次元数であり温度差に起因する。また、
微小重力実験の目的の一つである、カオス・乱流化過程を調べるため、 M a を徐々に上昇させ対流場の遷移過程を調べた。

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

図１ HZモデルの図

3. 結果条件で異なる M a 数で遷移する。この領域では臨界マランゴニ数、対流場遷移過程のアスペ
クト比依存性を示したグラフを図２に示す。アスペクト比の違いによって、二つの異なる対流と、共存する場合とがあるようである。

T W 1 aとT W 1 bそれぞれが独立に存在する場合

場構造（TW1a,TW1b）が観察された。この対
流場は無限に長い液柱内対流場の線形安定性解
析[1]の結果から予測されており、後に有限液柱である微小実験[2]によって観察されたものと同一のものである。二つの対流場はアスペクト比に依存して現れる。T W 1 aは長液柱にのみ存在
し、微小重力実験でのみ観察する事が出来る。グラフを見ると、アスペクト比を増加してていく時、TW1aからTW1bに臨界モードが切り替わると臨界値が一度減少し、その後徐々に上昇していく。

図３ 液柱表面の温度変動分。

4. 結言
Pr = 28の流体、宇宙実験で実施された長液柱を対象とし数値シミュレーションを実施した。
臨界マランゴニ数とアスペクト比の依存性、対流場の遷移過程を調べた。その結果、宇宙実験同様に異なる対流場が観察され、臨界マランゴ
ニ数のアスペクト比依存性を明らかにした。
また、対流場の遷移過程はアスペクト比によっ
図２ 臨界マランゴニ数と対流場遷移過程のアて異なり、高アスペクト比（2.5以上）では対流

スペクト比依存性。

場の切り替わりにおいてヒステリシスが観察された。

臨界マランゴニ数を超えて、さらにMaを大き

今後は新規テーマであり、2018年から予定さくしていくとTW1aからTW1bに遷移することを明らかにした。その時の温度の変動分を図３に示す。マランゴニ対流の不安定性によって現れる図の様な回転波をHydrothermal waveと呼ぶ。
表面の温度分布から異なる対流場パターンを持つ事が分かる。(a),(c)は、「く」の字に曲がった温度分布を持つのに対し、(b),(c)はライン状の構造をしている。
さらにマランゴニ数を大きくするとT W 1 aとTW1bが混在する準周期流(グラフではQ-Pとした)へと遷移する。TW1aからTW1bへの遷移領
域はヒステリシスをもち、Ma数上昇条件、減少

れている日欧共同宇宙実験に向けた予備計算を実施していく予定である。具体的には、液柱の周囲対流を制御する事によって臨界値の制御を目指す。

5. 参考文献[1] I. I. Ryzhkov, ”Thermoapillary instabilities inliquid bridges revisited,’’ Phys. Fluids 23(8),082103(2011)[2] T. Yano, K. Nishino, H. Kawamura, I. Ueno,
and S. Matsumoto, ’’Instability and associated roll
structure of Marangoni convection in Prandtl
number liquid big with large aspect ratio,’’ Phys.
Fluids 27, 024108 (2015)

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 5スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告

ࡵ࠵ՃՕஞ⎬ୀ⏐⏌ үࠦച‫ה‬ຟ⎯‫ۿ⎯ޔظ‬ਙшSuddhasattwa Das1 , Yoshitaka Saiki2∗ , Evelyn Sander3 , James A Yorke1
1
University of Maryland, 2 Hitotsubashi University , 3 Geroge Mason University.

⎰⎙ ⏂⎬

1

⎨⎃⏌ . ຢӁ‫ ⎩ذ‬f ⎍ C ∞ ⎨ , w ⏓ ࣘ⎨ߟ⎘ ⎠ड़߸⏓ ⏃ ⎥ຢӁ‫⎯ذ‬इ⏌ ೮⎅⎰ , Α଑࠵ՃՕஞ (೑ஓ୅⏓ ԥ⏁ ), ࡵ࠵ՃՕஞ, ⏧ ⏦⏵ ⎯‫ ⎬⎥މ‬੒⎎ ⎐ ഇຳ⎖ ⏍⏌ . ⎞ ⏍⏊ ⎰ຶ⎉ ⎱␦ ⏞ ␓␊ ␑ ‫⎯ ⎪ ⎫ॅޱ‬ച‫ה‬ຟ⎬ ⏉ ⎤ ⎧ ஦૎⎦⎒ ⎖ ⏍⏌ ⎍ , ⎞ ⎯࠴ਖ⎰ಒ࣢⎬ ઢ⎐ , ߹‫ ⎬ݘ‬Օஞ⎍ิ⎉ ⏊ ⏍⎠ ⎩⎎ ⎬ ஦⎬ ࡵ࠵ՃՕஞ⏓ ౜ധ⎚⏌ ⎔ ⎩ ⎰௔⎘ ⎅ . ⎿ ⎠ , ࡵ࠵C ∞ ࡚⏀ԟॅ⎩ ⎘ ⎧ , ۖ⎬ ࠚӺ⎯щ଑⏓ ю⎉ ⏌ ⎩ , ௥Τ⎯ߢীॅ m ⎬ ਻⎘ ⎧ , ⎃ ⏌ Cm > 0 ⎍ਣ‫ ⎧ ⎘ݚ‬,
�
�
�
�
�
�WBN (f )(x) − f dµ� ≤ Cm N −m��⎨⎃ ⏌ ⎔ ⎩ ⎍ߟ⎖ ⏍⏌ [1].

ՃՕஞ⏓ ஦૎⎦⎒ ⏌ Ѿୃॅ⏓ ‫ۿ‬३୞⎨ஔ଑⎚⏌ ⎔ ⎩ ⏃ ஔโ⎯๲ศ⎨‫ܟ‬௔⎨⎃⎤ ⎠ . ඕ٥֋⎨⎰ࡵ࠵ՃՕஞ⎬ୀ⏐⏌⎔ ⏍⏊ ⎯ຟ⏓ ‫ۿ‬ਙ⎨३୞ຜ⎐ ‫ ⏌⎚ޔظ‬൘൚⏓ ଚΕ⎘ ⎠ . ⎞⎔ ⎨ , क़‫މځ‬ਹ෻੕⎯⎃⏌ ࡵ࠵ՃՕஞ⎯Ѿୃॅ⎯‫ຶ ⏓ޔظ‬⎬ ⎘ ⎧ , ⎞ ⎯ࠡ൚⏓ ࣅѻ⎚⏌ .

Ιч⎨⎰ , ಘഇ൘ଟߧ (4) ⎨⎃⏊ ⏐⎖ ⏍⏌ क़‫މځ‬ਹ෻੕[2] ⎯ࡵ࠵ՃՕஞ⎬ ਻⎘ ⎧ , ࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ ⏓ה‬Ыๆ⎚⏌ .

dq1 /dt࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ה‬

2क़‫މځ‬ਹ෻੕

3␌⍬⏯ ␑ ⎯⏤␧⏰⍬␅ ଑๲⎰ , ⏤␧⏰⍬␅ ଯ⎫ ຢӁ‫⎯ذ‬Օஞ (xn ) ⎬ ਻⎘ ⎧ ଑զ⎖ ⏍⏌ ъ⏊ ⎌⎯ຟ f ⎬ ਻⎘ ⎧ , ⎞= p1 + q2 ,dq2 /dt = p2 − q1 ,−3
dp1 /dt = p2 − µ(q1 − 1 + µ)d−3
m − (1 − µ)(q1 + µ)dp ,dp2 /dt−3= −p1 − µq2 d−3
m − (1 − µ)q2 dp .

(4)−1⎯␌⍬⏯ ␑ ച‫ ה‬BN (f )(x) = ΣN
n=0 f (xn )/N ⎰ , Օஞૠ

N → ∞ ⎨‫׸‬ԝച‫ ה‬f dµ ⎬࠴ਖ⎚⏌ ⎩ ⎅⎇ ⏃ ⎯⎨⎃⏌ . ⎔ ⎔ ⎨ ,੫⎘ , µ ⎰ຢӁ‫⎯ذ‬೑ബҷ๻ਗ୞⎨⎃ ⏌ . ѱゖ ⎰ , ⎔ ⏍⎬dm = ((q1 − 1 + µ)2 + q22 )0.5 , dp = ((q1 + µ)2 + q22 )0.5 .
࠹।⏓ ю⎉ ⏌ ⎔ ⎩ ⎬⏉ ⎤ ⎧ , ࡵ࠵ՃՕஞ⎬ୀ⏐⏌ үࠦ⎯ຟ⎨⎃ ⏋ , µ = 0.1 ⎩ ⎚⏌ . ‫ ␯ ⏞␇ ␄␧ ␛␋⎯ذ‬H ⎰ ,⎯‫ۿ ⏓ޔظ‬ਙш⎖ ⎜⏌ .
�࡚⏀ԟॅ w : R → [0, ∞) ⏓ ෌‫ځ‬Ѿಘഇя௾⎨⎌⎥ (0, 1)⎨ w > 0, ⎞ ⏍Ιҙ⎨ w = 0 ⎩ ⎫ ⏌ ⏃ ⎯⎩ ⎚⏌ . Ιч⎨−1
H = [(p21 +p22 )/2]+[q2 p1 −q1 p2 ]−[µ d−1
p +(1−µ) dm ] (5)⎨⎃ ⏌ . H ⏓ ⎋ ⏉ ⎞ -2.63 ⎯ખ⎬ ঎଑⎘ ⎧ , q2 = 0 ⎌⎥⎰ , ஦⎬

w(t) :=


�
exp

1
t(t−1)

�dq2 /dt > 0 ⎬ ⏉ ⎤ ⎧ ଑⏂⎠ ␙⏞ ␯ ⏧ ␨ ઐ෣ࣘ⎯๸‫ޏ‬Օஞ,,

0,for t ∈ (0, 1)
for t ∈/ (0, 1).

(1) ⎥⎿ ⏋ 2 ߘ‫ࠃٴ‬ਉ⎬ ⎥⎅⎧ ૛⎺⏌ ⎔ ⎩ ⎬ ⎚⏌ .
Ιч⎨ߟ⎚ॅખ‫ ⎰ޔظ‬4 ఝ३୞⎨⎋ ⎔ ⎫ ⎅ , ॅખঀഇ⎩ ⎚⏌ . ⎔ ⏍⏓ ๆ⎅⎧ f ⎯࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ ה‬WBN ⏓WBN (f )(x) :=N−1�ŵn,N f (xn ),n=0⎫ ⎤ ⎠ . ⎿ ⎠ , ␙⏞ ␯ ⏧ ␨ ઐ෣ࣘ⎯Օஞ⏓ ३୞ຜ⎐ ֆ⏂⏌⎠ ⏂⎬ , (4) ⏓ ॅખঀഇ⎘ ⎧ ઐ෣⏓ Юউ⏌ ૥ঽ⎯୅⏓ ֆ⏂(2) ⎠‫ڤ‬, ࡓਛബॅ q ⎩ ஫๼ബॅ t ⎯ฉӍ⏓ ௡⏍੃⎉ ⎧ , q ⏓22஫๼ബॅ⎩ ⏀⎫⎘ ⎧ ␙⏞␯ ⏧␨ ઐ෣⎬⎢ ⏈ ⎇ ⎪ ஃ੬⎚⏌ ⏉⎩ ଑զ⎚⏌ . ੫⎘ ,ŵn,N =

⎰␧␯ ⏮ =⏫ ⏿ ⏻ ⎯ 8 ߘ‫[ ߧڷ‬3], ߗԝ‫ ⎰⏀܌‬10−5 ⎨⎋⎔⎇ ⎬ ␧␯ ⏮ =⏫ ⏿ ⏻ ⎯ 8 ߘ‫ ⏓ߧڷ‬1 ⏵ ␂ ⏿ ␓ഇ଱ๆ⎘ ⎠ .

w(n/N )−1
ΣN
j=0 w(j/N )(3)⎔ ⎯⏉ ⎇ ⎬ ⎘ ⎧ ண⏊ ⏍⎠ ࡵ࠵ՃՕஞ (Օஞૠ N ) ⎯ߗ‫܌‬n ⎬⎋⎒ ⏌ ୅⏓ xn ⎩ ⎘ ⎠⎩ ⎎ , य 1 ⎨ߟ⎚⏉ ⎇ ⎬ (q1 , p1 )ച෣ࣘ⎨ q1 ߫⎩ ⎯⎫ ⎚Һ φn ⏓ ଑⏂ , ߗ‫ ⎍܌‬1 ⎡⎒ छ⏁ඡ⎯Һ⎯छ⏀‫⎯܄״‬ച‫( ה‬Ѿୃॅ) ⏓ ࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ה‬∗E-mail:yoshi.saiki@r.hit-u.ac.jp

6 − Academic Center for Computing and Media Studies, Kyoto University⏓ ๆ⎅⎧ ֆ⏂⎠ .

य 1: ␙⏞ ␯ ⏧ ␨ ઐ෣ࣘ⎯ࡵ࠵ՃՕஞ. ߗ‫ ܌‬n ⎬ ⎋⎒ ⏌ Օ

य 2: य 1 ⎬ߟ⎘ ⎠ࡵ࠵ՃՕஞ⎯Ѿୃॅ ρN ⎯‫( ܴڬ‬Օஞૠஞࣘ⎯୅ xn ⏓ q1 ߫⎌⏊ ⎯⎫ ⎚Һ φn ⎨ਕ⎉ ⏌ .

N ⎩ ⎯਻ಅ). ‫| ܴڬ‬ρN − ρN ∗ | Ƅ (੫⎘ , N ∗ = 400, 000)
⎍ 4 ఝ३୞‫ޔظ⎯ޔظ‬३୞⎯‫ځ‬ҏ, ⎚⎫ ⏐⎢ Ԥ⏂‫⎯ܴڬ‬य 2 ⎰⎔ ⎯‫⎨ޔظ‬ண⏊ ⏍⎠ Ѿୃॅ ρN ⎯২଑⎖ ⏍⏌ ‫ٶ‬වખ (Օஞૠ N ∗ = 400, 000 ⎯ࣟ‫⎯܄‬ખ⎨੐ๆ) ⎌⏊ ⎯‫ڬ‬⏦⍬⏼⍬ 10−30 ଟ୞⎬ ⎿ ⎨‫ॅ ⎩ ⏌⎚ࢢٹ⎍ܴڬ‬ખ‫ࣘޔظ‬⎰࠴ਖ⎘ ⎠ ⎩ ⏀⎫ ⎜⏌ .

ܴ⏓ Օஞૠ N ⎯ԟॅ⎩ ⎘ ⎧ ߟ⎘ ⎠ ⏃ ⎯⎨⎃⏌ . ⎔ ⏍⎌⏊ ,N = 20, 000 ଟ୞⎬ ⎋ ⎅⎧ , 4 ఝ३୞⎯‫ޔظ‬३୞⎯‫ځ‬ҏ⎬ ⎅⎠⏂ , ࠴ਖ⏵␐ ⍬␅ ⎯ܴ⏓ ๆ⎅⎧ ⏧ ⏦⏵Օஞ⎩ ࡵ࠵ՃՕ‫ ⎅פ‬1030 ⎬⎿ ⎨ , ‫ ⏌⎌⏐⎍ ⎩ ⎔ ⏌⎅ ⎧ ⎘ࢢٹ⎍ܴڬ‬. ⎿ ⎠ , ஞ⎯ഇຳ⎬⏃ ๭ๆя௾⎨⎃⏌ [4]. ⎿ ⎠ , ࡚⏀೒⎎ ␌⍬⏯ ␑‫࠴⎯ܴڬ‬ਖ␨ ⍬␄ ⎰ 1/N 15 ଟ୞⎨⎃ ⏌ ⎔ ⎩ ⎍⏐⎌⏌ . N ച‫⎯ה‬଱ๆ൘൚⏅஦ड़⎬ԟ⎚⏌ ࣏‫ݔ‬, क़‫މځ‬ਹ෻੕Ιҙ⎯⎍ࡒഇ੒⎎ ⎐ ⎫⏊ ⎫⎅⎇ ⎢ ⎬ 4 ఝ३୞⎯‫ޔظ‬३୞⎬ஃ੬⎘ຢӁ‫ٌ ⏌⎚਻⎬ذ‬ї୾⎰ [5] ⏓ , ๲໱ଯఙ‫⎬ث‬ԟ⎘ ⎧ ⎰ [1]⎧ ⎘ ⎿ ⎇ ⎠⏂⎬࠴ਖ␨ ⍬␄ ⎰๲໱ખ (N ⎍ࡒഇ੒⎎ ⎅⎩ ⎎⏓ ‫ ⎅ ⎠ ⎎⎡ ⎠⎅ࢻދ‬.

⎬ ⎰ N ⎬ ԟ⎚⏌ ௥Τߘ⎯ਪ‫ ⏋ ⏉ߧ۽‬ਙ⎐ ࠴ਖ⎚⏌ ) ⎼⎪⎬⎰ਙ⎐ ⎫⎤ ⎧ ⎅⎫ ⎅⏃ ⎯⎯ , ␌⍬⏯ ␑ ച‫ ה‬BN ⏓ ๆ⎅⎠ࣟ‫࠴⎯ܴڬ⎯܄‬ਖ␨ ⍬␄ ⎨⎃ ⏌ 1/N ⏉ ⏋ ⏃ ΄ୢଯ⎬ ‫ۿ‬ਙш⎖ ⏍⎧ ⎅⏌ ⎔ ⎩ ⎍⏐⎌⏌ . ߹‫ݘ‬, ⎃⏌ ३୞ (ຶ⎉ ⎱ 30 ‫ن‬⎯३୞) ⎨Ѿୃॅ⏓ ֆ⏂⏌ ⎔ ⎩ ⏓ ۪⎉ ⏌ ⎩ , ࡚⏀೒⎎ ␌⍬⏯ ␑ ച‫ ה‬W BN ⏓ ๆ⎅⎠ࣟ‫ ⎰⎬܄‬, ␌⍬⏯ ␑ ച‫ ה‬BN ⏓ ๆ⎅⎠ ࣟ‫ ⎬܄‬ಅ⎺⎧ , 1/1025 ⎌⏊ 1/1020 ଟ୞⎯ૠ⎖ ⎯Օஞ⎨‫ ⏁݌‬, ⎚⎫ ⏐⎢ ⎞ ⏍⎡⎒ ‫ۿ‬ਙш⎖ ⏍⎧ ⎅⏌ ⎔ ⎩ ⎩ ⎫ ⏌ .

4⎿⎩⏂ඕ٥֋⎨⎰ , ⏤␧⏰⍬␅ ଯ⎫ ຢӁ‫ ⏌ ⎒⎋⎬ذ‬үࠦച‫ה‬ຟ‫۪ދ‬ഓ٤
[1] S. Das and J. A. Yorke, Super convergence of ergodic averages for quasiperiodic orbits, preprint,
arXiv:1506.06810.
[2] J. B. Greene, Poincaré and the Three Body Problem, Amer. Math. Soc, 1996.
[3] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Willy, 2nd edition, 2008.

⏓ ‫ ⎬⏂ ⎠ ⏌⎚ޔظ‬ๆ⎅⏊ ⏍⏌ ␌⍬⏯ ␑ ⎯⏤␧⏰ ⍬␅ ଑๲[4] S. Das, C. B. Dock, Y. Saiki, M. Salgado-Flores,⎬࠹।⏓ ю⎉ ⎧ , ⎫ ⏂⏊ ⎌⎫ ࡚⏀⏓ ๆ⎅⎠࡚⏀೒⎎ ␌⍬⏯

E. Sander, J. Wu and J. A. Yorke, Measuring
quasiperiodicity, preprint, arXiv:1512.07286.

␑ ച‫ ⏓ה‬஖௡⎘ ⎠ . ⎞⏍⏓ ๆ⎅⏌ ⎩ ࡵ࠵ՃՕஞ⎬ୀ⏐⏌ Ѿୃॅ⎯ஔ଑⎫ ⎪ ⎍‫ۿ‬ਙ⎌⎥‫ۿ‬३୞⎨⎋ ⎔ ⎫ ⎉ ⏌ ⎔ ⎩ ⎍ҷ௨⎖ ⏍⎠ . ⎔ ⎯‫ۿ‬ਙш⎰⏧ ⏦⏵Օஞ⎬਻⎘ ⎧ ⎰߹‫⎫⏍ ⎖ټ‬[5] S. Das, Y. Saiki, E. Sander and J. A. Yorke, Quantitative Quasiperiodicity, preprint, arXiv:1601.06051.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 7スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告

フラグメント分子軌道法と非経験的に最適化した長距離補正
密度汎関数法による電荷移動パラメータの高精度計算
鬼頭(西岡)宏任 a)、安藤耕司 b)a) 名古屋大学大学院理学研究科、b) 京都大学大学院理学研究科

1 INTRODUCTION

ルトニアンを射影することで、ホール伝導バンドを記述するタイトバインディング・ハミルトニア

有機 EL や有機薄膜太陽電池などの有機半導体

ンが得られる。そのハミルトニアンの対角項がサを用いたデバイスの効率は、その有機半導体が持イトエネルギー、非対角項が移動積分になる。

つキャリア移動度に大きく依存する。キャリア移フラグメントの電子状態計算に用いた LC−DFT

動度を理論から見積もるには、分子間電子移動積

法は、交換ポテンシャルが短距離で減衰するとい分やサイトエネルギーなどの主要な電荷移動パラ

う DFT の欠点を克服するために、誤差関数とパラメータを、”分子環境下”における構成有機分子のメータμを用いて以下の(1)式のように電子間の電子状態から決定する必要がある。

クーロン反発項を短距離と長距離に分割し、長距本研究[1]では、フラグメント分子軌道(FMO)法を用いて生体・溶液系の電子移動反応を ab initio
精度で効率的に解析する手法[2]を基に、DNA、液離部分には Hartree-Fock の正確な交換ポテンシャ
ル式を使って計算している。
1 / r12 = (1− erf(µ r12 )) / r12 + erf(µ r12 ) / r12(1)体結晶や分子結晶、導電性ポリマーなどの電荷輸

LC-DFT 法は、電子移動に関与する HOMO と送の理論解析へ応用することを目指した。特に電LUMO の軌道エネルギーが、イオン化ポテンシャ

子相関効果を適切に取り込むために、非経験的にル(IP)と電子親和力(EA)の計算値を再現するよう最適化した(NET)長距離補正(LC)密度汎関数法にμを自己無撞着に最適化することが出来る[5]。

(DFT)を FMO 法のフラグメント計算に採用した。

こうして得られた気相中単分子の HOMO や

最近報告された高精度波動関数理論ベンチマーク

LUMO の軌道エネルギーは、IP や EA の実験値や計算結果[3]と比較することで、我々の手法が現実

高精度波動関数計算結果を良く再現する[5]。従っ的な系の“分子環境下”においても高精度かつて分子集合体に対しては、非経験的に最適化した

著しく低計算コストで電荷移動パラメータを計算

LC-DFT 法(NET-LC-DFT)と FMO 法と組み合わせすることができることを示した。[1]ることで、分子環境下の電子分極効果を取り込んだサイトエネルギー計算が可能になる。

2 METHODS

3 RESULTS & DISCUSSION

FMO 線形結合(FMO-LCMO) 法[4]を用いると、他のフラグメントからの静電場環境下で解いたフ

まずイミダゾールとフランを”face-to-face”にスタラグメント単量体と二量体の Khon-Sham (KS)ハ

ックしたダイマーカチオン間の移動積分に対して、

ミルトニアンから全系の KS ハミルトニアンを構

NET-LC-DFT の精度検証を行った(図 1)。論文

築することができる。p 型有機半導体に注目する[3(a)]から、MRCI+Q と NEVPT2 を用いた高精度

場合、各フラグメント単量体の HOMO を基底分波動関数計算結果を参照値として採用した。

子軌道(MO)とする。その軌道空間に全系 KS ハミ

LC-BLYP が NET-LC-DFT を用いた計算結果であ8 − Academic Center for Computing and Media Studies, Kyoto Universityる。基底関数として重原子には cc-pVTZ を水素原から、前者で無視された電子分極の効果を考慮す

子には cc-pVDZ を用いた。最適化されたパラメー

ることで〜0.1 eVの違いが生じることが分かった。

-1

タμは両方とも 0.33 bohr であり、汎用量子化学
計算プログラム GAUSSIAN のデフォルト設定値
0.47 bohr-1 と大きく違っている。図 1 から、標準的に用いられる B3LYP 汎関数計算と比較して、
我々の LC-BLYP 計算は高精度参照値を定量的に再現していることが分かる。

図 3 アデニン 8 量体の移動積分
図 3 には隣接アデニンサイト間の移動積分の計
図 1(a)イミダゾールと(b)フランのダイマー
カチオン間移動積分の分子間距離依存性

算結果をプロットしている。FMO 法と他の結果を比較すると、特に両端(A1/A2, A7/A8)において移動積分値の顕著な違いが現れていることが分かる。

次に核酸を構成する塩基の一つであるアデニン

環境静電相互作用による軌道緩和(電子分極)は、8 量体を用いて、分子環境効果の重要性を示す。

基底 MO と全系 KS ハミルトニアンに影響し、構造は、理想的な B-DNA 螺旋を取るように作成

FMO 法は両方とも適切に考慮されている。図 2,3-1

した。パラメータμは 0.29 bohr を、基底関数はから、分子環境がこれら２つの因子を通して、ど6-31G(d)を用いた。図 2 にサイトエネルギーの計のようなバランスで電荷移動パラメータに影響を算結果をプロットした。Dimer in Vac は二量体の与えているかを解析することが可能である。

みの計算、Dimer in MM は二量体計算において残りのフラグメントを点電荷で近似したもの

4 ACKNOWLDGEMENTS(QM/MM 計算)である。FOA は、電子分極効果をKS ハミルトニアンに考慮したが、基底 MO には

本研究は、京都大学学術情報メディアセンターのスーパーコンピュータ共同研究制度(若手奨励考慮しなかった結果である。

枠)と科研費新学術領域研究「π空間」(領域番号2007)による。
[1] H. Kitoh-Nishioka and K. Ando, Chem.

Phys. Lett. 621, 96 (2015).
[2] (a) H. Nishioka and K. Ando, J. Chem. Phys.
134, 204109 (2011).; (b) H. Kitoh-Nishioka and
K. Ando, J. Phys. Chem. B, 116, 12933 (2012).
[3] (a) A. Kubas, et al., J. Chem. Phys. 140,
104105 (2014).; (b) A. Kubas, et al., Phys. Chem.
図 2 アデニン 8 量体のサイトエネルギー
環境効果を無視した Dimer in Vac と他の結果を比較すると、0.1〜0.4eV の違いが現れていることがわかる。また、Dimer in MM と(FMO, FOA)の結果

Chem. Phys. 17, 14342 (2015).
[4] (a) S. Tsuneyuki, et al., Chem. Phys. Lett.
476, 104 (2009).; (b) T. Kobori, et al., J. Chem.

Phys. 139, 094113 (2013).
[5] M.E. Foster and B.M. Wong J. Chem.

Theory Comput. 8, 2682 (2012).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 9スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告

射影された時系列データに対する経験的確率微分方程式モデリングとアンサンブル軌道予測
中野直人 ∗∗

1

国立研究開発法人科学技術振興機構さきがけ，北海道大学大学院理学研究院

はじめに
昨今の情報社会はビッグデータ時代を迎え，大量のデータが取得できるようになり，そのデータから(u, v)-平面に射影し，得られた {X tn = (utn , vtn )}10n=1
(tn = 100 + n∆t) を解析する時系列データとする(図 1 右)．
750

対象とする現象の法則を経験的にモデル化する試み
がなされている．データにあらゆる現象の振る舞いが含まれているとみなし，機械学習的なアプローチ
でモデル方程式を構成するのである．一方で，観測できるデータが現象の時間発展をつかさどる変数全

40

50
40
30
20
10
0

30
20

-20

てを網羅しているとは限らず，そのとき経験的モデリングでは観測できる一部分の変数だけでその現象を記述することになる．これは全て自由度を考慮する相空間ではなく，ある部分空間において予測モデルを構築することとなり，その部分空間への射影による情報欠損がそのモデルの記述力に重要な関わり

-10

0

10

20

-20

-10

0

10

20

10
0
-30

-20

-10

0

10

20

30

図 1: (左) Lorenz system の数値解による軌道，(右)平面に射影された Lorenz system の軌道．

3

手法をもつことになる．ここでは，ローレンツ方程式を用いる手法は確率微分方程式（SDE）モデリング

題材に取り，その時系列データをある 2 次元平面にであり，上記の時系列データから，次のように SDE

射影したデータから，どこまで元の時系列軌道を再現できるモデル方程式を構築できるか考察する．

2

データ
ここで用いるデータはローレンツ方程式

x′ = p(y − x), y ′ = −y + x(r − z), z ′ = −bz + xyの平面射影データである．ここで，パラメータは p =10, r = 28, b = 8/3 とする．方程式の数値解は，初期値を (x0 , y0 , z0 ) = (1, 1, 1) とし，時間間隔 ∆t =1.0 × 10−4 で 4 次ルンゲ＝クッタ法で計算する．0 ≤t ≤ 100 までは遷移状態とみなし，アトラクタ上の軌道である t0 := 100 < t ≤ 1100 までの 107 個のデータ√を採取する (図 1 左)．さらにこれを u = (x + y)/ 2,v = z で定まる二次元射影 (x, y, z) �→ (u, v) によって10 − Academic Center for Computing and Media Studies, Kyoto University

dX t = A(X t )dt + S(X t )dW t を経験的に構成する．

平面を長さ 1 の正方形のセルに分割し，各セル x̃m

を通過するデータ点に対して，ある基準時刻 τ 後のデータ点の平均と共分散の変化率を計算する：⟨Xtin +τ − Xtin ⟩x̃m,
τii⟨(Xtn +τ − Xtn )(Xtjn +τ − Xtjn )⟩x̃m
B̃ ij (x̃m ; τ ) =2τ
⟨Xtin +τ − Xtin ⟩2x̃m.
−
2τÃi (x̃m ; τ ) =ただし，⟨·⟩x̃m はセル xm を通過するデータ点に対するアンサンブル平均である．これは，平均的な状態変化度をドリフト項，軌道のばらつきを拡散項として持つような SDE をデータから経験的に構成する手法である．射影による情報欠損はこの拡散項によって補填されることとなる．この Ã = (Ã1 , Ã2 ) と

50

50

40

40

30

30

20

20

10

10

0
-30

-20

-10

0

10

20

30

0

2

1

0

-30

-20

-10

0

10

20

30

-1

(a) Projected orbits(b) SDE simulations

40

40

30

30

20

20

10

10

50

0
-30

-20

-10

50

0

10

20

30

0
-30

-20

-10

(c) Projected orbits(d) SDE simulations

景の緑の軌道は射影データ），(右) 拡散の大きさ trB．

40

40

30

30

B̃ = (B̃ ij ) が離散セル上に構成されたドリフトと拡散係数となる．SDE のノイズ係数の S と拡散係数 B

20

20

10

10

図 2: τ = 0.10 のときの (左) ドリフトベクトル（背との間には B = SST /2 という関係があるため，S を

2B の平方根行列として与える．

50

0
-30

-20

-10

0

10

20

30

0
-30

-20

-10

(f) SDE simulations

40

40

30

30

係数の評価時間を τ = 0.10 として上記の方法で2020構成した Ã と B̃ を図 2 に示す．図 2(左) の緑の軌道

10

10

50

結果

で表された射影されたアトラクタの外縁部では，軌
道の速度が早いためにドリフトベクトルは軌道の接ベクトルとは異なるベクトルとして求められているが，用いる評価時間 τ における平均的な変動成分を表している．また，アトラクタの中央部分では u-成
分は非常に小さく，v-軸が負の方向の成分しか持たない．これは，元の軌道は左右へ交差するデータがあるため，アンサンブル平均で相殺されているからである．すなわち，この部分ではドリフトベクトルだけでは軌道の変動成分を表し得ず，振る舞いが決定論的でないことがわかる．これがデータの 2 次元
射影による情報欠損の顕著な例であると言える．これを補填するのが拡散係数 B̃ であり，軌道交差部においてその大きさ trB̃ の値が大きくなっており，一方でアトラクタ外縁部の決定論的領域ではその値が小さくなっている（図 2 右）．
この経験的に構成された SDE を用いて実際に数値シミュレーションをおこなう．SDE の数値解の計算にはオイラー＝丸山法を用いる（時間刻み幅 0.10）．
このとき，SDE の各係数は離散分割セル上にしか定義されていないため，平面上の任意の点における係数の値は双一次補間を用いて計算する．中心を(u, v) = (20, 25), (20, 40), (0, 20) とするサイズ 1 × 1のセル内に初期値に取った SDE の数値解をそれぞれ
図 3(b), (d), (f) に，同じセルを通過する射影された
データ軌道を図 3(a), (c), (e) に示す．それぞれの軌

0
-30

-20

-10

10

20

30

0

10

20

30

0

10

20

30

50

(e) Projected orbits

4

0

50

0

10

20

30

0
-30

-20

-10

図 3: (a), (c), (e): 1 × 1 のセルを通過するデータの時間 1.0 までの軌道，(b), (d), (f): 同じセル内に初期値を持つ SDE シミュレーションの軌道．その中心はそれぞれ，(a), (b): (20, 25)，(c), (d): (20, 40)，(e), (f): (0, 20) とし，ランダムに選ばれた 20 軌道ず
つをプロット．
道は時間 1.0 の長さでプロットする．図右側はデータから離散的に求められた統計的 SDE の数値解にも関わらず，実際のデータの挙動と定性的に一致す
る．また，平面上での軌道交差は自励系では起こし得ないが，このモデリング手法では適切な確率論的成分を抽出することでそれを達成しており，適切な予測モデルの構築法としての期待が持たれる．この解析手法は射影演算子法との関連もあり，理論的枠組みと統計計算的枠組みをつなぐ良いモデル問題となっており，さらなる発展が期待される．
参考文献[1] N. Nakano, M. Inatsu, S. Kusuoka, Y. Saiki,Time-series analysis and predictability estimates by
empirical SDE modelling, Proceedings of the 47th
ISCIE International Symposium on Stochastic Systems Theory and Its Applications, accepted.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 11スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告

相分離による自己組織化構造形成に与える乱流影響łÒ æ«
ĵĸý§üăÿ1 諸言

 ėõ5FX8] ,á(.^ËçkÅıËĀ®!ĬĐ ® Navier-Stokes ËĀ®

eāŀf_ Às ,ÚÀ/.ÕÉçk!ł

-Navier-Stokes ®c}NC\<SX

ëìěð½!oè--|ı+.úsĺbľ1xľy.

 èçk1ºr+îxøa

+úsĺçk17BKW\:.

ð½.Àsúsĺ. ús/, ÅıËĀ®!Ðķ¨sã+Ąĳĺ

ĺòħ!=JG]AXsģ"/úsĺ+Æ}/Runge-Kutta ã# Crank-Nicholson²À/.ė©Ďē}Úī Ĩ²ð1Ã

ã1ôÎĳĂs/FBD[]9Úī. ,ÝĖ»ÕÉ¢ħç/!^ÄÂ/dç1~

lĞKZ?=¹ô/.ė©Ďē}Úī

Č}^ÛĆËdç-öÑö

ï¸1â.ğ!}øĎÀ*îxÖg ĥčĦ%,/. úsĺ÷óÎ

Ög1įôąËĿ1 2563

Ùísw

dç!dçð½1đÃ.(°uçkxø³Ľ!ÓģÍİs 

x1ånªð½ C4V]O49Z=;

Ôüă!
úsĺ+.çkxø³Ľ

]X Y4GX>Ç! 21 .6]Aïdç ³ĽeÀsċçk1¢]HVQ]@!tÑÖgM[4DG4>1

ħdç<PTY]<U\1ÌĦ%

yúsĺ

µĭĜ.Ög

Cahn-Hilliard ® 

2 解析手法&/.HVQ]@!

3-methylpentane  nitroethane 1ô Pineet al.[1]Ł¢¹.m1ôģ×!

50%-50%.ĻĒ» eCahn-Hilliard ® Navier-Stokes ËĀ®1ĬÀÀsèçk(ĎÀ AB)¢Às ĨĲ

.Ög(SPTU)fCahn-Hilliard ® '1tÑÀsàñ
¬ρAρB,f` 6]A]HVQ]@1Ĕ

.

ψ=

ģ (SP) Navier-Stokes ËĀ® '1ģ
(TU))àĩ (Ĝ

ρ A − ρB
ρ0

3 結果と考察ρ0 !èð½ « ¬. 6

úsĺ+.ė©Ďē}Úī ÷ò1Ġ.(]A]HVQ]@¢úsĺ ÅıËĀ®!

6]A]HVQ]@ h¼Êļ .Îĳ

Cahn-Hilliard ËĀ®

÷¥1 Fig. 1 þ&çk1ĕ¿Ö

∂2 µ
∂ψ
∂ψ
+ uj
=Γ 2
∂t
∂x j
∂x j

g SP

ďØåù.Èíð úsĺEQ

4\ ò/Ċ}(FBD[]9}).pø=JG]AXsģ Çmøqò ..t !Îvxi !­Üui !Ī¬

 0 ./¢çk ³Ľ1ĕ¿

Γ!RIWC3HVQ]@µ!}NC\<SX.(Ög SPTU)sĺĿ!ò/. dç.}NC\<SX! Landau-Ginzburgè ³Ľ1 ¯ i"/ĚÚī n

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

Phys. Rev. A, 29 (1984), pp.308-313.

-Ď2-dç ³Ľ úsĺ°
ò/. 0 .Þúsĺ ÷ó

2.

Horiuti,K.

andTakagi,Y.

.+.dç$ ³Ľ1Ħ%.(ú“Identification method for vortex sheetsĺ ĭĜÎv .çk 5FX8]

structures in turbulent flows”, Phys.

=L9DX1 Fig. 2 þàĩ (

Fluids, 17 (2005), 121703.

dç ' Ög TU

ďØ)Îþ.Čĉçk ' 5FX8]=L
9DX!äÇ k~101-102 ¾»£Ŀ=;]X ,£=;]X$ 5FX8]7=;]E Ġ,/. úsĺ1ĬÀ=;]X(k<101) ,5FX8] ę é¤dç 5FX8]7=;]E úsĺ ÷ó+Ķ/. 0 .
dçc Úī!àĩø=;]X  <]Dð ê¦Ï£=;]Xú±.êĈ 

t = 0.1

t = 5.0

t = 10.0(b) SPTU û,/.[2]  5FX8]=L9D
X é¤!úsĺ zØ ·čêĈÚī!
 =;]X ê¦¢ò/.
1þ. ďØ! Fig. 1

ġ}ďØ

¦ð ¯ i"/ė©Ďē}ÚīĢ¡/)^Ę.

t = 0.1

4 結言

t = 5.0

t = 10.0

Fig. 1 Coarsening process of phase separation: (a)with only phase separation, (b) turbulentÔüă!úsĺ(=JG]AXsģ)1ĝò

flow with phase separation.

. Cahn-Hilliard ËĀ®1ĬÀdçĤć1eÀsċ^ÛĆËdç1¢ħĜúsĺ
+.ė©Ďē}Úī ²Àdç ĴĬ»Ħ%ÔüăôĤćÖg!ús
ĺ Ċ}ĮĀ dç +Ķ/ė©Ďē}Úī <]Dð dçÚī(ê¦)ŀj ²½. 0 &dç 5
FX8]7=;]E Ģí ,úsĺ ³Ľ1ĕ
.àĩø =;]X5FX8]7=;]E tÑßĹ úsĺ+Áu/- úsĺ ³Ľ.=;]XĿ1
Ė{øu¶.+ė©Ďē}Úī1
çku¶ Ģí ,}. Ñ´ .
Fig. 2 Effect of phase separation on energy

参考文献
1.

spectra.

Pine, D. J., Eswar, N., Maher, J. V. andGoldburg,

W.

I.,“Turbulent

suppression of spinodal decomposition”,

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

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

液柱内の温度差と濃度差に起因するマランゴニ対流の不安定性の解明
水口 尚
琉球大学 工学部 機械システム工学科

1 緒言で、Si/Ge 結晶作製を対象とし、温度差と濃度差両者に起因するマランゴニ対流の共存効果について

半導体材料であるシリコンとゲルマニウムの合一連の研究を行ってきた(1)。

金（Si/Ge）は従来のシリコン製チップに比べ大幅な

本研究では温度差マランゴニ対流と濃度差マラ

消費電力の低減が可能なため、エレクトロニクス製

ンゴニ対流の二つの対流が共存した場合の対流構品の更なる小型化、安定化を可能にする。しかし、造に及ぼす各対流の大きさの影響について調査しシリコンとゲルマニウムの比重差が極めて大きいた。

ため通常の条件で融液から結晶を作製すると重力偏析が生じ、均一でかつ大型の結晶を作製すること

2 解析方法が困難になる。そこで重力偏析の生じない微小重力環境においてフローティングゾーン（FZ）法を用い解析には FZ 法の半分を模擬したハーフゾーンモて Si/Ge バルク単結晶を作製することが提案されて

デル（Fig.2(a)）を用いた。融液の物性値には Si/Ge

いる。FZ 法の概略図を Fig.1 に示す。

（Pr 数：6.37×10-3, Sc 数：14.0）をアスペクト比（As = L/a）には 0.5 を用いた。また、無重力を想定するとともに、気液および固液界面形状は平坦とした。解析基礎式には連続式、Navier-Stokes 式、エネルギー方程式、拡散方程式を用い、これらの式
を有限体積法によって離散化し PISO 法を用いて解析した。温度の境界条件は、上面を最低温度 Tc、下面を最高温度 Th で一定とし、側面の自由表面ではFig.1 Schematic diagram of floating zone method.

断熱とした。また、濃度の境界条件は、模擬的に温度の境界条件と同様にした。速度の境界条件は、壁FZ 法は、原料となる多結晶を上部から供給し、多面上で no slip 条件とし、自由表面上では温度差と結晶と単結晶の間に融液を作り、その全体を下方に濃度差による表面張力差に起因するマランゴニ対移動させ融液を冷却して単結晶を得る方法である。

流を考慮した。対流の強度を表すマランゴニ数を次融液の自由表面に温度差が大きいので不安定なマに示す。

ランゴニ対流が発生し、均一な結晶の成長を阻害する。そこで FZ 内マランゴニ対流に関して多くの研Ma T   T TL ,


MaC 

 C CL

究がなされているが、そのほとんどが温度差に起因対流の向きは、温度差マランゴニ対流では温度がするマランゴニ対流を対象にしたものである。しか低い方が表面張力は大きいため、自由表面上を上昇

し、Si/Ge の結晶作製を対象とすると、結晶化の際する流れ、Si が Ge に比べ表面張力が大きいため、に偏析係数が大きいため、融液内に濃度差が生じ、濃度差マランゴニ対流も自由表面上を上昇する流Si と Ge には表面張力に差があるので、濃度差に起れになっており、両対流とも同方向に流れている。

因するマランゴニ対流も考慮する必要がある。そこ

初期条件として、速度なし、最低温度 Tc、濃度 CSi

14 − Academic Center for Computing and Media Studies, Kyoto University（Si:100%, Ge:0%）で一様とした。解析に用いた格

定常流となる。これは、濃度差マランゴニ対流に対子(Fig.2 (b))は、側面の自由表面と上下の固液界面付して温度差マランゴニ対流が支配的になったため

近の格子が細かくなるような 14.4 万の格子を用いと考えられる。

た。これらの計算にはフリーの熱流体解析ソフトである OpenFOAM を使用した。

Fig.2 Schematics (a) and computational grid (b) of
Fig.3 Concentration contours at 1000s ((a): MaT =the liquid bridge.

2857, MaC = 0, (b): MaT = 0, MaC = 1786, (c): MaT =3 結果と考察714, MaC = 714, (d): MaT = 1071, MaC = 1786, (e):MaT = 2143, MaC = 1786, (f): MaT = 2857, MaC =Fig.3 には(a)に温度差マランゴニ対流のみ、(b)に1786)

濃度差マランゴニ対流のみ、(c),(d),(e),(f)に温度差と濃度差のマランゴニ対流が共存した場合の融液内の Si 濃度分布を示す。垂直断面図では、自由表面付
近にマランゴニ対流による大きな渦ができており、その渦の大きさが周方向に変化し、水平断面図のよ
うに、周方向に規則的な模様の濃度分布になっている。(a),(c)は、定常流（(a)は非軸対称定常流）、(b),(d),(e),(f)は、m = 4,5,6,7（m:周方向の対称面数）の振動流となった。これらを温度差と濃度差のマランゴニ数でまとめた対流構造のマップを Fig.4 示す。
縦軸上の温度差マランゴニ対流のみの場合は軸対
称定常流から m = 3 の非軸対称定常流に遷移し、横軸上の濃度差マランゴニ対流のみの場合は軸対称
定常流から m = 4 の振動流へと遷移した。このマッ

Fig.4 Flow mode map.

4 結言

プより、両対流が共存した場合は、温度差・濃度差

温度差と濃度差に起因するマランゴニ対流の共マランゴニ数を大きくしていくと波数（m）が大き

存効果について数値解析を行ったところ、下記の知くなった。濃度差マランゴニ対流だけでは、m = 4

見を得た。

の振動流だけであったが、温度差マランゴニ対流も

共存効果により波数（m）が変化する。

共存することで波数の異なる振動流が得られた。



MaC の臨界値付近では MaT を大きくすると

MaC = 714 以下では MaT ＝ 2857 まで大きくして一度振動流になり、さらに MaT を大きくす

も波数は大きくならず、定常流となることから、濃ると定常流となる。

度差マランゴニ数がある程度大きくないと振動流にならないと考えられる。濃度差マランゴニ数の臨

引用文献

界値付近（MaC = 893,1071）では MaT を大きくす(1) H.Minakuchi et al., J. Crystal Growth, 385

ると一度振動流になり、さらに MaT を大きくすると(2014) 61.

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

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

多孔質壁面による乱流構造の変化と抵抗低減効果
S½ æþ
Īĭëéðí1

緒言ĆąĄĉ İ,ÒÐâÝ Wç īÌ

3

結果と考察

ę ),UÑª© ÔnÂ.ē_,

ª©[ÔnÂ.]¬ĆąãÊĀ¤¾¶'ª©[Ôà(DR).Ù ĳ¿X(Cf0)

D?M7K¶ è*-,D?M7K),ę ĳ¿X(Cf p) ª©`µ Í *ó,ª©[Ô!7K L2M9PÞ¤òÈi

ċ| ¯Å-, ĔöF3@;E ď⎛ Cp ⎞
DR[%] = ⎜⎜1− f0 ⎟⎟ ×100
⎝ Cf ⎠¸!ă *ºĥk¡¨Ï êñ½éð!D?M7Kā ¤ę.
«irõÇ§·î)G=Kp,

ª©`µ Cf !fåĤĉÀĩ<HA

µbĐó)+D?M7Kā ª©

KUÑ , FIK ¥ò)+ói[2]v¿[ÔF3@;E.¸* D?M7K Ñ\X , CfDR Ù<HAK , ĳĜģ¤),UÑª©$ Ĵ.ėĖ,

±²ĝ)ÜËdp-ę s
δp+. Table 1 ìę\ s& Ñ\ĵ

2

計算手法

3% 5%  ª©[ÔnÂ.ì ę\ s!<HAKUÑ äďÁĘ!ĉÀĩ<HAKUÑ,

s& y+~10 wî  , ěc

D?M7K!qjÓÎ,ÃÉÇĞgÎ,TÑ·x ĝo(u’) učp.Fig. 2

e+Ě&ĂÚ~)gĦ Î.­iìFigure 2(a)(c)x ·xY#ûß ĵ!tü,½éð! ¤ę.ę\ĹĝO[ĝ:>JP5z"-,Figure 2(b)(d)ÓĜÖ i áĘý °ĸÄ Z.

ù(ĺOÛć)!UÑÕ1AK4P´ġ òÿ¦ Darcy l)+G=Kp[1].ĉbù.--Ċ,ę\Q!:ĻFig. 1ļ}\ Wě!D?M7Kā s>JP5ÇĞ ę\gĦ * Ö i)&.ÿ¦ę\ ,³ħ·î1AK4P´ġ aĠ-, 

! Navier-Stokes ğ÷, ĳä¿

3%D

X! no-slip ¿XTÑ·x:CN·x!

ø¬-, j ,%[ĝ:>JP5

ĝ~m-(y¼ä¿X.ãÇĞ!

³ħ·î!»Ĭ\ïÏ)+ı´p¹

HAKS¢$ Q*-,§-,

ĩå! 1 Ëô20IPĮďÏ.ãı'[ĝ:>JP5Wě!ÕÇĞ´ġ´p-·î! PISO /K8J;E)

Ď-,Flat !´ġĵ T[ĝ:>Jď<HAKgĦ.Ñ-,ÑĨ× L0B

P5RĦj, 3%D ![ĝK;µ!5600 ,:>JP5 QĦ´ġĵ Č*-,
-!

6P:!:>JP5ÇĞ

FlatwòĴ)ÔĝÑ\ h#<ę\ äg ĲĢÑ-®+gĦ * Ö i  ' ě

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

c[ĝ Ñ\ Ø%+ [ĝ Ñ\
}\ 
4

(a)

ĳ±².ú{,ÿ*-,

結言
Ñ\Ĝģ¤.»,D?M7K!D?M7

Kā UÑäwî s&.( Wě UÑÇĞ pķĈª©[ÔnÂ(b).ì j ½ďÁì-)
ª©[ÔnÂ º,)ās& ,'L0BK;µ^¤.ÆĒį ÊĀ¤¾¶$ £ã.Æđ, V

ĕĶ

,

(c)

参考文献
1.

Suga, K., Mori, M., Kaneda, M., Int. J. Heat Fluid
Flow, Vol. 32, pp. 586-595 (2011).

2.

Fukagata, K., Iwamoto, K. and Kasagi, N., Phys.
Fluids, Vol. 14, pp.73-76 (2002).

(d)Fig. 2 Snapshots of the flow field around streak structurenear the lower wall for the cases of Flat (a,b) and
3%D (c,d) in the x-z plane of y+~20 (a,c) and theFig. 1 Schematic diagram of a wall coated with hydrogel
painting.

y-z plane of x+~1500 (b,d). Black lines represent
the contour of the second invariant of velocity

Table 1 Drag reduction ratio and friction coefficient.
Case

Cf x 103

Flat

7.25

-

-

1%D

7.31

2

- 0.2

3%D

6.47

5

10.6

5%D

6.55

9

9.6

10%D

7.16

18

1.3

20%D

7.56

37

- 4.3

δp+

DR[%]

gradient tensor, Q+ = 0.2 and dark gray lines thedissipation rate, ε+ = 2.8 x 10-3.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 17スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告

有機太陽電池応用を目指した新規光機能性有機分子材料の構造と電子構造の解明²¯~aùāËÐÍnª

1 緒言ċ 2,6-I=(59@V57<)J3EVÚ ¨ ZnPBAT-o-C8 /èçČFigure 1čßÕ¢Ăą·ČDSSCčÄ+WDA9=>V Ý ć/_+[ÞÇ

2 結果と考察S2ÕWDA9=4FV8Z½eP
0.7-0.8V Ï=OA@(*ý«ąČVOCč

pĀ,#   ĉ
Grätzel )S2ÕWDA9=(*&õoúv÷jąe/&;HVCWDA9=ºÈ
YD2-o-C8 Ö$z.+ VOC = 0.965 V¨tÁ η = 11.9%/ó¤ČFigure 1č1;H
VCWDA9=/Ä+ z[ÞÅĆïąÎu ð +%YD2-o-C8 (ÒgÇċÚ¨/lßÕ/Ä+ä + °ËÐcq:VZ
LýÆ ZnPBAT2-4 R<@V/ÒgÇ

O
C8H17

N
N

Zn

N

N

Zn

N

COOH

ZnPBATN

ZnPBAT {´Ê¹üü¹üinâ

)KXZD*450-500 nm bîk¦ĄÝ©,+#Q & YD2-o-C8 µ!ü¹ü<JCü¹üĈk¦ĄÝ&|\+  .

N

N

C6H13

N

ZnPBAT-o-C8 }y=M9CV Soret 

2.0

COOH

N

V¾¼`}y=M9CV/ Figure 2 ÌN

ZnPBAT-o-C8N

ZnPBAT-o-C8ZnPBATYD2-o-C8 4?GZ

ε / 105 M–1cm–1

C8H17
O

2.1 新規ポルフィリン色素の諸物性と色
素増感太陽電池への応用

C8H17
O
O
C8H17

1.5
1.0
0.5

C6H13
N
N
N

Zn

N
N

0

C8H17
O
O
C8H17

Figure

400

2.

500

600
700
Wavelength / nm

UV–vis–NIR

absorption

800

spectraof

ZnPBAT-o-C8 (solid line), ZnPBAT (dashed line) and

YD2-o-C8 (dotted line) in EtOH.
COOH
YD2-o-C8

Figure 1. Structures of porphyrin compounds.

18 − Academic Center for Computing and Media Studies, Kyoto University;HVCąæ¼/Ä>Vhã/à-ZnPBAT-o-C8 ÊØą» JSC = 13.2

mA cm–2VOC = 0.83 V­Ù ff = 0.75, η = 8.2%

 YD2-o-C8 >V ÝôKXA7Y:

×±t±ô(*úv@?Yd)ąæ ) ,   ZnPBAT-o-C8 

YD2-o-C8 µíJSCVOC mf]+
2,6-I=Č59@V57<čJ3EV¬

YD2-o-C8   2    +      êúvg ïąÎu ë *' Ã
%§,+ (°çÓ×±Ċ×±ÂéÇæû[r¿¡Û

ZnPBAT-o-C8  1 %KXA7Y

+

:t± ^wnúv@?Yd)ąæ
êúvg ïąÎu ë *' >V

3 参考文献Ý f]Ü),+
[1] A.Yella, H.-W. Lee, H. N. Tsao, C. Yi, A. K.

Chandiran, M. Nazeeruddin, E. W.-G. Diau, C.-Y. Yeh,

2.2 理論計算
NVJ1UY®v³ñ(JXY
B10ìòą³ñÉå/+%¶ÿ¬¸ČDFTč(+ÂéçÓČB3LYPč/àÿ¬Ô 6-31G*/ö¥Gaussian

09 LX:TQ/Ä
Figure 3  ZnPBAT-o-C8

ZnPBAT ( YD2-o-C8 ®ċáxìò
ČHOMOč®fÑìòČLUMOč +ąn/Ì[Þ0Y6Zøe LUMOą  "sëÀ£ßÕúv
@?Y 3d ìòþą6ALUY: 
+  É),+ZnPBAT-o-C8 
ZnPBAT ' YD2-o-C8 {´0Y6Zøe
LUMO ą +  .*

ZnPBAT-o-C8  ą  º l t Á  ZnPBAT 'YD2-o-C8 {ÏċÜ),+%;HVCąæ¼/Äă ZnPBAT-o-C8

S.

M

Zakeeruddin“Porphyrin-Sensitized

Solar

and

Cells

M.

with

Grätzel,Cobalt(II/III)–Based Redox Electrolyte Exceed 12 Percent
Efficiency”, Science, 334, 629-634 (2011).

[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] H. Hayashi, A. S. Touchy, Y. Kinjo, K. Kurotobi, Y.
Toude, S. Ito, H. Saarenpäa, N. V. Tkachenko, H.
Lemmetyinen,

and“Triarylamine-Substituted

H.

Imidazole-

Imahori,andQuinoxaline-Fused Push-Pull Porphyrin for High

Performance

Dye-Sensitized

ChemSusChem, 6, 508-517 (2013).

SolarCell”,[4] T. Higashino, H. Imahori, “Porphyrins as Excellent
Dyes

for

Dye-Sensitized

Solar

Cells:

RecentDevelopments and Insights”, Dalton Trabs., 44,448-463 (2015).

Figure 3. Some sets of molecular orbital diagrams for(a) ZnPBAT-o-C8 (left), (b) ZnPBAT (middle), and (c)YD2-o-C8 (right) obtained by DFT calculations with

B3LYP/6-31G(d). To simplify the computations, alkyl
chains on the diarylamino groups were replaced withmethyl one, whereas octyloxy groups on the phenylgroup were replaced with methoxy ones.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 19

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

急縮小・急拡大流路における高分子溶融体流れのマルチスケールシミュレーション
原田 浩平 ∗ , 谷口 貴志 ∗∗1

京都大学大学院工学研究科 化学工学専攻

はじめにり現実的である多分散高分子流体流れにおいて MSS

高分子製品の製造では，成形加工の段階で高分子法が適用可能であるか検討を行った．

溶融体（高分子流体）の流動を制御する必要がある．
このため，流体力学方程式 (式 (1)) を用いた流動シ
ミュレーションが行われてきた．
ρ{∂v+ (v · ∇) v∂t}

= ∇ · σ − ∇P + F(1)

ここで，v ，ρ，P ，σ ，F はそれぞれ，速度，密度，
圧力，応力，外力場である．高分子流体はマクロス

計算手法の概略

2

ミクロ・マクロ各スケールの概要とそれらを組み合わせた MSS 法について以下で述べる．

2.1

マクロスケールの計算手法の概要

ケールの流動とミクロスケールの状態（高分子鎖同計算手法として粒子法 [3] を採用した．粒子法は士の絡み合いや配向状態）が強く相関するため非常

流体を流体粒子と呼ばれる要素に離散化して解く方に複雑な振る舞いを示すことが知られている．この法であり, ラグランジュ的描像に基づいている．応ような流動特性は，式 (1) 中の応力 σ に集約される．

力は高分子鎖の微視的な状態から決定され, その値

従来は，この応力 σ を速度勾配と少数のパラメー

は高分子鎖の過去の変形履歴に依存する．そのため，

タを用いた構成方程式から求め，シミュレーションを過去の履歴を自然に取り入れることのできるラグラ

行ってきた．しかし，構成方程式を用いた方法では，ンジュ的手法を用いた．

用いた構成方程式そのものが必ずしも対象としてい2.2

る流体の記述に対して常に正しいとはいえない．まミクロスケールの計算手法の概要再現しようとも，ミクロスケールの状態とマクロスモデルとして，Slip-link モデル [4] を採用した．
Slip-link モデルは，絡み合いによる束縛のために各ケールの流動を直接関係づけることはできないし，製高分子鎖がチューブ状の領域を運動するというレプ品の品質や特性に大きな影響を与えるミクロスケー

テーション理論に基づくモデルである．Slip-link モ

た，たとえ広い流動領域で正しくレオロジー特性をルの情報を得ることはできない．そこで，このようデルでは，各高分子鎖の絡み合い点 (Slip-link) とそな欠点を解消するためにミクロモデルとマクロモデれを結んだ線 (Primitive path) を考え，それらの運ルを相互に組み合わせて解くマルチスケールシミュ

動を求める．また，絡み合いをバーチャルな情報としレーション (MSS) 法が提案されている．MSS 法は，て取り扱うため (図 1)，計算負荷を大幅に低減するこ

平行平板間や円柱周りといった比較的簡単な流路形とができる．具体的には以下の手順で計算を行った．

状における単分散高分子流体流れについて成果を挙げている [1, 2]．しかし，工業的な問題に MSS 法を(1) 流れ場によるアフィン変形

応用するためには，より複雑な系での検討が必要で(2) 式 (2) に基づく Primitive path 長の変化ある．そこで本研究では，工業的に用いられる典型(3) 高分子のレプテーション運動

的流路形状として急縮小・急拡大流路を採用し，よ(4) Slip-link の開放・生成

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

(5) 式 (3) による応力算出

( )
dL
dL
1
(L (t) − Leq ) ++ g (t)=−dt
τR
dt aﬃne⟩⟨rα rβ, α, β ∈ {x, y, z}
σαβ = σea |r|空間単位 b : 流体粒子の初期間隔
時間単位 τe : ミクロスケールの単位時間

(2)
(3)

応力単位 η0 /τe ，ここで，η0 :ゼロせん断粘度
ミクロスケールの無次元化では，

ここで，L，Leq はそれぞれ Primitive path の全長，

空間単位 a : 平衡時の Slip-link 間の平均距離

バルク平衡時の全長，a は後述する空間単位である．

時間単位 τe : 長さ a の Primitive path のτR は長さ Leq の Primitive path の Rouse 緩和時間，
g(t) はランダムな熱揺動力による L の変動，r は隣
接する Slip-link 間の結合ベクトル，σe は後述する応力単位である．また，< (· · ·) > は全 Primitive path

での (· · ·) の統計平均を表す．

2.3

Rouse 緩和時間
応力単位 σe : Slip-link モデルの単位応力とした．ここで，Slip-link モデルは十分粗視化されているため，マクロスケールの時間単位はミクロスケールと同じ時間単位を用いた．以下に示す結果は，上記の値を用いて無次元化されており，無次元化さMSS 法の概要 (図 1)

れた変数は*を付けて示す．

本研究では，マクロスケールの各流体粒子中に，

Slip-link モデルを用いた高分子シミュレータを埋め込むことで MSS 法を実行した．以下の手順で計算を
行った．

4

結果と考察
急縮小・急拡大流路での流動に対して MSS 法を適用

(1) マクロスケールの速度勾配場を計算(2) 求めた速度勾配場の下，各流体粒子中で Sliplink モデルを用いたシミュレーションを行いマクロの応力場を計算(3) その応力場を用いて粒子法で式 (1) を解くした．ここで，拡大部高さ，縮小部高さはそれぞれ，

H1 ∗ = 48，H2 ∗ = 12 とし，圧力勾配 (dP ∗ /dx∗ =−5.0 × 10−4 ) を印加した．また，高分子流体は平衡

時絡み合い数 Zeq = 8.74 の高分子鎖からなる溶融体に，Zeq = 65.7 の高分子鎖を 7.3w% 添加した多分散MSS 法では応力場を計算する際，ミクロスケールのの系である．また，統計精度向上のため各流体粒子

シミュレーションを行う必要があり莫大な量の計算中に高分子鎖数を p = 1000 と設定した．

が必要となる．しかし，本研究で用いた MSS 法は各

図 2 に t∗ = 250 での (i) ひずみ速度テンソルの xy

流体粒子でミクロスケールの計算が独立であるため，∗成分の絶対値 (|Dxy|) と (ii) 応力テンソルの xy 成分

並列計算を行うことにより計算時間を大幅に低減す∗|) を示す．またミクロスケールの情報の絶対値 (|σxyとして，各流体粒子中の Primitive path を，重心を原ることができるという特徴がある．

点として重ね合わせた図を図 3 に示す．ここでは見易さのため Zeq = 8.74 の Primitive path のみを表示している．この図から大まかな高分子鎖の配向の様子がわかる．ここで，図 3 の (a)～(c) は図 2 (i) 中の (a)～(c) とそれぞれ対応している．図 4 は，(i)Zeq = 8.74の高分子鎖の絡み合い数 Zshort ，(ii)Zeq = 65.7 の高分子鎖の絡み合い数 Zlong の空間分布である．中心部

図 1: Multi-scale simulation method using the par-ticle method and the slip-link model.

(y ∗ = 0 付近) の流れについて，縮小部 (a) では x 方向の流速が増大するため高分子鎖は x 軸方向に伸張し配向した．拡大部 (b) では y 軸方向に広がる流れが発生するため，高分子は y 軸方向に配向した．図 4(i) より，Zshort は，x∗ > 6 の領域で平衡時絡み合い数よりも低下していることがわかる．これは，縮小部

3

無次元化

マクロスケールの無次元化では，の強いせん断流動によって引き伸ばされた Primitivepath(図 3 (c)) がこの領域では緩和するため，末端のSlip-link が開放され生じたと考えられる．これは高京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 21

分子の履歴効果であり，適切に取り扱うことができ

本研究で用いた手法は，ミクロスケールの計算をているといえる．一方，Zlong (図 4 (ii)) ではそのよ行うため計算量が膨大となる．そのため，本共同研

うな低下は起きていないが，これは，Zeq = 65.7 の究制度で並列計算機を利用することで並列計算のプ

高分子鎖は Zeq = 8.74 の高分子鎖に比べ Rouse 緩和

ログラム開発を滞りなく進めることができた．ここ

時間が非常に長いためである．

に感謝の意を表したい．

参考文献∗
図 2: Distributions of (i) strain rate |Dxy| and (ii)∗| at t∗ = 250.
stress |σxy

図 3: Configurations of primitive paths.

図 4: Distributions of entanglement number of (i)short polymers Zshort and (ii) long polymers Zlong
at t∗ = 250.

5

結言
MSS 法を急縮小・急拡大流路内の多分散高分子流体流れに適用し，マクロな流動とミクロな情報の関係について考察した．MSS 法は履歴効果を取り扱うことができ，複雑な系であっても有効であることを明らかにした．今後の方針としては，温度分布の導入や枝分かれのある高分子鎖への応用が考えられる．
将来的には，MSS 法が工業的な問題に応用されることが期待される．

< 謝辞 >

22 − Academic Center for Computing and Media Studies, Kyoto University[1] T. Murashima and T. Taniguchi, Europhysics Letters., 96, 18002 (2011).
[2] T. Murashima and T. Taniguchi, ”Multiscale simulation of history-dependent flow in entangled
polymer melts”, Europhysics Letters, 96, 18002(2011).
[3] J. J. Monaghan, Annu. Rev. Astron. Astrophys.,
30, 543 (1992).
[4] M. Doi and J. Takimoto, Phil. Trans. R. Soc.
Lond. A, 361, 641 (2003).

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

周囲気体を考慮した HZ 液柱内温度差マランゴニ対流における

粒子集合現象の数値解析

中村 翔一
東京理科大学大学院理工学研究科機械工学専攻

1 緒言アプローチが必要である．また，液柱だけでなく周囲気体も着目しているため，計算コストが非常

自由表面上に温度差を与えると，それに起因しに高くなりスーパーコンピュータの利用が不可欠た表面張力差を駆動力とする対流が発生する．こである．本研究では，実験に近い条件で従来よりの流れは，マランゴニ対流と呼ばれており，特にも広範囲な領域での数値解析を実施し，周囲気体

微小重力下やマイクロスケールにおいて支配的にが HZ 液柱内流れの不安定性に与える影響を調べなることが知られている．

ることを目的としている．

Smith と Davis [1]は，薄液膜において自由表面に対して水平方向の温度勾配が発生している場合，
マランゴニ対流により hydrothermal wave(以下2 数値解析手法

HTW)が生じる流れの不安定性を線形安定性解析

数値解析で用いた計算対象を図１に示す．これ用いて発見した．現在，図 1 のような Half-Zone

は，液柱(Pr=4)の周りにシールドを取り付けその液柱(以下 HZ 液柱)はマランゴニ対流の不安定性中にアルゴンガス(Pr=0.67)を充満させている．支を研究するための代表的な系となっている．HZ

配方程式は円筒座標系非圧縮非定常の３次元ナビ

液柱においてマランゴニ対流の強さを示す無次元

エストークス方程式，エネルギー方程式，連続の数としてマランゴニ数(以下 Ma)がよく用いられ式であり不等分割スタッガード格子を用いて解いている．Kamotani ら[2]は HZ 液柱において高プた．また，有限差分法を用い 2 次精度中心差分にラントル数流体(Pr=24～49)を用いて実験的に自より離散化した．連続の式，ナビエストークス方由表面での熱輸送を示すパラメータであるビオ数

程式のカップリングにはフラクショナルステップ(以下 Bi)が臨界値に与える影響を示した．2014

法，ポアソン方程式の解法として Chebyshev

年から実施されている日米共同宇宙実験では，

SOR 法を用い，打ち切り誤差は 10-6 とした．時間

Kamotani ら[2]の実験を元に Bi に着目している．

進行方法は，液柱には 2 次精度のさらに，
2018 年から実施予定の日欧共同宇宙におAdamus-Bashforth 法，周囲気体は，軸方向の粘いても，Bi と可視化用トレーサ粒子が一本の線内性項に 2 次精度の Cranck-Nicolson 法その他の項に集合する現象である粒子集合現象に焦点を当て

は Euler 法 を 用 い た ． 粒 子 の 計 算 方 法 はた研究が行われる予定である．しかし，実験装置

Hoffmann と Kuhlmann [3]の手法を用いている．

の制約上得られない情報も多く，数値解析によるこの計算コードは，Open MP と MPI を用いたハ

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

イブリット並列化を行った．

場合さらに計算負荷が増加するため，今後はさら3 結果と考察に並列化の効率を上げる必要がある．

今までの先行研究では、Bi を考慮し液柱内流れの不安定性を調べる際計算コストが増加するため、周囲気体を無視し簡略化した境界条件を定義し，周囲気体の数値解析はほとんど行われこなかった．
周囲気体を考慮した液柱と周囲気体を無視し自由表面に Bi に関する境界条件を定義した液柱のBi 数と臨界値の関係を図 2 に示す．まず，二次元
図 1 計算対象

定常流から三次元振動流へ遷移する臨界マランゴ
ニ数 Macr の算出方法は Imaishi ら [4]の方法を用8000いた．局所的なビオ数を軸方向に平均したビオ数7000の条件においても低 Biave(≦0.3)では Macr が減少
し， Biave≧0.3 では Macr が増加している．このMacr

Biave は Macr 上での値である．これを見ると，どwith amb gas
without amb gas(T amb=Tc)
without amb gas(T amb=Tc+Tz)

6000
5000
4000

傾向は，Simic-Stefani ら[5]の実験結果と類似している．彼らは，Γ(=H/R H 液柱高さ R 液柱半

30000径)=0.55-1.28 の範囲で蒸発による自由表面でのBi(=3.2,12.7)を比較し，Bi を約 4 倍増加させたこ

1

2

Biave

図 2 Macr と Biave の関係とで Macr が約 2.5 倍増加したことを示している．
また，彼らの実験結果では Γ=1.05,Bi=3.2,Macrが約 6500-7500 となっており計算結果 Γ=1.0,
Biave=2.1,Macr=6272 と近い値になっていることから，今回の数値計算の妥当性が高いと言える．
さらに，周囲気体を考慮した場合 Biave=0.8 で従図 3 粒子集合現象

来の境界条件と比べMacr が約1.5倍変化している: 左 実験結果[6] 右 数値解析結果

ことからこの計算が必要不可欠であると言える．
図 3 は粒子集合現象の実験[6]及び数値解析結

参考文献

果を示している．
計算条件は Γ=0．
66，
Ma=7200，
無次元粒子径 dp=0.028 であり，周囲気体を考慮[1]132, 119-144 (1983)

していない液柱内に粒子を混入し計算している．
実験結果と数値解析結果を比べると粒子集合現象

[2]

4 結言

Y. Kamotani_et_al, Int. J. Heat Transfer,46, 3211-3220 (2003).

の形が類似していることから計算の妥当性が高いと考えられる．

M.K.Smith and S.H.Davis, J Fluid Mech,[3]

E.Hofmann, H.Kuhlmann, Phys. Fluids,
23, 072106 (2011)[4]

N.Imaishi_et_al, Crystal Growth 230,
164-171 (2001)

より実験に近い条件である周囲気体を考慮した[5]

HZ 液柱の場合，従来の簡易的な境界条件と比べ，Macr が最大で 1.5 倍変化していることからこの計算手法が必要不可欠であると言える．また，粒子

S.Simic-Stefani_et_al, Int. J.Heat Transfer
49, 3167-3179 (2006)[6]

D.Schwabe_et_al,

Phys.

Fluids,19,072102 (2007)

集合現象の数値解析を行い，実験結果と近い値を最後に，本研究制度（若手奨励枠）を活用得ることができた．しかし，周囲気体を考慮したさせて頂きましたことに大変感謝致します．

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

プログラム高度化支援事業研究報告

ディレクティブベースプログラミングによる

FIT の高速化と性能検証
中畑 和之 ∗∗

愛媛大学大学院理工学研究科

1 はじめにバイス間の通信に MPI を用いることが多い．

NVIDIA 社による統合開発環境 CUDA によって，GPU を用いた汎目的計算が一般技術者にも身近に択して頂き，有限積分法 [4](Finite Integration Tech-

なって久しい．しかし，ソースコードを大幅に変更しなければ性能が発揮できないため，プログラミングの初学者にとっては敷居が高い．これに対して，ディレクティブベースにコードを挿入する OpenACC[1]の規格が発表されている．これは，ソースコードに並列化を指示する行を挿入するだけで GPU 計算が可能となり，手軽で可搬性の高いソフトウェアの開発環境として注目されている．GPU 計算の台頭を受けて，Intel 社から Xeon Phi[2] が発売され，このコプロセッサ上で並列計算を行う環境 (Many Integrated

Core：MIC) も整備されつつある．Xeon Phi は，HPC
計算機の標準プロセッサとしても普及している Intel
Xeon プロセッサの浮動小数点演算コアを利用しており，これをプロセッサ・チップ上に多数搭載した

並列処理用の専用プロセッサである．また，GPU な

どのアクセラレータとは異なり，OS(Linux) を搭載

し，IP アドレスをもつ独立した分散ノードとして動作することが特徴の 1 つである．MIC Ofﬂoad[3] はOpenACC と同様にディレクティブベースにコードを

平成 24，25 年度にプログラム高度化支援事業に採nique：FIT) にイメージベース処理を組み合わせた超音波伝搬シミュレーション [5] の高速化についてア

ドバイスを頂いた．FIT は，支配方程式 (偏微分方程式) を微小四角形領域 (セル) で積分した後に離散化する手法であり，粒子速度と応力を時間ステップをずらして交互に計算して時間更新を行う．陽的に解の更新を行うため，マルチ CPU 使用時の並列計算効率がよいことが示された．平成 26 年度の高度化支援事業では，MIC Ofﬂoad によるディレクティブベース

プログラミングでイメージベース FIT の高速化を検討して頂いたので，それを報告する．また，比較のため OpenACC を用いて FIT を高速化したものを別途作成し，MIC Ofﬂoad のコーディングと比較を行った．本支援事業では 3 次元 FIT の MIC 版をチューニングして頂いたが，申請者自身のコーディングによる OpenACC 版と比較するため，2 次元 FIT に書き直したもので両者を比べることとする．FIT のベク

トル化する箇所や，ホスト-デバイス間のメモリ転送量・タイミングは MIC Ofﬂoad と OpenACC で同じとし，両者の速度の違いを調べた．

追加するだけで並列演算を行うことができる．さらに，Xeon プロセッサと開発環境や数値演算ライブラリなどが共用できることから，既存の HPC 環境との親和性が高く，パフォーマンスの向上が見込める．

OpenACC や MIC Ofﬂoad は大量のスレッドを利用し

2 アクセラレータの仕様とディレクティ
ブプログラミング
2.1

アクセラレータて並列処理することができる一方で，ホストとデバ

GPU や Xeon Phi デバイスは，ホストと独立した

イスのデータ転送がボトルネックになることがある

メモリ空間を持っており，基本的にはディレクティブ

ため，コーディングする際に配慮する必要がある．まで指示した部分をデバイスで高速に実行するというた，デバイス側のメモリ容量に制限があることから，意味でアクセラレータと呼ばれることが多い．MIC

必要に応じて複数のデバイスを用いる．このときデOfﬂoad 計算を行うために，京都大学スーパンコンピ

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 25ュータシステム E に搭載されている Xeon Phi 5210D

コプロセッサを用いた．Xeon Phi 5120D は，コア，キャッシュ，メモリ・コントローラ，PCIe クライア

ント・ロジック，双方向リング・インターコネクトから構成される．Intel 社の x86 アーキテクチャがベースのスカラユニットと 512 ビット幅の積和算ベクト
ルユニットを備えている．Xeon Phi 5120D は 60 個のプロセッサコアで，1 コアあたり 4 ハードウェア・
スレッドをサポートしているため，全体では 240 ス

レッドが利用可能である．デバイスは 8GB のメモリを持ち，各デバイス間はリング・ネットワークで接続され，これらは PCIe 経由で外部からコントロール
される．
一方，OpenACC による計算を行うために，GPUボードとして NVIDIA 社の Tesla C2075 を用いた．
自身の研究室で所有するワークステーションに Tesla

C2075 が実装されており，デバイスとホストのイン

ターフェースは PCI Express 2.0 である．Tesla C2075には 448 個のストリーミング・プロセッサ (SP) と呼ばれる最小単位の演算処理ユニットが存在する．この SP が 32 個集まって 1 つのストリーミング・マル

program testMIC
implicit none
real, dimension(:) allocatable :: x, y
integer :: n, i
!dir$ attributes offload : mic :: x,yn = 10allocate( x(n), y(n) )x(:)=1; y(:)=2!dir$ offload_transfer target(mic) &in (y,x,n) : alloc_if(.true.) free_if(.false.)!dir$ offload begin target(mic) &nocopy(y,x,n)!$ omp parallel
do i = 1, ny(i) = y(i) + x(i)end do!$ omp end parallel

チプロセッサ (SM) を形成する．SM にはワープスケ!dir$ end offload

が 2 個あり，2 つのワープを並列に実行することがで

deallocate(x, y )ジューラが 2 個と命令を行うディスパッチユニット

きる．GPU1 基あたりの DRAM は 6GB であり，大規模な計算を行う場合には，複数の GPU を使用することになる．特殊関数ユニット (SFU) では，正弦関数，余弦関数，相反関数，平方根などの命令を実行でき，各 SFU が 1 クロックで 1 スレッドの命令 1個を実行するため，8 クロックで 1 ワープが実行できる．

2.2 ディレクティブプログラミング
ディレクティブプログラミングの利点は，これら

MIC ǁȇȸǿ᠃ᡛǪȕȭȸȉር‫׊‬Ʒਦ‫ܭ‬

end program
図 1: MIC Ofﬂoad を用いたプログラムの例• ofﬂoad ディレクティブ
まず，attributes ofﬂoad ディレクティブは，モジュール
変数をホストとデバイスの両方のメモリに割り当てることを示すために用いる．次に，ofﬂoad transfer ディ
レクティブは，データをホストとデバイス間でやりとが使用できない環境では，その指示行は無視されるりするために用いる．また，デバイス上でのメモリ確点であり，並列環境と非並列環境でほぼ同一のソー保も兼ねている．図 1 のプログラムの ofﬂoad transfer

スコードを使用できることである．MIC Ofﬂoad はOpenMP などで並列化されたコードにディレクティ
ブを挿入し，その並列領域をアクセラレータ上で実行するモデルである．ここでは，MIC Ofﬂoad における Ofﬂoad ディレクティブについて説明する．MIC

Ofﬂoad によるプログラミング例を図 1 に示す．主に用いるディレクティブは，以下の 3 つである．
• attributes ofﬂoad ディレクティブ• ofﬂoad transfer ディレクティブ

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

ディレクティブは，配列 x，y と変数 n をデバイス側にコピーし，メモリの割り当てを行い，メモリ解放
はせずに常駐させるという指示である．ofﬂoad ディ
レクティブは，処理をデバイス側で実行させるために用いる．通常は，OpenMP 指示行によってループが並列化されている箇所を指定する．また，signal 節と wait 節と組み合わせることで非同期データ転送を行うこともできる．
次に，OpenACC について述べる．これは，アクセラレータ上で並列プログラミングを行うための標
準規格であり，CUDA などの独自言語を用いることなく，手軽に並列プログラミングを行うことができ
る．既存のソースコードに OpenACC ディレクティ

ブを追加すると，OpenACC 対応コンパイラは，ホス

示的に指示することができる．図 2 にある!$acc loopgang(32) vector(64) は，gang 分割 32 でベクトル長を 64 で並列分割を行うという意味である．worker は

1 つの実行ユニット内で細粒度の並列性の実行を行う単位であり，各 gang は 1 つ以上の worker を有す

ト (CPU 側) コードとデバイス (GPU 側) コードを自る．vector は，worker 内部で実行される SIMD(Single

文，kernels 構文，Loop 構文の 3 つからなる．図 2 に位である．gang，worker，vector をそれぞれ CUDA

動生成する．OpenACC のディレクティブは Data 構

OpenACC を用いた簡単なプログラミング例を示す．

Instruction Multiple Data) またはベクトル型の並列単で用いる計算単位で表すと，スレッドブロック，ワー

Data 構文はホストとデバイス間のデータ転送を明示

プ，スレッドとなる．

では，この処理の導入が最も重要となる．また，Data

3 ディレクティブを導入した FIT

構文によって既にデバイス上で存在するデータに対3.1

的に指示するためのディレクティブである．OpenACC

して，ホスト上のデータと明示的に同期をとりたい
場合は update 指示行を使用する．kernels 構文では，
デバイス上に並列計算するループ対象部分を指定す
る．これは，ベクトル化・並列化可能なループを対象としており，ループの直前にディレクティブを挿入することにより，コンパイラが並列化コードを生成する．指定した領域は GPU で実行されるカーネ
ルとして生成される．この構文では，ループのベク
トル長や並列分割の方法を節を使用してユーザが明program testACC(x3 ) の粒子速度を v3 (x, t)，せん断応力を τ31 (x, t)，
τ32 (x, t) とする．FIT は，運動方程式と構成則を領域
V で積分することによって，離散化式を導く．V の大きさを画像データのボクセル（ピクセル）サイズと等しくすれば，例えば CT 画像などを簡単に数値シミュレーションに取り込めるのが特徴である．時
間方向については，以下のような中心差分法で近似する．

implicit none
real, dimension(:) allocatable :: x, y
integer :: n, i
n = 10allocate( x(n), y(n) )x(:)=1; y(:)=2有限積分法 (FIT)FIT の詳細な定式化については，過去の論文 [5] で記載されているので，ここでは簡単に述べる．直交座標系の位置ベクトルを x=(x1 , x2 ) とし，面外方向

1

1

z

z−1

{τ }z+ 2 = {τ }z− 2 + ∆t{τ̇ }z ,{v} = {v}ȇȐǤǹǁȇȸǿ᠃ᡛ

(1)
(2)

ここで，z は整数，∆t は時間ステップ幅である．上式
は，速度 v3 と応力 τ31 ，τ32 を交互に時間更新する陽解法であることを示す．これらは，実際には τ31 (i, j)!$ acc data copyin(x(1:n)) copy(y(1:n))!$ acc kernels!$ acc loop gang(32) vector(64)do i = 1, ny(i) = y(i) + x(i)end do!$ end kernels+ ∆t{v̇}z− 12ȇȐǤǹܱᘍƷር‫׊‬や v3 (i, j) といった 2 次元配列を確保して，時間更新
する毎にメモリを上書きしながら利用することにな
る．ここで，i は x1 方向，j は x2 方向のインデック
スを表す．図 3 に FIT のフローを示す．更新の過程

で，速度 v3 は数値積分によって変位 u3 に変換してファイル出力している．また，MPI 通信では，v3 の1 レイヤー分の袖部を相互に転送する．

!$ acc end data

3.2

deallocate(x, y )3.2.1

高速化・効率化のための改良プログラムの基本的項目の修正

MPI と OpenMP を併用するため，デバイス－ホス

end program
図 2: OpenACC を用いたプログラムの例

ト間のメモリコピーのタイミングとスレッドの同期に注意する必要がある．高度化の中で分かったことは，4 バイト程度の小さなスカラー変数のコピーで京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 27

Startࣇ࢓࢖ࣝධຊ

v3 ࡢィ⟬
v3 ࡢ MPI ㏻ಙฎ⌮t3j ࡢィ⟬ (j=1,2)u3 ࡢィ⟬
Update

ることになる．

3.2.3

FIT は陽解法であるので，定期的に計算結果を

write する必要がある．非破壊検査への応用を考えて，粒子速度 (v3 ) を数値積分して変位 (u3 ) に換算し，それを出力する．全セル分の変位を出力すると膨大なファイル量になるため，実際には間引いて出力を行っている．高度化前のコードでは，この出力処理がネックとなっていた．出力する総数を L とす

ると，次のように出力バッファOutput にコピーしていた．

L=0
do j = 1, M, Kࣇ࢓࢖ࣝฟຊ

do i = 1, N, K
L=L+1

End

図 3: FIT のフローチャート

出力処理の最適化！x2 方向ループ (間引き：K)！x1 方向ループ (間引き：K)

Output(L) = u3(i, j)

である．また，計算のセル領域で超音波の波源を探

end doend do
ここでは，出力バッファの番地 L について依存性が索するためのフラグ配列があったが，領域数に比べあるため，ベクトル化の対象外となっていた．そこあっても相応のオーバーヘッドの原因となったことてほとんどヒットしないフラグについては，検索方法を変えるなどの必要があることがわかった．また，Xeon Phi ではベクトル化が重要であるが，間接参照
等を含む場合にはベクトル化できていても，あまりで，以下のように L でループを駆動し，L から出力する位置の u3 をコピーすることとした．i = 1 ∼ Nのうち Ni 個を取り出し，j = 1 ∼ M のうち Nj 個を取り出して間引くとする．このとき L = 1 ∼ Ni × Nj

良い性能が得られないこともわかった．

とし，依存性を解消するには，以下のように計算す3.2.2

ることにした．

分岐処理の見直し

do L = 1, Ni * Nj
i = mod(L − 1, Ni ) * K + 1
j = mod((L − 1) / Ni , Nj ) * K + 1

FIT では，自由反射境界 (空洞) や PML などの境界条件を扱うことができる．これは，v3 (i, j)，τ31 (i, j)，τ32 (i, j) のすべてのカーネル計算に組み込まれている．これらの更新計算は i と j の 2 重ループとなっ

Output(L) = u3(i, j)

理，PML の 3 つを場合分けしており，それぞれフend do
この方法であれば，コンパイラは依存性がないことをラグをチェックすることで計算するかどうかを判定判断できるため，より効率的なコードが生成できる．

している．空洞については，計算領域の任意の場所

計算順序およびデータ転送のフローており，これらのループ中では，内点処理，空洞処にあるため，その数もケースバイケースである．しかし，PML については常に使用する訳では無く，しかも PML は計算領域の四隅のレイヤー部であるため，この処理を主要ループから分離する手法を採用した．FIT は陽解法であるので，メモリを節約するため

v3 (i, j)，τ31 (i, j)，τ32 (i, j) の更新時は，バッファ配列などを設けずに自身の配列に上書きする．しかし，

PML 領域だけは別途配列を用意し，このバッファ配
列を v3 (i, j)，τ31 (i, j)，τ32 (i, j) にコピーする．従って，PML 領域は 1 ステップあたり 2 回計算が行われ

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

3.2.4

プログラム高度化後のフローチャートを図 4 に示

す．高度化前のプログラム (図 3) では，粒子速度と応力を計算した後に変位 u3 の計算を行っていた．高
度化後は，図 4 に示すように，変位 u3 をホストに転送しファイル出力する作業と，せん断応力のデバイ
ス側での計算をオーバーラップさせている．

Startࣇ࢓࢖ࣝධຊ
ධຊࢹ࣮ࢱ → ࢹࣂ࢖ࢫ㌿㏦

v3 ࡢィ⟬u3ࢹࣂ࢖ࢫฎ⌮1ࢹࣂ࢖ࢫ㸫࣍ࢫࢺ
࣓ࣔࣜ㌿㏦
࣍ࢫࢺฎ⌮

✵Ὕ
0.5

5.0 μsv3 ࡢ MPI ㏻ಙฎ⌮u3 ࡢィ⟬

10.0 μs
0.0

u3 → ࣍ࢫࢺ㌿㏦-0.5

t3j ࡢィ⟬ (j=1,2)ྠᮇࣇ࢓࢖ࣝฟຊ

-1.0

Update15.0 μsEnd

図 5: FIT による SH 波伝搬の計算例．

図 4: プログラム高度化後のフローチャート．

4 高速化性能の検証
計算例として，FIT による SH 波の伝搬の様子を

図 5 に示す．ここでは，図 5 のような正方形の領域を対象として，セル数が 10002 ，20002 ，40002 ，50002の場合における OpenACC と MIC Ofﬂoad の実行時間を比較する．なお，時間更新は各セル数において2000 ステップである．MIC Ofﬂoad 計算を行うために，京都大学スーパンコンピュータシステム E に搭載されている Xeon Phi 5210D コプロセッサを用い

20.0 μs

要している．セル数が増えると，OpenACC も MIC
Ofﬂoad も CPU 計算を上回るようになる．OpenACCの性能はセル数が 20002 以上では CPU 計算の 1.3 倍程度で頭止まりであるが，MIC Ofﬂoad はセル数の増加と共に性能が向上している．

5 まとめ
本事業では，デバイス－ホスト間のメモリコピーのタイミングを意識したプログラミングフローや，分た．また，OpenACC による計算を行うために，当方の研究室で有するワークステーション（NVIDIA 社

2の Tesla C2075 の GPU を搭載）を用いた．
ここで，性能(P) =CPU10 スレッドの計算時間
アクセラレータの計算時間2と定義したものを指標として用いる．京都大学のスパ
コンと，研究室のワークステーションではノード間のࢭࣝᩘ
2

通信環境が大きく異なるため，MPI 通信の比較は行っていない．CPU の 10 スレッド計算は，同じシステム

E のコア (Intel Xeon E5 系，クロック 2.5GHz) を用いた．図 6 に計算性能の比較を示す．OpenACC では，
vector(128)，gang(4) と指定し，MIC 実行時には，59
コア 236 スレッドを利用して計算を行った．セル数が1000 のとき，MIC Ofﬂoad も OpenACC も CPU 計
2

算よりも性能が小さいことがわかる．セル数が 20002のとき，OpenACC は CPU 計算よりも高速化しているが，MIC Ofﬂoad は CPU 計算の 2 倍程度の時間を0,&2IIORDG
2SHQ$&&

2SHQ032ᛶ⬟ 3
図 6: OpenMP(10 スレッド) の計算時間を基準とした
場合の MIC Ofﬂoad と OpenACC の計算時間の比較京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 29

岐処理の見直し，出力処理の最適化についてご指導を頂いた．これ以外にも，カーネルループのブロック
化，ロードバランスの改善など，様々な取り組みをして頂いた．3 次元版の FIT では，オリジナルコードを

CPU のみで実行した場合と比較して，MIC Ofﬂoad
では約 8.5 倍の高速化が実現できた．また，自身でOpenACC コードを作成し，MIC Ofﬂoad コードと比較を行った．当方の所属する研究室の GPU ボードのハードウエア上の制限もあって，2 次元版 FIT でこれらを比較した．セル数が増えると，OpenACC もMIC Ofﬂoad も高速化の効率が向上する．OpenACCの性能は，セル数が少ないときでも発現するが，MIC

Ofﬂoad はセル数の増加と共に性能が向上している．
今後の課題は，アクセラレータを複数使用した場合（CPU 間の通信が発生する場合）の高速化についても検討を行うことである．
謝辞： 有益なご助言を頂いた「スーパーコンピュー
タシステム共同研究企画委員会」の皆様，並びに，チューニングコードを提供して頂いたクレイ・ジャ
パン・インクの武田大輔様に感謝を申し上げます．

参考文献[1] 松岡聡，OpenACC を用いた大規模流体アプリケーショ
ンの高速化，学際大規模情報基盤共同利用・共同研究
拠点 平成 25 年度共同研究 中間報告書，2013．
[2] ジム・ジェファース，ジェームス・レインダース著/すがはらきよふみ訳，インテル Xeon Phi コプロセッサー
ハイパフォーマンス・プログラミング，カットシステ
ム，2014.
[3] Developer Zone，
https://software.intel.com/en-us/articles/theheterogeneous-programming-model/[4] P. Fellinger, R. Marklein, K.J. Langenberg and S. Klaholz, “Numerical modeling of elastic wave propagationand scattering with EFIT –elastodynamic ﬁnite integration technique”, Wave Motion, 21, pp.47-66 (1995).
[5] K. Nakahata，J. Chang，M. Takahashi，K. Ohira and Y.
Ogura，Finite integration technique for coupled acoustic
and elastic wave simulation and its application to noncontact ultrasonic testing, Acoustical Science and Technology, Vol.35, No.5, pp.260-268, 2014.

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

プログラム高度化支援事業研究報告

動的／静的水～土骨格連成有限変形解析コードの高度化～領域分割法の適用の試み～
野田利弘
名古屋大学 減災連携研究センター

1 平成 26 年度事業の概要

2 領域分割法の GEOASIA への適用

筆者らが開発した解析コード（GEOASIA® [1],図１は本解析コードに適用した際の処理手順フロ[2]）に対して、平成 24 年度（下半期）[3],・平成ーを示す。今回適用した領域分割法では、分割条件

25 年度[4],に続き、平成 26 年度も京都大学学術情報の指定に伴う、領域割当、付番変更、境界条件など

メディアセンターのプログラム高度化支援事業に採の付け替えを行う部分【前処理コード】、領域分割に択していただいた。この概要を報告する。

伴い実際に並列計算を行う部分【GEOASIA 計算】

本事業の目的は、大規模・高速化計算が実現できるよう当該解析コードの高度化・高性能化を図るこ

および出力のための集約を行う部分【前処理コード】の 3 つに分けて計算を実施した。

とであった。平成 24 年度においては、係数行列生
成部（KGRMAT）における OpenMP と MPI とのハイブリッド並列化を実現し、連立一次方程式の求解部分（ソルバ）に直接法として MUMPS を、反復法として PETSc などを組み込み、当該解析コードを拡張した。平成 25 年度においては、KGRMATにおける零対角要素の適切な処理のために PETSc
ソルバにおける前処理を施すとともに、性能向上に繋げるため、実行時スレッド数や Affinity 設定など
を確認・検討した。
平成 26 年度は、更なる並列化による高度化のた
め、KGRMAT およびソルバ関連サブルーチンに対する領域分割法に基づくプロセス並列化を適用した。
また、プロセス並列化のための領域分割方式、領域分割に対応した疎行列ソルバ、プロセス間通信方式

図１ 領域分割関連ファイルなど処理手順フローなどの検討を行うとともに、領域分割に基づくプロ
セス並列化の適用を行うに際しての技術調査を実施

領域分割法に関する動作検証を行うため、表１にした。ここでは、領域分割法の検証結果を中心に紹示す係数行列の特徴を有する 3 次元データを 4 種類

介する。なお、領域分割法による並列処理は、計算使用した。1) 3D（要素数：12,000，節点数：13,671）,

対象領域を複数の領域に分割し、各分割領域をプロ

2) 3DL（要素数：96,000，節点数：102,541）, 3) 煙セス並列化の手法により並列処理し、領域境界間で突 half（要素数：5,500，節点数：6,894）, 4) 煙突

必要な値の送受信をプロセス間通信により行う計算

full（要素数：11,000，節点数：13,239）である。

手法である。

図２はこれらの分割の一例を示す。2)は 1)の xyz 方京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 31

向をそれぞれ 2 倍の分割にしたもので、地表面中央に集中荷重が鉛直に与えられる。3)は 4)の半分の断面である。また、煙突部分は弾性体で地盤部分は弾塑性体で異なる要素になっていて、地震動が地盤底面で水平に与えられる。なお、一定時間間隔 DT で収束計算を行いながら逐次時間を進めているが、表１の「計算ステップ数」とはその時に回数である。
また、「制約条件」とは向かい合う側面で同一変位が生じるように土の変位場に課した条件を示す。
表１ テストデータ(係数行列：接線剛性マトリクス)
テストデ
ータ名

係数行列非
零成分数

解ベク
トル
成分数

計算ス
テップ
数

制約
条件

3DL

3,612,682

55,38010

有3D

28,756,582

413,16010有／無煙突 half

1,543,042

24,1821000

有

煙突 full

3,141,416

49,0841000

有3) 煙突 half（2 分割：x2×y1×z1 の例）

4) 煙突 full（4 分割：x4×y1×z1 の例）
図 2 動作検証に用いた 3 次元データ

3 実行性能計測結果
図３～図６には各データケースについて、全体とソルバ部の計測時間の一例を示す。図３と図４の 1）

1) 3D（要素数：8 分割：x2×y2×z2 の例）と 2)のデータについては PETSc でのソルバ部計測も抜粋も示す。縦軸はソルバライブラリ（棒グラフは PARDISO：PA、MUMPS ：MU、PETSc：PE
の違い）
と並列度
（PARDISO はスレッド並列のみ、
MUMPS、PETSc を用いた領域分割で z 方向以外の分割を含む場合は括弧内に示す）、
上横軸は計測時
間（秒）を表す。PETSc ソルバ部では下横軸にプロファイル計測を行った際のソルバ部計測箇所
SOLVE(2)（PETSc 演算部）の（PETSc の p1t1 に対する）速度向上比を示している（折線）。なお、本
計測で使用可能な計算リソースが８ノードまでであったため、p8 以下は１ノード１プロセス、p16 以上

2) 3DL（32 分割：x2×y2×z16 の例）

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

はノード内複数プロセスとしている。また pmtn は
m ヶのプロセス並列と n ヶのスレッド並列を示す。

1)と 2)のデータについては、特に PETSc による並列化効果が見られている。各サブルーチンでの速度向上比の詳細は省略するが、2)では p64t2(x と y
方向 2 分割、z 方向 16 分割)時に、シリアル（p1t1）
比で、係数行列作成ルーチンKGRMAT で約120 倍、
主要ソルバで約 12 倍、全体で約 23 倍速度が向上した（なお、PARDISO の p1t1 比ではソルバ部で 1）が 222 倍、2)が 1092 倍の速度向上比である）。
3DL データで制約条件無のケースは、計測結果を省略するが、全般的に制約条件有のケースに比べ早い傾向があり、有のケースの場合に収束しなかった
ケース（1）,2)の x,y 方向 4,4 分割のケース）でも、収束が見られた。
次に、3)と 4)のデータについては、PETSc を使用
した場合、反復法収束計算における収束が得られなかったため、MUMPS を適用している。MUMPSは PETSc に比べ求解処理に時間がかかるものの、
4)のデータの場合、p4t16（水平方向 4 分割（16 ス
レッド並列）、鉛直方向分割なし）時に、p1t1 比で
KGRMAT が約 68 倍、主要ソルバが約 8 倍、全体で約 13 倍速度が向上した。

図４ 性能計測結果（3DL データ）

図５ 性能計測結果（煙突 half データ）

図３ 性能計測結果（3D データ）

図６ 性能計測結果（煙突 full データ）

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

4 おわりにえている。
謝辞 京都大学学術情報メディアセンターの中島

平成 26 年度では、当該解析コードの係数行列作

浩教授および牛島省教授ならびにクレイ・ジャパ

成やソルバ関連のサブルーチンに加えて、プロファン・インクの方々には、上記以外の点においてもごイル計測での上位のサブルーチンなどに対して領域尽力いただきました。ここに深く謝意を表します。

分割法を適用し、比較的シンプルな形状を有する解
析対象の入力データを用いて実行性能解析を行った。

参考文献

領域分割法を適用した場合とオリジナル解析コー[1] Asaoka, A. and Noda, T. (2007): All soils all states all

ドのPARDISOシリアルの場合の解の精度の比較に

round geo-analysis integration, Int. Workshop on

ついては省略するが、いずれの入力データ、並列化

Constitutive Modelling - Development, Implementation,

条件のもとで殆ど差がないことを確認した。

Evaluation, and Application, Hong Kong, China, 11-27.

一方で、入力データの内、煙突データについては[2] Noda, T., Asaoka, A. and Nakano, M. (2008):

PETSc での収束は得られなかったこと、中間処理コ

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),

平成 27 年度は、筆者のミスで本事業に申請でき

771-790.

なかったが、
上述のような課題が残っている。
また、[3] 京都大学学術情報メディアセンター全国共同研

平成 25 年度終了時点での解析コードによっても、

究利用版広報，Vol.12, No.2, 33-35.

平成 25 年度の研究成果[4]に加えて、図７と図８の[4] 京都大学学術情報メディアセンター全国共同研

ような解析結果[5]などが新たに得られている（計算究利用版広報，Vol.13, No.2, 44-47.

条件を含む詳細は文献[5]を参照）。これは、解析コ[5] 豊田智大・野田利弘・山田正太郎・浅岡顕(2015):

ードが高速化されていて、系統的な計算ができる点

横ずれ断層に伴う地表面のリーデルせん断帯の発生

が大きい。
今後、
更なる研究を進めていくためにも、に及ぼす材料的初期不整の影響，第 50 回地盤工学会

改めて平成 28 年度に向けて応募してゆけたらと考研究発表会，2013-2014.

(a) 弱点 1 箇所(b)

弱点 1 箇所(c)

弱点 1 箇所(d)

強点 2 箇所

図７ 初期不整に応じたリーデルせん断帯発生状況の違い(a) 弱点 1 箇所(b)

弱点 1 箇所(c)

弱点 1 箇所

図８ 初期不整に応じたフラワー構造発生状況の違い

34 − Academic Center for Computing and Media Studies, Kyoto University(d)

強点 2 箇所プログラム高度化支援事業研究報告

破壊力学に基づく損傷モデルを用いた鉄筋コンクリートの 3 次元破壊シミュレーション
車谷 麻緒 ∗∗

1

茨城大学 工学部 都市システム工学科と剛性低下により，近似的に 3 次元ひび割れ進展解

はじめに
土木構造物の設計や維持管理の高度化において，

鉄筋コンクリートの力学性状の評価は重要である．
鉄筋コンクリートには大小様々なひび割れが発生す
る．鉄筋コンクリートに生じるひび割れは，部材の力学性能に影響を与えるだけでなく，力学性能に影響しないような微細なひび割れであっても，水やイオンの移動経路となって，部材の耐久性の低下に起因する．
コンクリートに発生するひび割れをモデル化およびシミュレートするには，時々刻々と移動する不連
続な境界を追跡する必要があり，一般に困難な数値
シミュレーションとなる．コンクリートの場合は，多数のひび割れが複雑に発生・進展することになるので，その挙動を 3 次元で詳細に再現するのは容易ではない．実際に，鉄筋コンクリート内部におけるひび割れ進展挙動を 3 次元で精度よく詳細に再現した

研究例はほとんど見当たらない．ひび割れ進展を 3次元で詳細に追跡でき，かつ実験結果を定量的に再現できる数値解析手法の開発は，土木工学と計算力
学分野における重要な研究課題であり，鉄筋コンク

リートに生じる内部ひび割れを 3 次元で詳細にかつ定量的に再現できれば，土木構造物の維持管理を高度化することもできる．
そこで本研究では，準脆性材料の破壊力学に基づいて，材料の剛性を低下させる損傷モデルを定式化
し，これを有限要素解析に導入することにより，鉄筋コンクリートのひび割れ進展挙動を 3 次元で詳細に再現できる数値解析手法を構築する．さらに，解析プログラムの並列化・高度化を行い，大規模シミュ
レーションと組み合わせることにより，損傷モデル

析を行うことができる．通常の FEM がベースである
ので，FEM で培った技術はすべて利用することができ，信頼性のある数値解析を行うことができる．これにより，通常は目視できない鉄筋の周囲に発生するコンクリート内部のひび割れを 3 次元でシミュレー
トすることができ，土木構造物の維持管理に対する新たな知見を得ることができる．開発した方法を鉄筋コンクリートの 3 次元ひび割れ進展解析に適用し，大規模数値解析への適用性を検討するとともに，有限要素解析におけるメッシュサイズの依存性についても検討する．

2

数値解析手法
コンクリートのひび割れ進展挙動のモデル化には，準脆性材料の破壊力学を考慮した修正 von-Mises モ

デルに基づく等方性の損傷モデル [1] を適用する．このモデルは，コンクリートの破壊進行過程における結合力と開口変位の関係を連続体損傷モデルに応用したもので，コンクリートの破壊力学の特徴である
破壊エネルギーの定義も満足するモデルである．これを通常の有限要素法に導入することにより，コンクリートのひび割れ進展挙動を再現する．
多次元での損傷を評価するに際して，ひずみテン
ソルの不変量を用いた等価ひずみを用いる．具体的には，次式で与えられる修正 von-Mises 基準に基づ
く等価ひずみを採用する．
√
(
)2
k−1 ′
12k
k−1
1
′
I1 +
I1 +
J′
εe =
2k(1 − 2ν)2k
1 − 2ν
(1 + ν)2 2(1)

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

450
400

3D view
450(Plane strain condition)

150

Principal strain ε 2

1200
200 160

y

0

x
z

D16

k=1

110

Material parameter of concrete
E

ν

κ0

Gf

k

30 GPa0.2

0.0001

0.1 N/mm

20

Material parameter of steel

k=2
E

ν

200 GPa

0.2

0Principal strain ε 1

図 2: せん断補強筋のない RC はりの解析条件と材料

図 1: 2 次元および 3 次元での主ひずみ空間におけるパラメータ

等価ひずみの等値面上式において，ν はポアソン比，k は圧縮に強く引張に弱い特性を付与できる圧縮引張強度比である．I1′

3

数値解析例
本研究で構築した破壊シミュレーション手法では，はひずみテンソルの第 1 不変量，J2′ は偏差ひずみテ

等方性の損傷モデルを用いて，コンクリートのひび

3 次元および 2 次元の主ひずみ空間に εe の等値面をプロットした例を図 1 に示す．圧縮引張強度比 kルを用いた有限要素解析では，解析結果がメッシュ

ンソルの第 2 不変量である．

を用いて引張破壊と圧縮破壊を考慮したことにより，圧縮に強く引張に弱いというコンクリートの力学特性を表現できることが分かる．
損傷は回復しないこととし，変形履歴における等価ひずみの最大値を κ ≥ 0 で表すことにより，等方性の損傷モデルは次式で与えられる．

σ = (1 − D(κ)) c : ε(2)D(κ) は損傷の度合いを 0 ≤ D ≤ 1 で表す損傷変数であり，次式で表される．

D(κ) = 1 −()

Eκ0 heκ0(κ − κ0 )
exp −
κ
Gf

(3)

ここで，κ0 は損傷開始時の等価ひずみである．損
傷モデルによる非線形有限要素解析では，割線法と割れ進展挙動をモデル化している．既往の損傷モデの寸法や分割に依存することが知られており，これを回避するために非局所モデルを用いることが多い．
以下では，RC はりの破壊シミュレーションを対象に，メッシュサイズを小さくしていき，コンクリー
トのひび割れ進展挙動に対するメッシュサイズの影
響について検討する．

3.1

解析対象と条件

解析対象は，図 2 に示すような，せん断補強筋（スターラップ）のない RC はりの 4 点曲げ試験である．

主鉄筋は D16 の異形鉄筋とし，数値解析では異形鉄筋の幾何形状を忠実に有限要素メッシュに反映させることとする．メッシュ分割に対するシミュレーショ

ン結果の客観性を検証するため，全体を約 200 万，

400 万，800 万，1600 万要素に分割した 3 次元有限

Newton-Raphson 法を組合わせた陰的な解法により，要素モデルを数値シミュレーションの対象とする．有コンクリートのひび割れ進展挙動を解析する．

は，強制変位 2.5 mm を 200 ステップで載荷する変36 − Academic Center for Computing and Media Studies, Kyoto University

限要素には，1 次の四面体要素を用いる．解析条件

Load (kN)

100

2,000,000 elements
500: 2,000,000 elements: 4,000,000 elements: 8,000,000 elements: 16,000,000 elements
02

4,000,000 elements

4

Displacement (mm)

図 3: メッシュサイズの異なる有限要素モデルに対す

8,000,000 elements

る荷重－変位関係の比較

位制御試験とする．
鉄筋とコンクリートの境界面のモデル化について説明する．本研究では異形鉄筋の幾何形状を詳細に

16,000,000 elements(a) Strain distribution at loading step: 100/200

モデル化していることに加えて，節が表現できる程度の小さなメッシュサイズを与えている．そのため，特殊な界面要素や付着モデルを用いずに，境界面に位置する有限要素の応力集中によって，境界面から

2,000,000 elements

損傷が発生し，付着力が徐々に消失する過程を再現する．

一般に，せん断補強筋のない RC はりでは，鉄筋はほとんど塑性降伏せずに，コンクリートの破壊の4,000,000 elements

みが進行し，せん断破壊に至ることが知られている．
加えて，通常の有限要素解析においては，鉄筋の塑性変形はメッシュの寸法や分割に依存しないことから，ここでは鉄筋を線形弾性体として，鉄筋の降伏

8,000,000 elements

直前までを数値解析の対象とし，コンクリートのひび割れ進展挙動に対するメッシュ分割の依存性について比較・検討する．
16,000,000 elements

3.2(b) Strain distribution at loading step: 200/200

解析結果と考察

0.0001

解析結果として，荷重－変位関係を比較したもの

を図 3 に示す．メッシュサイズによって解析結果に若干の差異はあるものの，全体としての傾向はほぼ

0.1

図 4: メッシュサイズの異なる有限要素モデルに対す
るひび割れ分布の比較

一致しており，メッシュ分割の相違が荷重－変位関係に与える影響は非常に小さいといえる．
次に，ひび割れの発生と進展の相違について検討

現する方法であるので，ひずみの分布がひび割れのする．モデル表面での等価ひずみの分布を比較した分布に対応する．どのようなひずみを可視化してもものを図 4 に示す．提案手法による 3 次元破壊シミュ
レーションは，ひずみの局所化によりひび割れを再

よいが，ここでは数値解析で利用している等価ひず
みの分布を示している．メッシュサイズを小さくし京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 37

ていくと，ひび割れの本数や経路が変化しそうであるが，結果を見て分かるように，メッシュサイズを小さくしても，ひび割れの本数や経路はほとんど変化していない．その要因として，破壊エネルギーを考慮した損傷モデルを用いていることと，陰解法に基づく非線形有限要素解析を行っていることが挙げられる．

Damage distribution at load step: 100/200

生じたひび割れの分布に着目すると，どのメッシュパターンにおいても，最初に曲げひび割れが発生した後，せん断破壊を示す典型的な斜めひび割れに発展する様子がシミュレートされている．以上の数値解
析結果より，本研究でコンクリートの材料モデルに適用している等方性の損傷モデルと，これを有限要素解析に導入した破壊シミュレーション手法は，メッシュサイズに対する依存性が小さく，提案手法によDamage distribution at load step: 200/200

0.0

1.0

る破壊シミュレーションは客観性のある解析手法であることが分かる．

図 5 は，有限要素モデルのなかで最もメッシュサ
イズの小さい 1600 万要素モデルに対して，モデル内部の損傷分布を可視化した結果である．本研究では，異形鉄筋の幾何形状を忠実に考慮してメッシュ分割

図 5: メッシュサイズの小さい有限要素モデルの内部における 3 次元損傷分布

せた数値解析は，強い非線形性により不安定な数値を行っているため，提案手法による破壊シミュレー

解析になることが予想される．今後は，鉄筋の弾性れ）に加えて，RC はり内部の鉄筋の表面から発生・トの両方の非線形力学挙動を 3 次元で詳細に再現でションは，RC はりの表面から生じる損傷（ひび割

限界後の塑性変形をモデル化し，鉄筋とコンクリー

進展する損傷（ひび割れ）までを解析できる方法できる解析プログラムを開発し，数値シミュレーショ

あることが見てとれる．

4

おわりに
本研究では，主にコンクリートの 3 次元ひび割れ

進展解析プログラムの開発に重点をおき，鉄筋コン
クリートにおける鉄筋には弾性モデルを適用した．
鉄筋が塑性変形に達しない問題を対象に，本手法を適用した結果，実際の実験結果に近いひび割れ進展挙動を 3 次元で詳細に，かつメッシュサイズに影響されることなく再現することができた．本研究の成果 [2] は，学術雑誌に投稿中であり [3]，図 5 の数値
解析結果は，可視化ソフトウェアの可視化事例として掲載されている [4]．

実際の RC 部材の破壊挙動は，コンクリートのひび割れ進展挙動と鉄筋の塑性変形が組み合わさったものである．鉄筋とコンクリートはまったく異なる材料であり，ひび割れ進展挙動と塑性変形もまった
く異なる力学挙動であることから，両者を組み合わ

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

ンにより実験結果を定量的に再現できるかどうかを検討する予定である．

参考文献[1] 車谷麻緒，寺田 賢二郎，加藤準治，京谷孝史，樫山
和男，「コンクリートの破壊力学に基づく等方性損傷モデルの定式化とその性能評価」，日本計算工学会
論文集，pp.20130015，(2013).
[2] M. Kurumatani, S. Abe, Y. Nemoto, Y. Soma, N.
Henmi, A method for fracture simulation of reinforced
concrete using a damage model based on fracture mechanics for concrete, The 3rd International Workshops
on Advances in Computational Mechanics, (2015).
[3] 車谷麻緒，根本優輝，相馬悠人，寺田 賢二郎，
「コン
クリートの破壊力学を考慮した鉄筋コンクリートの3 次元破壊シミュレーションとその性能評価」，日本計算工学会論文集，投稿中.
[4] CYBERNET AVS 可視化事例（材料），
http://www.cybernet.co.jp/avs/example/category/

プログラム高度化支援事業研究報告

都市の人口規模・空間分布における秩序形成森知也
京都⼤学経済研究所
18

1 はじめに

17
16
15

先進諸国の多くでは 8 割以上の経済活動が都市

New York
Los Angeles
ChicagoDallasPhiladelphia

BostonDetroit
San Francisco
Pheonix
Seattle

3

5

14
13

に集積しており、これらの国々における地域・国際

12貿易を始めとする経済地理は、都市群の経済的構造

10

11

9

と空間パターンで概ね特徴付けられることを意味す0る。興味深いことに、都市の⼈⼝規模、数、位置、産

1

2

4

6

図2 都市規模分布のフラクタル構造

業構造の関係は顕著な秩序によって関連付けられることが知られている (Mori et al. [2])。図 1は、⽇
独⽶の都市規模分布、および、アメリカついては、4⼤都市 (ニューヨーク, ロサンゼルス, シカゴ, ダラ

2 モデルと数値解析における課題ス) 経済圏の都市規模分布を⽰している。地域経済
圏は物流センサスを⽤いて同定されたもので、アメ分析には、多産業・多地域に拡張した Pflüger [3]

リカのような⼤国では複数の同規模の地域経済圏、モデルを⽤いた。異なる産業は異なる種類の消費財

⽇本やドイツ程度の規模なら国全体が概ね⼀つの地を⽣産し、各種消費財はさらに多数の差別化された域経済圏で構成されている。図が⽰すように、地域

財バラエティで構成される。各産業の集積は、財の経済圏単位では、都市規模分布の上部は概ね共通の多様性を嗜好する消費者、⽣産における規模の経済、べき乗則に従っている。さらに、このように同定さおよび、輸送費⽤の存在から⽣ずるフィードバック

れた地域経済圏は階層構造を持ち、図 2に⽰される効果によって形成され、集積の規模と空間的周期は、ように、都市規模分布は地域経済圏間でフラクタル

フィードバック効果の⼤きさにより異なる。産業の構造を持つ (Akamatsu et al. [1] )。本研究では、こ集積はその⽣産財の消費者の集積を伴って都市を形のような都市規模分布の空間的フラクタル構造を再成する。⼤都市は、共通の消費者を介して、複数の現できるミクロ経済学モデルの構築に取り組んだ。

産業が集積することにより形成される。
仮想経済では、円周上に⼀定間隔で配置した地
点群からなる⽴地空間を仮定し、所与⼈⼝の労働者 (= 消費者) を各産業・地点にランダムに配置した初期状態から、労働者の産業間・地点間分布を所与とした “短期均衡” の達成と、労働者の効⽤⽔準改善を⽬的とした移住を繰り返し、最終的に各労働17Japan (2000೥)
Germany (2008೥)US (2007೥)New York (2007೥)Los Angeles(2007೥)Chicago (2007೥)

16
15
14

者が移住の動機を持たない “⻑期均衡” に到達する。
各短期均衡の導出は、産業数 × 地点数の未知数にਓ‫نޱ‬໛
ʢର਺ʣ13

対して、産業毎に (地点数 × 地点数) の⾏列演算を12繰り返す do-loop を複数含む。今回のプログラム⾼

11
10
9

1

2

3

4

5

6

7

度化は、これら do-loop の OpenMP と MPI を⽤

8ਓ‫نޱ‬໛ॱҐʢର਺ʣ

いたハイブリッド並列化・円周⽴地空間から⽣ずる

図1 ⽇⽶独の都市規模分布

対称性を活かした演算の効率化を中⼼に⾏われた。

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

160

3 結果

140
120
100

3.1 計算性能ਓ‫نޱ‬໛ 80

図 3は、プログラム⾼度化による計算効率の改善

6040を、図 4は、ノード/コア数と実⾏時間の関係を⽰し2000ている。いずれも、(地点数, 産業数)= (512, 256) の

100

200

300

400

500

Ґஔ

下で、5 万回の短期均衡計算の所要時間で評価して

図5 都市の位置と⼈⼝規模いる。筆者によるオリジナル･プログラムに対して、⾼度化プログラムでは、より広範囲での並列化・効図 6は、⻑期均衡における、代表的な 4 産業につ率化が⾏われ、4 〜 5 倍程度の速度改善があった。

いて各⽴地点での (潜在的な) 利潤⽔準を⽰してい並列化は産業⽅向にのみであるため、並列化効率はる。⻑期均衡で実現している各産業の利潤⽔準を 1

総コア数 < 産業数の場合に限られるが、1 ノード・に正規化しており、産業⽴地は利潤⽔準 = 1 の地

28 コアの場合に⽐べて最⼤で 3 倍程度⾼速化した。

点のみで起こり、それ以外には⽴地は無い。産業

1, 2, 3, 4 の順で、供給財の差別化程度は低く設定されている。図が⽰す通り、より差別化の程度が⾼い

600
ΦϦδφϧ
ߴ౓Խ‫ޙ‬

480

財を供給する産業の⽴地点には、必ず、より差別360化の程度の低い財を供給する産業も⽴地しており、240

産業間で集積の空間的同期が起こっている。特に、

120

産業 1 は図 5に⽰される 3 ⼤都市にのみ⽴地して࣮ߦ࣌ؒ ඵいる。

0ϊʔυ਺ ʷ ίΞ਺  ϊʔυ

1.00

図3 パフォーマンス

0.98
0.96

૬ରར५ 0.94
0.92
0.90

࢈‫ۀ‬1

0.88

࢈‫ۀ‬2

0

࢈‫ۀ‬3

100

࢈‫ۀ‬4
200

300

400

500

Ґஔ

図6 産業集積の空間的同期現象

図 7は、実際の地域経済圏と同じ⽅法で同定した

120

仮想経済における地域経済圏単位の都市規模分布を10080⽰している。図 2で⽰した現実の都市規模分布と同60

様に、空間的なフラクタル構造が発現する。(点線40

は経済全体の都市規模分布を⽰す。)

20
0

6
4

図4 並列化効率2ਓ‫نޱ‬໛ ର਺

0-2-4

3.2 シミュレーション結果

-6

0

1

2

3

4

5

6

7

ਓ‫نޱ‬໛ॱҐ ର਺

上記の計算効率⽐較に⽤いたシミュレーションの図7 都市⼈⼝規模分布

結果として、所与⼈⼝の労働者を産業 × ⽴地に対してランダムに配置した初期分布から導出した⻑
期均衡を⽰す。図 5は⼈⼝の空間分布を表しており

4 おわりに(但し、⽴地空間は円周であり、地点 0 と地点 512は同⼀である)、その空間パターンは近似的にフラ本共同研究では、現実の空間経済において観察さクタル構造を持っていることがわかる。

れる、都市規模分布の空間的フラクタル構造と共通

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

べき乗則の創発を再現する数値解析プログラムを開発した。特に、計算の並列化と効率化によって、⼤規模なモンテカルロシミュレーションが可能になった。今後は、⽴地空間を 2 次元に拡張し、より現実的な設定の下で、現実に観察される秩序形成の再現を試みる予定である。

参考文献[1] Takashi Akamatsu, Tomoya Mori, and Yuki
Takayama. Spatial coordination among industries and the common power law for city size
distributions. Technical report, Institute of
Economic Research, Kyoto University, 2015.
[2] Tomoya Mori, Koji Nishikimi, and Tony E.
Smith. The number-average size rule : A newempirical relationship between industrial location and city size. 48(1):165–211, 2008.
[3] Michael Pflüger.

A simple, analyticallysolvable, chamberlinian agglomeration model.
Regional Science and Urban Economics,
34(5):565–573, September 2004.

3
京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 41

スーパーコンピュータ共同研究制度（大規模計算支援枠）研究報告

高次精度差分法による高レイノルズ数チャンネル乱流場の大規模直接数値シ
ミュレーション

山本 義暢
山梨大学大学院総合研究部

1 はじめにラを仮定した温度場の輸送方程式である．
空間離散化手法は，主流及びスパン方向に 10

非定常・3 次元構造を有する乱流現象を解析す

次精度の中心差分法3)，
壁垂直方向に 2 次精度の中る手法としては，
その基礎方程式を数値的に解き，心差分を使用する．時間進行は，対流項及び粘性

数値解を求める直接数値シミュレーション

項に 2 次精度Adams－Bashforth法，圧力項にEuler(Direct Numerical Simulation, DNS)が有力であ陰解法を適用し，Fractional step法により解く．こる．このDNSは渦運動の最小スケールから最大スの際圧力ポアソン解法には，
高速フーリエ変換 ( x,

ケールを包括する格子が必要となるため，その計

z 方向 ) と 3 重対角行列解法 ( y方向 )による直算負荷は莫大となる．しかし計算機性能の発展に接解法を使用する．境界条件は，主流及びスパン伴いその適用範囲は拡大している．特に近年高レ方向に周期境界条件，壁面でno-slip条件を課した．

イノルズ数境界層流れにおける大規模構造1),2)が温度場に対しては壁温一定 (上壁温度:θtop > 下壁

注目を集めている．この大規模構造は境界層厚で

温度:θbed) とし，主流及びスパン方向には周期境スケーリングされる構造を有するとともに乱流統界条件を課した．

計量への影響2)も指摘されている．室内実験による研究ではこの大規模構造は，壁面摩擦速度(uτ)と境界層厚(h)に基づくレイノルズ数（摩擦レイノ

た．
表 1 に本研究でのレイノルズ数 ( Reτ )，主流方
向，壁垂直方向，スパン方向 ( x，y，z )の計算領

ルズ数:Reτ = uτh/ν, ν: 動粘性係数）において 7300

域 ( Lx，Ly，Lz )，格子点数 ( Nx，Ny，Nz )，及程度に達するとその影響が顕著になることが報告

び格子解像度 ( ∆x+，∆y+，∆z+ ) を示す．なお，されている2)．
本研究ではこの大規模構造解析を目的として，

上付き添え字 + は，摩擦速度 ( uτ ) と動粘性係
数 ( ν )で規格化した値であることを示す．

京都大学学術情報メディアセンターCray・XE6・
256 ノードで実行可能な最大レイノルズ数： Reτ=4000 のDNSを実行した．

2 直接数値計算の概要

本研究における直接数値計算の解析対象は，図1 に示す，一定の圧力勾配で駆動される 2 次元チ
ャンネル乱流場である．また基礎方程式は，非圧縮性流体の運動方程式，連続式及びパッシブスカ

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

図 1 計算体系と座標系

元計算の最外ループ ( y 軸 ) で領域分割を行っ

表 1 計算条件
Reτ
present
.4)
5)

4000
4200
4079

Domain(Lx/h, Ly/h,Lz/h)16, 2, 6,42π, 2, π
6π, 2, 2π

Grid number( Nx,Ny,Nz )Resolution( ∆x+, ∆y+, ∆z+ )ている関係上，最外ループをさらに OpenMP で分

5760, 2048, 3072
2048, 1081,2048
8192, 1024, 4096

11.1, 0.5-8.0, 8.3
12.8, 0.1-8.0, 6.4
9.4, 0.1-54 , 6.2

限されるため，OpenMP における collapse 機能を割すると総並列数が結局 1 次元方向の格子数に制使用して，これを回避した．また 2D-FFT 部分においては，FFTW3.3 の OpenMP 版を適用した．

また表 1 には参考のため，同規模レイノルズ数

・メモリ量削減
10 次精度差分法では，1 点の微係数を計算するの直接数値計算条件を併記している．

際に， 10 点（拡散項では，20 点）のステンシルを要する．これに対処するために境界部分にステ

3 並列化

ンシルを拡張すると大幅なメモリ増加をもたらす
ため，境界部分についてはリスト配列を使用し，

3.1 領域分割

メモリ量は高次精度差分法においても，2 次精度

本研究では，壁垂直方向 (y 方向 ) への MPI

差分法の場合と同等となるようにした．

( Message Passing Interface ) による領域分割による並列化を適用する．この際に否通用となる通信
は，
壁垂直方向の差分演算時に必要となる袖通信，
圧力ポアソン方程式の 3 重対角行列解法時に必要となる all-to-all 通信となる．袖通信においては，
隣接間通信となるが，通信データを連続データとするために，プログラミング上，y 軸方向は，3次元配列の最外ループとした ( fortran の場合 )．
さらにスーパーコンピュータ「京」のような 3 次
元トーラスでるが，隣接プロセスが必ずしも隣接ノードとはならない場合があるため，ノード配置
マップを手動で作成し，隣接プロセスが隣接ノードとなるように工夫した（ただし Cray XE6 では
ユーザがノード配置を指定できないため、本設定は使用していない）．
またall-to-all 通信は，3 重対角行列解法時に，分割方向 ( y 方向 ) に再起参照関係となるために，これを回避するために分割軸の入れ替えを行うが，この際に必要となる．従来の並列化におい

図 2 1 次元領域分割ては分割軸の入れ替えをz軸方向に行っていたが，（FFT 及び TDMA の場合）

チャンネル乱流場ではy 軸とz軸方向の格子数が異なり，各プロセス間での通信量及び演算量のア

4 並列性能ンバランスが生じる可能性があった．
そこで今回，
分割軸の入れ替え時に転送配列を(x, z, y) の 3 次

表 2 にCray・XE6/256 ノード及び東京大学

元配列から，(xz,y ) の 2 次元配列とし，NxNz の格

FX10，理化学研究所「京」での実行演算速度(Rmax)

子数をMPIプロセスにより領域分割した．

を示す．理化学研究所「京」では 21 Tflopsの実行

・OpenMP

演算速度が得られており，これは理論性能(Rpeak)

また通信量の削減と並列ノード数を稼ぐことをの約 8.2%に相当する．Cray・XE6 の場合は，実目的として，ノード内は，OpenMP による並列化

効演算速度が約 5Tflopsであり，実行効率は 6.1%を計算コードのほぼ全領域に適用した．また 3 次とやや低めの値となった．この原因はノード間通京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 43

信（特にall-to-all通信）の性能がやや落ちること

法に基づくBernardini et al. (4) の結果は良好に一致

及び演算部分の実行性能もやや低めとなっており，したものの空間離散化精度が高い，Duran &全体として 6%台の実行効率に留まってしまった

Jimenez (3) の計算とは，主流方向乱流強度においものと考えられる．

て壁遠方 (y+ > 800)にて違いが見られ，彼らの結果は本計算結果より大きな値となっている．

表 2 実行演算速度の比較

5 計算結果と考察

図 5 散逸スペクトルの比較（最小スケール渦の再現性）

図 3 平均速度分布の比較

図 6 バッファー層における主流方向乱流強度の比較
そこでy+ = 800 の高さにおける散逸スペクトル
分布を図 5 に示す．本計算結果においては，解像度不足の影響は見られないため，上記乱流強度分
布の違いは，離散化手法あるいは解像度の影響によるものではないと考えられる．この原因は，
3)

Lozano-Durana & Jimenez のケースは若干レイノズ

図 4 乱流強度分布の比較数が高い(Reτ =4200)こと，また計算領域が小さい
ことによる大規模構造の過剰評価が影響している

次に計算結果の比較を図 3 及び 4 に示す．図 3のではと，推測される．また図 6 に示すバッファ

における平均速度分布については，ほとんど差異ー層付近の主流方向乱流強度分布においては，2は見受けられなかった．また図 4 に示す乱流強度

次精度中心差分法が適用されているBernardini et

分布においては，本計算結果と 2 次精度中心差分

al. 4) の結果では，100 > y+> 30 の領域で過剰評価

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

が見られた．これは低次精度の差分法を適用した

Tflops (実行効率約 6%)が得られるとともに，

場合に見られる特徴である．一方，本研究（10 次Reτ=4000 の数値解を得ることに成功した．同レイ

精度差分法）の場合，上記の影響はないものの、

ノルズ数の結果と比較し，基本特性はよく一致す

ピーク値が過小評価されていることがわかる．こることが確認された．しかし主流方向乱流強度分

れは主流方向の格子分解能がやや低いためと思わ

布の過小評価及びスカラー場の未発達が生じておれる．

り，解像度の上昇と時間積分長を増やす必要性が大規模構造の特長である，スペクトル分布の第2)あることが判明した．時間積分長の増加は計算機

二ピーク を確認する．図 7 にpre-multiplied スペ

占有時間の増加となり，本研究のような大規模ジ

クトル（波数をかけたスペクトル）を示す．室内

ョブではさらなる増加は難しく，実効性能をより実験と類似した壁遠方領域での第 2 ピークの出現高める必要性があると考えられる．

が明確に確認できる．ただしその値はまだ弱く，
レイノルズ数を上げることにより強い影響が出現

参考文献

するものと考えられる．
最後に本研究では，速度場に加えて温度場（パッシブスカラー場）の計算を同時に行っているが，

1) Kim, K. C., & Adrian, R. J.(1999), Very large-scale motion

本研究で行った時間積分長では，発達場が得られin the outer layer, Physics of Fluids, 11(2), 417-422.

ることができなかった．この原因は大規模構造の2)Hutchins, N., & Marusic, I. (2007), Evidence of very long

揺らぎの影響により，低レイノルズ数の場合より

meandering features in the logarithmic region of turbulent

もより長い時間積分長が必要になるためと考えらboundary layers, J. Fluid Mech., 579, 1-28.

れる．

3) Morinishi, Y., Lund, T. S, Vasilyev, O.V., and Moin, P.
(1998), Fully Conservative Higher Order Finite Difference
Schemes for Incompressible Flow, J. Comp. Phys. 143,
90-124.
4) Lozano-Durana, A. and Jimenez, J. (2014), Effect of thecomputational domain on direct simulations of turbulent
channels up to Reτ = 4200, Phys. Fluids, 26, 011702.
5) Bernardini, M., Pirozzoli, S., and Orlandi, P. (2014),
Velocity statistics in turbulent channel flow up to Ret =4000, J.
Fluid Mech., 742, 171-191.

図 7 主流方向 pre-multiplied スペクトルのコンター（Outer peak が第 2 ピークに相当）

6 まとめ
本研究では，京都大学学術情報メディアセンタ
のCray XE6/256 ノードを使用した大規模直接数値計算を実行した．その結果，実効演算速度で５京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 45

サービスの記録・報告システム A 運転状況 ( 2015 年 4 月 ～ 2015 年 9 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]障害発生日時2015/04/01

0:00

2015/04/02

9:00

33.00

2015/06/09

9:00

2015/06/10

9:30

24.50

2015/08/04

9:00

2015/08/05

9:30

24.50

サービス再開日時

ダウン時間[h]２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

687.0026,532

61,037

9,796,650

7,396,260

782.1

48

%

5月

744.0035,56650,344

8,287,300

7,049,130

751.9

47

%

6月

695.5044,51876,5699,282,730

7,696,580

751.8

54

%

7月

744.0030,534

63,981

11,240,500

9,893,430

752.0

63

%

8月

719.5030,07551,451

9,152,680

7,614,200

753.8

53

%

9月

720.0030,59159,057

9,958,610

7,842,360

751.8

59

%

計

4310.00

197,816

362,440

57,718,470

47,491,960

757.2

54

%

ノード
利用率

ノード利用率(棒グラフ)＋コア利用率(折れ線グラフ)
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

ノード利用率
コア利用率

4月

5月

6月

7月

8月

9月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

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

システム B 運転状況 ( 2015 年 4 月 ～ 2015 年 9 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]2015/04/01

0:00

2015/04/02

9:00

33.00

2015/06/09

9:00

2015/06/10

9:30

24.50

2015/08/04

9:00

2015/08/05

9:30

24.50

障害発生日時

サービス再開日時

ダウン時間[h]２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

687.0061,517

170,148

2,926,430

213,727

511.3

61

%

5月

744.0096,526

175,549

3,922,710

267,074

536.3

64

%

6月

695.5073,642

240,809

3,977,930

307,175

525.6

66

%

7月

744.0068,874

295,033

4,152,510

308,559

539.0

68

%

8月

719.50103,781257,457

3,599,500

292,702

536.3

64

%

9月

720.00102,741

287,304

3,734,530

296,901

538.9

65

%

計

4310.00

507,081

1,426,300

22,313,610

1,686,138

531.3

64

%

ノード
利用率

ノード利用率(棒グラフ)＋コア利用率(折れ線グラフ)
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

ノード利用率
コア利用率

4月

5月

6月

7月

8月

9月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

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

システム C 運転状況 ( 2015 年 4 月 ～ 2015 年 9 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]2015/04/01

0:00

2015/04/02

9:00

33.00

2015/06/09

9:00

2015/06/10

9:30

24.50

2015/08/04

9:00

2015/08/05

9:30

24.50

障害発生日時

サービス再開日時

ダウン時間[h]２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

687.00672

7,396

84,59755,498

12.3

47

%

5月

744.00

1,096

9,079

96,45559,520

12.0

48

%

6月

695.50

1,031

8,386

79,40544,914

12.0

37

%

7月

744.00

1,006

4,410

56,80334,965

12.0

26

%

8月

719.50762

1,655

45,38933,046

11.9

19

%

9月

720.00992

4,100

34,873

16,765

12.0

20

%

計

4310.00

5,559

35,027397,522244,707

12.0

33

%

ノード
利用率

ノード利用率(棒グラフ)＋コア利用率(折れ線グラフ)
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

ノード利用率
コア利用率

4月

5月

6月

7月

8月

9月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

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

システム D 運転状況 ( 2015 年 4 月 ～ 2015 年 9 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]2015/04/01

0:00

2015/04/02

9:00

33.00

2015/06/09

9:00

2015/06/10

9:30

24.50

2015/08/04

9:00

2015/08/05

9:30

24.50

障害発生日時2015/05/05

サービス再開日時

15:00

2015/05/05

ダウン時間[h]

22:45

7.75

２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

687.005,509

31,657

4,078,360

187,950

357.3

56

%

5月

736.25

10,062

41,0405,105,330

269,159

398.5

64

%

6月

695.5015,94258,826

5,465,200

310,326

415.8

66

%

7月

744.0027,41863,207

6,358,920

291,873

416.0

70

%

8月

719.5012,54673,9555,359,440

303,489

416.0

65

%

9月

720.0010,77869,7054,830,910

249,961

416.0

56

%

計

4302.25

82,255338,389

31,198,160

1,612,758

403.3

63

%

ノード
利用率

ノード利用率(棒グラフ)＋コア利用率(折れ線グラフ)
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

ノード利用率
コア利用率

4月

5月

6月

7月

8月

9月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

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

システム E 運転状況 ( 2015 年 4 月 ～ 2015 年 9 月 )１）保守作業に伴うサービス休止およびシステムダウン障害発生状況
保守作業に伴うサービス休止
保守開始日時システムダウン障害発生状況

サービス再開日時

保守時間[h]障害発生日時

サービス再開日時

ダウン時間[h]2015/04/01

0:00

2015/04/02

9:00

33.00

2015/04/07

9:00

2015/04/07

10:30

1.50

2015/06/09

9:00

2015/06/10

9:30

24.50

2015/04/09

8:40

2015/04/09

10:30

1.83

2015/08/04

9:00

2015/08/05

9:30

24.50

２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

683.67

2,545

7,802

2,046,000

59,226

363.7

76

%

5月

744.00

1,764

5,090

2,145,080

56,508

359.6

86

%

6月

695.50

2,012

5,712

2,200,800

75,703

360.0

81

%

7月

744.00

1,416

6,854

2,300,000

57,220

360.0

85

%

8月

719.50

2,410

9,395

1,969,270

126,659

359.8

76

%

9月

720.002,050

15,912

2,207,210

58,038

360.0

84

%

計

4306.67

12,19750,764

12,868,360

433,353

360.5

81

%

ノード
利用率

ノード利用率(棒グラフ)＋コア利用率(折れ線グラフ)
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

ノード利用率
コア利用率

4月

5月

6月

7月

8月

9月



占有時間 = 合計(経過時間×占有コア数)

平均稼動ノード数 = 電源 ON 状態のノード数の月平均 (10 分間隔のサンプリングデータより算出)



ノード利用率 = 稼動ノードに対するジョブが実行されているノードの割合

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

汎用コンピュータシステムのサービス状況
1. ホ ス テ ィ ン グ ・ ホ ー ム ペ ー ジ サ ー ビ ス 利 用 状 況

VM�����������
6���
55��

556�

554�

WEB�����������

562�

554�

552�

549�

������������

5���
45��
4���
�5��
����

282�

28��

27��

269�

268�

267�

25��
2���
�5��
����
5��
��

��

��

2�

2�

2�

2�

平成27年4⽉�

平成27年5⽉�

平成27年6⽉�

平成27年7⽉�

平成27年8⽉�

平成27年9⽉�

(平成27年4⽉から平成27年9⽉)

大型計算機システム利用承認件数について

平成27年9月末現在、大型計算機システムの利用件数は、2,782件となっています。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 51

サービスの記録・報告

センター利用による研究成果―平成 26 年度計算機利用結果報告書から―Computing System and Infrastructure，2015【1000 情報学】
1.

2.

Tatsuya Abe and Toshiyuki Maeda ： A【1500 生活科学】

General Model Checking Framework for

8.

Various Memory Consistency Models：IEEE

Takaaki Shimada, and Hiromi Teraoka ：28th International Parallel & Distributed

Prediction of Diaper Shape While Worn by

Processing Symposium Workshops，pp. 332 -Using Finite Element Method - Focus on341，2014

Posture of Open Legs Wearing Tape-Ty pe

Tatsuya Abe

andBaby Disposable Diaper：SEN-I GAKKAISHI，

Toshiyuki Maeda ：Vol.70，No.8，pp. 180 - 186，2014

Optimization of a General Model Checking

3.

4.

Framework for Various Memory Consistency【4300 物理学】

Models：Proceedings of the 8th International

9.

Charge-Transfer

Space Programming Models，2014

FMO-LCMO Approach: Hole Transfer in

Masatoshi Kawai, Takeshi Iwashita, andDNA

Hiroshi Nakashima：about 48hr to the review

Range-Separated DFT ： Chemical Physics

due ： Proc. Intl. Mtng. High-Performance

Letters，Vol.621，pp. 96 - 101，2015

7.

MatrixwithElements

Parameter

by

Tuned

Computing for Computational Science，2014

10. Yosuke Harashima, Kiyoyuki Terakura, Hiori

Takeshi Minami, Motoharu Hibino, Tasuku

Kino, Shoji Ishibashi and Takashi Miyake：

Hiraishi, Takeshi Iwashita, and Hiroshi

First-Principles Study of Structural andNakashima：Automatic Parameter Tuning ofMagneticThree-Dimensional Tiled FDTD Kernel：Proc.

R(Fe,Ti)12N (R=Nd, Sm, Y)：JPS Conf. Proc.，Intl. WS. Automatic Performance Tuning，Vol.5，pp. 011021-1 - 011021-8，2015

Properties ofR(Fe,Ti)12and

11. Ryuma Sato, Hirotaka Kitoh-Nishioka, Koji

Shingo Okuno, Tasuku Hiraishi, Hiroshi

Ando, Takahisa Yamato ： Computational

Nakashima, Masahiro Yasugi, and Jun

study on the roles of amino acid residues in

Sese：Parallelization of Extracting Connected

the active site formation mechanism ofSubgraphs with Common Itemsets ： IPSJ

blue-light photore：Chemical Physics Letters，Trans. Programing，Vol.7，No.3，pp. 22 - 39，

Vol.633，pp. 247 - 251，2015

2014
6.

Hirotaka Kitoh-Nishioka and Koji Ando：

Conference on Partitioned Global Address

2014
5.

Yosuke Horiba, Shigeru Inui, Yuki Maeda,

12. Shun Sakuraba, Kiyoshi Asai, Tomoshi

Takeshi Iwashita, Naokazu Takemura, andKamedaHiroshi Nakashima：A Fill-In Strategy for

Dimerization Free-energy Changes Upon

Fast ICCG Solver with SIMD Vectorization：

Mutations

Annual Meeting on Advanced Computing

Simulations ： The Journal of Physical

System and Infrastructure，2015

Chemistry Letters，Vol.6，pp. 4348 - 4351，

Shingo Okuno, Tasuku Hiraishi, Hiroshi

2015

：

Predicting
Using

RNA

Molecular

DuplexDynamics

Nakashima, Masahiro Yasugi, and Jun【4400 地球惑星科学】

Sese ： Reducing Redundant Search using

13. Kazem Ardaneh, Dongsheng Cai, Ken-Ichi

Exception Handling in a Task-Parallel

Nishikawa,

Language ： Annual Meeting on Advanced

COLLISIONLESS WEIBEL SHOCKS AND

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

andBertrand

Lembège ：IN

21. 森田早紀，西川雅章，北條正樹：有限要素解析

GAMMA-RAY BURSTS：The Astrophysicalを用いた角部を有する構造のための織物プリフ

Journal，Vol.811，pp. 57 (9pp)，2015

ォーム賦形における変形の評価：第５回日本複ELECTRONACCELERATION

合材料会議 講演論文集，pp. 3B01 -，2014【4500 プラズマ科学】
14. Kevin Obrejan, Kenji Imadera, Jiquan Li and22. 内藤悠太，西川雅章，北條正樹：形状記憶複合

Yasuaki Kishimoto：Development of a global

材料の展開のための加熱温度分布の検討：第５toroidal gyrokinetic Vlasov code with new

回日本複合材料会議 講演論文集，pp. 3C-09 -，real space field solver：Plasma and Fusion

2014
23. 西川雅章，上田桂，北條正樹，武田一朗：キンResearch，2014
15. Kenji Imadera, Yasuaki Kishimoto, Kevin

ク損傷予測に基づく一方向繊維強化複合材料のObrejan Takuya Kobiki and Jiquan Li：Global

圧縮強度解析法の検討：2014 年度 JCOM 若

Profile Relaxation Coupled with E × B

手シンポジウム 講演論文集，pp.27 -，2014

ITG

24. 橋爪良輔，新玉重貴，西川雅章，北條正樹：ひTurbulence ： Proceedings of 25th Fusion

ずみ不変量破壊理論を用いた CFRPL 字曲り梁Energy Conference，2014の層間はく離強度解析：第 39 回複合材料シン

Staircasein

ToroidalFlux-Driven

ポジウム 講演論文集，pp.79 -，2014【4600 基礎化学】
16. 河合信之輔：多自由度系ダイナミクスにおける25. 北口尚紀，金崎真人，西川雅章，北條正樹：熱

隠された自由度の抽出：日本物理学会第 70 回融着を目的とした CF/PA6 複合材料の平行平板

年次大会，2015間における繊維・樹脂流動特性評価：第 39 回

17. 河合信之輔，田中明徳，久邇実加子，関根理香：水溶液中のイオンの拡散挙動における環境ダイ
ナミクス：日本化学会第 95 春季年会，2015

複合材料シンポジウム 講演論文集，pp.65 -，2014
26. 福造晃，西川雅章，北條正樹：射出成形された【4800 材料化学】

CFRTP の非線形力学特性の解析：第 39 回複合

18. K. Shizu, M. Uejima, H. Nomura, T. Sato, K.

材料シンポジウム 講演論文集，pp. 81 -，2014

Tanaka, H. Kaji, and C. Adachi：Enhanced

27. 森田早紀，西川雅章，北條正樹：角部を有する

Electroluminescence from a Thermally構造のための織物プリフォーム賦形における繊Activated Delayed-Fluorescence Emitterby Suppressing Nonradiative Decay：Phys.

維束変形の評価：第 39 回複合材料シンポジウ

Rev. Applied，Vol.3，No.1，pp. 014001-1 -28. 内藤悠太，西川雅章，北條正樹：展開時間を考014001-7，2015

ム 講演論文集，pp. 193 -，2014
慮した形状記憶複合材料の展開モデルの検討：

19. K. Shizu, H. Tanaka, M. Uejima, T. Sato, K.
Tanaka, H. Kaji, and C. Adachi：Strategy for

第 58 回宇宙科学技術連合講演会 講演論文集，pp. 2M07 -，2014

Designing Electron Donors for Thermally

29. 永田恵輔，西川雅章，北條正樹：シェル要素をActivated Delayed Fluorescence Emitters：J.

用いた展開膜の展開挙動の解析：第 58 回宇宙Phys. Chem. C，Vol.119，No.3，pp. 1291 - 1297，

科学技術連合講演会 講演論文集，pp. 2M08 -，

2015

2014

【5000 機械工学】【5200 土木工学】

20. T.Kitano, J. Nishio, R. Kurose, S. Komori：

30. Taiki

Fumuro,Seizo

Tanaka,Kazuo

Effect of ambient pressure on soot formation

Kashiyama ： Stabilized Finite Elementin oxy-fuel spray jet flame：In Proc. of The

Method15th
International
Heat
Transfer
Conference (IHTC-15) ， Vol. DIGITAL

Simulation：The 1st International Conference
on Computational Engineering and Science

LIBRARY，pp. IHTC15-8744 (15 pages) -，for Safety and Environmental Problems，pp.

2014

210 - 213，2014

forLarge-Scale

Tsunami京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 53

31. 不室太希，田中聖三，川辺赳史，樫山和男：安Sciencein

Space(ISPS6)

Conference

and

on

10th

定化有限要素法による大規模津波シミュレーシ

International

Two-Phase

ョンとその可視化：計算工学講演会論文集，

Systems for Space and Ground Applications

2014(ITTW2015)，pp. 16FlMa-13 -，2015

32. 不室太希，田中聖三，樫山和男：VOF 安定化有40. 山本卓也、高木洋平、岡野泰則：円形液膜内振限要素法による大規模並列津波解析：第 28 回

動 Marangoni 対流に対する数値解析：日本流

数値流体力学シンポジウム 講演要旨集，2014

体力学会年会 2015 講演要旨集，2015

33. Yasuaki Ito, Hiromichi Shirato, Masaru【5600 総合工学】

Matsumoto, Kojiro Nozawa ： Span-wise

41. Kazushi Sato and Katsuyoshi Fukiba：Effectcorrelations of fluctuating lift forces on 2D

of Wing Shape and Reynolds Number in Drag

rectangular cylinders：Journal of Wind andReduction Method by Using Wing withEngineering，Vol.11，No.2，pp. 1 - 10，2014

Rotation

34. Yasuaki Ito, Hiromichi Shirato, Masaru

Shaft ： 11th

International

Conference on Flow Dynamics，2014

span-wise

42. 佐藤一志，吹場活佳：自由回転軸を持つ翼型を

correlations of fluctuating lift forces on

用いた抗力低減法における翼形状と Re 数の影

rectangular cylinders with various side

響：第 46 回流体力学講演会，2014

Matsumoto ： Researchonratio ： Proceedings of 6th Computational【6200 林学】

Wind Engineering，2014

43. 渡辺靖崇・鈴木保志・後藤純一・酒井寿夫：施【5400 材料工学】

業方法の違いによる人工林における土砂流出量

35. Y. Hinuma, A. Gruneis, G. Kresse, and F.

の変化：森林利用学会誌，Vol.29，No.3，pp. 149

Oba：Band alignment of semiconductors from- 156，2014

density-functional theory and many-body【6800 薬学】

perturbation theory：Physical Review B，

44. Nobuya Tsuji, Yusuke Kobayashi, Yoshiji

Vol.90，No.15，pp. 155405-1 - 155405-16，2014

Takemoto：Electrophilic iodine(I) compounds

36. Y. Kumagai, M. Choi, Y. Nose, and F. Oba：

induced semipinacol rearrangement via C-X

First-principles study of point defects inbond cleavage：Chem. Commun.，Vol.50，pp.

chalcopyrite ZnSnP2：Physical Review B，Vol.90，No.12，pp. 125202-1 - 125202-12，2014
37. Y. Kumagai and F. Oba：Electrostatics-based

13691 - 13694，2014
45. Chihiro Tsukano, Masataka Okuno, Hiromi
Nishiguchi,

YoshijiTakemoto

finite-size corrections for first-principles point

Palladium(0)-Catalyzed

defect calculations：Physical Review B，Vol.89，

Carbon[BOND]Hydrogen

No.19，pp. 195205-1 - 195205-15，2014

Functionalization

for

：
Bond

the

Synthesisof【5500 プロセス工学】

Indoloquinazolinones：Adv. Synth. Catal.，

38. T. Yamamoto, Y. Takagi, Y. Okano and S.

Vol.356，No.7，pp. 1533 - 1538，2014

Dost：Numerical investigation of the effect of

46. Yusuke

Kobayashi,

ShanjiLi,

Yoshijifree surface shape on the direction of

Takemoto

thermocapillary flow in a thin circular pool：

Hydrogen-Bond-Donating Organocatalyst for

J. Chem. Eng. Jpn，Vol.48，pp. 407 - 417，2015

Enantioselective Nucleophilic Epoxidation of

39. Takuya Yamamoto, Youhei Takagi, Yasunori

α,β-Unsaturated Amides：Asian. J. Org.

Okano, Sadik Dost：Liquid film volume effect：Oxidant-ResistantChem，Vol.3，No.4，pp. 403 - 407，2014

on unsteady thermocapillary flow due to【9999 その他】

hydrothermal wave in a thin circular liquid

47. X.-Y. Jin, Y. Kuramashi, Y. Nakamura, S.

film：Proc. of The joint conference of 6th

Takeda, A. Ukawa：Critical endpoint of theInternational

finite temperature phase transition for three

Symposiumon

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

flavor QCD：Physical Review D，Vol.91，pp.
014508 -，2015

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 55

資料

大型計算機システム利用負担金
別表１ スーパーコンピュータシステム
コース

タイプ

セット

タイプA
タイプB
タイプC
パーソナル タイプD
タイプE

基本
基本
基本

エントリ

タイプG
タイプA1
タイプA2
タイプA3
タイプB1
タイプB2
タイプB3
タイプC1
グループ タイプC2

基本
基本
基本

タイプA
タイプB

基本

100,000 円/年
200,000 円/年
200,000 円/年
240,000 円/年
120,000 円/年
600,000 円/年
300,000 円/年
250,000 円/年
250,000 円/年
300,000 円/年
150,000 円/年
750,000 円/年
375,000 円/年
400,000 円/年
200,000 円/年
240,000 円/年
120,000 円/年
300,000 円/年
150,000 円/年
360,000 円/年
90,000 円/年
900,000 円/年
450,000 円/年
280,000 円/年
140,000 円/年
336,000 円/年
84,000 円/年
840,000 円/年
420,000 円/年
250,000 円/年
250,000 円/年
20,000 円/週(7日）
5,000 円/週(7日）
24,000 円/週(7日）
6,000 円/週(7日）
20,000 円/週(7日）
10,000 円/週(7日）
30,000 円/週(7日）
7,500 円/週(7日）
28,000 円/週(7日）
7,000 円/週(7日）
24,000 円/週(7日）
12,000 円/週(7日）
750,000 円/年
375,000 円/年
20,000 円/年

追加単位

最小
追加単位

最小
追加単位

最小
追加単位

最小
追加単位

最小
追加単位

最小
追加単位

最小

最小
追加単位

最小
追加単位

最小

追加単位
最小
タイプC
追加単位
大規模ジョブ
最小
タイプD
追加単位
最小
タイプE
追加単位
最小
タイプG
追加単位

専用クラスタ―ライセンスサービス

12,600 円/年
100,000 円/年
100,000 円/年
100,000 円/年
100,000 円/年
100,000 円/年

最小

追加単位
最小
タイプD1
追加単位
最小
タイプD2
追加単位
最小
タイプD3
追加単位
最小
タイプE1
追加単位
最小
タイプE2
追加単位
最小
タイプE3
追加単位

タイプG1

利用負担額

最小
追加単位

提供サービス
経過時間
（時間）

ディスク
（GB)

利用者
番号

最大1ノード相当((16コア、64GBメモリ)×1)
最大4ノード相当((32コア、64GBメモリ)×4)
最大4ノード相当((16コア、64GBメモリ)×4)
最大2ソケット相当((8コア、384GBメモリ)×2)
最大2ノード相当((28コア、64GBメモリ)×2)
最大2ノード相当((10コア、32GBメモリ + 1MIC)×2)

1
168
168
168
168
168

60
1,000
1,000
1,000
1,000
1,000

-

共有 最大2ノード相当((16コア、64GBメモリ + 1GPU)×2)

168

1,000

-

4ノード((32コア、64GBメモリ)×4)
336
4ノード((32コア、64GBメモリ)×4)
8ノード((32コア、64GBメモリ)×8)
336
A
準優先
4ノード((32コア、64GBメモリ)×4)
8ノード((32コア、64GBメモリ)×8)
336
占有
4ノード((32コア、64GBメモリ)×4)
4ノード((16コア、64GBメモリ)×4)
336
優先
4ノード((16コア、64GBメモリ)×4)
8ノード((16コア、64GBメモリ)×8)
336
B
準優先
4ノード((16コア、64GBメモリ)×4)
8ノード((16コア、64GBメモリ)×8)
336
占有
4ノード((16コア、64GBメモリ)×4)
4ソケット((8コア、384GBメモリ)×4)
336
優先
2ソケット((8コア、384GBメモリ)×2)
C
4ソケット((8コア、384GBメモリ)×4)
336
準優先
2ソケット((8コア、384GBメモリ)×2)
4ノード((28コア、64GBメモリ)×4)336
優先
2ノード((28コア、64GBメモリ)×2)
8ノード((28コア、64GBメモリ)×8)
336
D
準優先
2ノード((28コア、64GBメモリ)×2)
8ノード((28コア、64GBメモリ)×8)
336
占有
4ノード((28コア、64GBメモリ)×4)
4ノード((10コア、32GBメモリ + 1MIC)×4)336
優先2ノード((10コア、32GBメモリ + 1MIC)×2)
8ノード((10コア、32GBメモリ + 1MIC)×8)
336
E
準優先
2ノード((10コア、32GBメモリ + 1MIC)×2)
8ノード((10コア、32GBメモリ + 1MIC)×8)
336
占有
4ノード((10コア、32GBメモリ + 1MIC)×4)
2ノード((16コア、64GBメモリ + 1GPU)×2)
336
B
優先
(GPU)
2ノード((16コア、64GBメモリ + 1GPU)×2)
8ノード((32コア、64GBメモリ)×8)A
占有
2ノード((32コア、64GBメモリ)×2)
8ノード((16コア、64GBメモリ)×8)
B
占有
2ノード((16コア、64GBメモリ)×2)
4ソケット((8コア、384GBメモリ)×4)
C
占有
2ソケット((8コア、384GBメモリ)×2)
8ノード((28コア、64GBメモリ)×8)
D
占有
2ノード((28コア、64GBメモリ)×2)
8ノード((10コア、32GBメモリ + 1MIC)×8)
E
占有
2ノード((10コア、32GBメモリ + 1MIC)×2)
4ノード((16コア、64GBメモリ+
1GPU)×4)B
占有(GPU)
2ノード((16コア、64GBメモリ + 1GPU)×2)
8ノード((16コア、64GBメモリ)×8)
B
4ノード((16コア、64GBメモリ)×4)
可視化ソフト（AVS,ENVI/IDL)およびプリポストウェアの1ライセンスにつき

8,000
8,000
9,600
4,800
16,000
8,000
8,000
8,000
9,600
4,800
16,000
8,000
8,000
4,000
4,800
2,400
8,000
4,000
9,600
2,400
16,000
8,000
8,000
4,000
9,600
2,400
16,000
8,000
4,000
4,000
16,000
8,000

8
8
16
8
16
8
8
8
16
8
16
8
16
8
16
8
8
4
16
4
16
8
8
4
16
4
16
8
8
8
16
8

システム資源システム バッチ
B
A
B
C
D
E
B
(GPU)

共有
共有
共有
共有
共有
共有

優先

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。
2． 大型計算機システムの全ての利用者は、上記表のサービスの他、次のサービスを受けることができる。
1） 大判プリンタサービス
2） その他、大型計算機システムが提供するサービス、機器の利用
3． 上記表の大規模ジョブコース、ライセンスサービスの申請には、大型計算機システムの利用者であることが必要である。
4． 「共有」： 当該カテゴリのユーザ間で一定の計算資源を共有するベストエフォートのスケジューリングを行う。
「準優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/4の計算資源が確保されることを保証する。
「優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/2の計算資源が確保されることを保証する。
「占有」： 稼働状況によらず記載値（以上）の計算資源が確保されることを保証する。
5． ディスク容量はバックアップ領域（最大で総容量の1/2）を含む。

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

6． グループコース及び専用クラスタコースのシステム資源は、下記の負担額を支払うことにより増量することができる。
なお増量は各月１日に実施し、増量した資源は当該年度末までの期間にわたって利用されるものとする。
コース

システム資源増量単位

タイプA120,000 円/月

4ノード((32コア、64GBメモリ)×4)

8,000

タイプA2

12,000 円/月

4ノード((32コア、64GBメモリ)×4)

4,800

タイプA3

30,000 円/月

4ノード((32コア、64GBメモリ)×4)

8,000

タイプB1

25,000 円/月

4ノード((16コア、64GBメモリ)×4)

8,000

タイプB2

15,000 円/月

4ノード((16コア、64GBメモリ)×4)

4,800

タイプB3

37,500 円/月
20,000 円/月
12,000 円/月

4ノード((16コア、64GBメモリ)×4)
2ソケット((8コア、384GBメモリ)×2)
2ソケット((8コア、384GBメモリ)×2)
2ノード((28コア、64GBメモリ)×2)
2ノード((28コア、64GBメモリ)×2)
4ノード((28コア、64GBメモリ)×4)
2ノード((10コア、32GBメモリ + 1MIC)×2)
2ノード((10コア、32GBメモリ + 1MIC)×2)
4ノード((10コア、32GBメモリ + 1MIC)×4)
2ノード((16コア、64GBメモリ + 1GPU)×2)
4ノード((16コア、64GBメモリ)×4)

8,000
4,000
2,400

グループ タイプC1
タイプC2
タイプD1
タイプD2
タイプD3
タイプE1
タイプE2
タイプE3
タイプG1

専用クラスタ

ディスク増量（ＧＢ）

追加負担金額（増量単位あたり）タイプ―

15,000 円/月
9,000 円/月
45,000 円/月
14,000 円/月
8,400 円/月
42,000 円/月

25,000 円/月
37,500 円/月

4,000
2,400
8,000
4,000
2,400
8,000

4,000
8,000

7． グループコース及び専用クラスタコースを通年でなく利用する場合には、下記の負担額を支払うものとする。
ただし、利用期間は当該年度内に限るものとする。
利用期間
最小
追加単位
最小
タイプA2
追加単位
最小
タイプA3
追加単位
最小
タイプB1
追加単位
最小
タイプB2
追加単位
最小
タイプB3
追加単位
最小
タイプC1
追加単位
グループ
最小
タイプC2
コース
追加単位
最小
タイプD1
追加単位
最小
タイプD2
追加単位
最小
タイプD3
追加単位
最小
タイプE1
追加単位
最小
タイプE2
追加単位
最小
タイプE3
追加単位
最小
タイプG1
追加単位
最小
専用クラス
―
タ コース
追加単位

タイプA1

3ヶ月
80,000 円
80,000 円
96,000 円
48,000 円
240,000 円
120,000 円
100,000 円
100,000 円
120,000 円
60,000 円
300,000 円
150,000 円
160,000 円
80,000 円
96,000 円
48,000 円

6ヶ月
120,000 円
120,000 円
144,000 円
72,000 円
360,000 円
180,000 円
150,000 円
150,000 円
180,000 円
90,000 円
450,000 円
225,000 円
240,000 円
120,000 円
144,000 円
72,000 円

9ヶ月
180,000 円
180,000 円
216,000 円
108,000 円
540,000 円
270,000 円
225,000 円
225,000 円
270,000 円
135,000 円
675,000 円
337,500 円
360,000 円
180,000 円
216,000 円
108,000 円

75,000 円
37,500 円
90,000 円
22,500 円
225,000 円
112,500 円
70,000 円
35,000 円
84,000 円
21,000 円
210,000 円
105,000 円

150,000 円
75,000 円
180,000 円
45,000 円
450,000 円
225,000 円
140,000 円
70,000 円
168,000 円
42,000 円
420,000 円
210,000 円

225,000 円
112,500 円
270,000 円
67,500 円
675,000 円
337,500 円
210,000 円
105,000 円
252,000 円
63,000 円
630,000 円
315,000 円

100,000 円
100,000 円
300,000 円
150,000 円

150,000 円
150,000 円
450,000 円
225,000 円

225,000 円
225,000 円
675,000 円
337,500 円

8． グループコース及び専用クラスタコースの利用者番号は利用者あたり年額5,000円を負担することで追加できる。
9 ． 機関・部局定額制度
他機関又は学内における部局（『国立大学法人京都大学の組織に関する規程』第3章第2節から第11節で定める組織をいう。）の組織が、その組織単位でグループコースサービス（年間）の利用を申請する場合、料金表(年間）に掲載額の1.5倍を利用負担金とする。なお、利用負担金額が150万円未満の場合は100人、150万円を超える場合は、150万円毎に100人までの利用者を認める。
10 ． スパコン連携サービス
学術情報メディアセンターのスーパーコンピュータシステムと密な連携により、学内における部局の組織が計算サーバ等を設置する
場合、下記の負担額を支払うものとする。
冷却方式
水冷
空冷

利用負担額
10,300 円/月
12,900 円/月

利用負担額算定単位
水冷冷却方式の計算サーバ等の定格電力 1kWにつき
空冷冷却方式の計算サーバ等の定格電力 1kWにつき

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 57

別表２（汎用コンピュータシステム）
区分
VMホスティングサービス
ホームページサービス
ストリーミングサービス

利用負担額
72,000円／年
6,000円／年
6,000円／年

単位
1仮想マシンにつき1ドメイン名につき1申請につき

備考
１．利用負担額は、総額表示である。
２．上記表の汎用コンピュータシステムのサービスを利用するためには、大型計算機システムの利用者であることが必要である。
３．VM ホスティングサービスにおいて、下記の負担額を支払うことにより CPU、メモリ、ディスクを増量することがで
きる。
区分
CPU増量
メモリ増量
ディスク増量

利用負担額
18,000円／年
18,000円／年
18,000円／年

単位
2コアにつき（最大8コアまで）
8GBにつき（最大64GBまで）
200GBにつき（最大1,000GBまで）

４．VM ホスティングサービスにおいて VMware を用いる場合は、下記の負担額を支払うことにより VMware の利用及びCPU、メモリ、ディスクを増量することができる。ただし、システム資源が限られているためサービスの提供を限定することがある。
区分
VMware利用
CPU増量
メモリ増量
ディスク増量

利用負担額
72,000円／年
36,000円／年
36,000円／年
18,000円／年

単位
1仮想マシンにつき2コアにつき（最大8コアまで）
8GBにつき（最大64GBまで）
200GBにつき（最大1,000GBまで）

５．ホームページサービス及びストリーミングサービスにおいて、下記の負担額を支払うことにより公開スペースの上限を拡大することができる。
区分
公開スペース20GBプラン
公開スペース50GBプラン

利用負担額
3,000円／年
9,000円／年

６．利用負担額は、当該年度（4月から翌年3月まで）の利用に対して年額として算定するが、年度途中から利用を開始する場合には月数に応じて減額する。

別表３ スーパーコンピュータシステム（民間機関利用）
システム資源

システム

A

B

D

E

8ノード(32コア、64GBメモリ)×8)
12ノード(32コア、64GBメモリ)×12)
16ノード(32コア、64GBメモリ)×16)
8ノード(16コア、64GBメモリ)×8)
12ノード(16コア、64GBメモリ)×12)
16ノード(16コア、64GBメモリ)×16)
8ノード(28コア、64GBメモリ)×8)
12ノード(28コア、64GBメモリ)×12)
16ノード(28コア、64GBメモリ)×16)
8ノード(10コア、32GBメモリ + 1MIC)×8)
12ノード(10コア、32GBメモリ + 1MIC)×12)
16ノード(10コア、32GBメモリ + 1MIC)×16)

経過時間 ディスク
（時間）
（GB)

336
336
336
336
336
336
336
336
336
336
336
336

利用者
番号

9,600
14,40019,2009,600
14,40019,2009,600
14,40019,2009,600
14,40019,200

16
24
32
16
24
32
16
24
32
16
24
32

利用負担額
960,000 円/年
1,440,000 円/年
1,920,000 円/年
1,200,000 円/年
1,800,000 円/年
2,400,000 円/年
1,440,000 円/年
2,160,000 円/年
2,880,000 円/年
1,344,000 円/年
2,016,000 円/年
2,688,000 円/年

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。
2． ディスク容量はバックアップ領域（最大で総容量の1/2）を含む。
3． 通年でなく利用する場合には、下記の負担額を支払うものとする。
ただし、利用期間は当該年度内に限るものとする。
システム

A

B

D

E

システム資源
8ノード
12ノード
16ノード
8ノード
12ノード
16ノード
8ノード
12ノード
16ノード
8ノード
12ノード
16ノード

3ヶ月
240,000 円
360,000 円
480,000 円
300,000 円
450,000 円
600,000 円

利用期間
6ヶ月
480,000 円
720,000 円
960,000 円
600,000 円
900,000 円
1,200,000 円

9ヶ月
720,000 円
1,080,000 円
1,440,000 円
900,000 円
1,350,000 円
1,800,000 円

360,000 円
540,000 円
720,000 円
336,000 円
504,000 円
672,000 円

720,000 円
1,080,000 円
1,440,000 円
672,000 円
1,008,000 円
1,344,000 円

1,080,000 円
1,620,000 円
2,160,000 円
1,008,000 円
1,512,000 円
2,016,000 円

58 − 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



システム D：magnolia.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



ホームページ・ホスティングサービス関する問い合わせ先【クラウドコンピューティング掛】
E-mail : whs-qa@media.kyoto-u.ac.jp / Tel : 075-753-7494
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/whs/

コンテンツ作成支援サービスに関する問い合わせ先【コンテンツ作成室】
E-mail : cpt@media.kyoto-u.ac.jp / Tel : 075-753-9012
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/content/

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.14 No.2（ Mar. 2016）− 59

編集後記

編 集 後 記情報環境機構のウェブページ
http://www.iimc.kyoto-u.ac.jp/ja/services/comp/supercomputer/coming_2016.html
等でもアナウンスされている通り，現行のスーパーコンピュータシステムのサービス（一部除く）は2016 年 8 月中に終了し，10 月より次期システムによるサービスを開始する予定です．本広報誌でも新システムの紹介記事を予定しておりますのでご期待ください．それに合わせてというわけでもないのですが，表紙などのデザインの刷新も計画しております．その際，この編集後記を残すのかということをちょっと迷っております．編集委員の日記的な文章という意味での編集後記は今回か次回くらいが最後になるかもしれません．楽しく読んでいるので無くさないでくれという意見がありましたらご一報を．
(副部会長)

京都大学学術情報メディアセンター全国共同利用版広報 Vol．14，No．2
2016 年 3 月 14 日発行
広報編集部会
編集者 京都大学学術情報メディアセンター

深沢 圭一郎（部会長）

全国共同利用版広報編集部会

平石 拓（副部会長）

発行者 〒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
株式会社エヌジーピー

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

表紙デザイン：谷 卓司（ティアンドティ･デザインラボ）

Vol.14, No.2 2015

【巻頭言】
・Vol.14，No.2 号の発刊にあたって ......................................................................................................................................................................................... 1
牛島 省【スーパーコンピュータ共同研究制度（若手研究者奨励枠）
研究報告】
・三次元非定常マランゴニ対流の解明 .................................................................................................................................................................................. 2
山本 卓也
・高プラントル数流体を用いた HZ 液柱内温度差マランゴニ対流の数値解析 ............................................................................................... 4茂木 孝介
・準周期軌道に纏わる各種平均量の計算の高速化 ......................................................................................................................................................... 6
Suddhasattwa Das, Yoshitaka Saiki, Evelyn Sander, James A Yorke・フラグメント分子軌道法と非経験的に最適化した長距離補正密度汎関数法による電荷移動パラメータの高精度計算 .. 8
鬼頭（西岡 ) 宏任，安藤 耕司
・射影された時系列データに対する経験的確率微分方程式モデリングとアンサンブル軌道予測.................................................. 10中野 直人
・相分離による自己組織化構造形成に与える乱流影響 ............................................................................................................................................ 12高木 洋平
・液柱内の温度差と濃度差に起因するマランゴニ対流の不安定性の解明.................................................................................................... 14
水口 尚
・多孔質壁面による乱流構造の変化と抵抗低減効果 ................................................................................................................................................. 16
中本 真義
・有機太陽電池応用を目指した新規光機能性有機分子材料の構造と電子構造の解明........................................................................... 18
梅山 有和
・急縮小・急拡大流路における高分子溶融体流れのマルチスケールシミュレーション...................................................................... 20
原田 浩平，谷口 貴志
・周囲気体を考慮した HZ 液柱内温度差マランゴニ対流における粒子集合現象の数値解析 .............................................................. 23
中村 翔一【プログラム高度化支援事業研究報告】
・ディレクティブベースプログラミングによる FIT の高速化と性能検証 ...................................................................................................... 25
中畑 和之
・動的／静的水～土骨格連成有限変形解析コードの高度化～領域分割法の適用の試み～................................................................. 31
野田 利弘
・破壊力学に基づく損傷モデルを用いた鉄筋コンクリートの 3 次元破壊シミュレーション............................................................. 35
車谷 麻緒
・都市の人口規模・空間分布における秩序形成 ........................................................................................................................................................... 39
森 知也【スーパーコンピュータ共同研究制度（大規模計算支援枠）
研究報告】
・高次精度差分法による高レイノルズ数チャンネル乱流場の大規模直接数値シミュレーション.................................................. 42
山本 義暢
【サービスの記録・報告】
・スーパーコンピュータシステムの稼働状況とサービスの利用状況.............................................................................................................. 46・センター利用による研究成果（平成 26 年度）............................................................................................................................................................... 52【資料】
・大型計算機システム利用負担金 別表 .............................................................................................................................................................................. 56
・サービス利用のための資料一覧 ......................................................................................................................................................................................... 59【編集後記】
・編集後記、奥付 ............................................................................................................................................................................................................................... 60(March,2 16)