巻頭言

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

本号では、学術情報メディアセンターの平成 27 年度共同研究報告について特集いたします．平成 27 年度では若手研究者の支援事業に関して 13 件，プログラム高度化支援事業について 5 件，大規模計算の支援事業について 2 件の共同研究が実施され，それぞれの報告が掲載されています．数値計算に関する様々な報告が掲載されていますので，計算機を利用される研究者や学生の方の良い参考になればと思います．
採択件数の最も多い若手研究者の支援事業は，学生を含む 40 歳未満の若手研究者を対象とした利用者向けの奨励研究制度で，優れた研究提案でスーパーコンピュータを利用することで学術的にインパクトがある成果を創出できると期待される課題に対して計算機利用負担金の全部または一部を本センターが補助しています．平成 27 年度 2 次募集から若手に加えて，女性枠（年齢は問わない）も追加され，若手・女性研究者奨励研究と変更になりました。今回報告のあった 13 件には化学物性や流体計算，数理モデル，運動モデル，物理など多様な研究分野が含まれています．また，採択された研究者の所属は北海道大学，一橋大学，東京理科大学，名古屋大学，京都大学、大阪大学と京都大学だけでは無く全国の若手研究者との共同研究となっているのも特徴です．これまでスーパーコンピュータを利用されてこなかった分野の研究もあり，挑戦的な研究の支援の一つになってきたと感じています．
次に，プログラム高度化支援事業では，スーパーコンピュータをグループコースまたは専用クラスタコースでご利用の研究グループを対象に，ユーザから提供いただいたプログラ
ムに実際に手を加えて，制御構造・データ構造の改良による性能チューニングや並列化手法の改良による高度化・高性能化，問題分割・負荷分散方式などの改良による高度化・高性能化といった大規模計算プログラムの高度化・高性能化を支援しています．平成 27 年度では，
流体計算、コンクリート・地面の構造解析，経済問題といった研究分野のアプリケーション
高度化を支援させていただきました．プログラムによって，高度化の内容は異なりますが，
主に高並列計算が可能になるようにスレッド並列，MPI 並列の実装，ノード間通信の効率化などが行われ，同様の計算を行う研究者にとっても良い情報となると思います．
そして，大規模計算の支援事業ですが，スーパーコンピュータをグループコースまたは専用クラスタ コースでご利用の研究グループを対象に最大で 128 ノード×2 週間の大規模ジョブコース利用支援を行っています．今回は，地震の前兆現象の可能性のある電磁プラズマ波の大規模解析とプログラム高度化にも採択された流体計算の大規模実行への支援報告となっています．近年のスーパーコンピュータは大規模計算を行えないとその性能の恩恵を受けることができないため，このような大規模計算の情報を知っておくことは，今後大規模京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 1

化を狙う研究者にも役立つと思います．
また，本号では昨年に Salt Lake City で開催された世界最大のスーパーコンピュータの会議である SC16 への参加報告があります．本センターでは例年この会議にブースを出しており，情報の収集，広報活動も行っています．SC17 でもブースを出す予定ですので，ご興味のある方は本センターへお知らせください．
最後に，今回報告にありますような共同研究、計算機支援を主として，これからも皆様の研究，教育にご活用いただけるようにセンター教職員も尽力していきますので，今後ともご利用，ご支援のほど、よろしくお願いいたします．

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

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

液膜内振動マランゴニ対流のメカニズム解明
山本 卓也
大阪大学大学院基礎工学研究科物質創成専攻化学工学領域

宇宙飛行士の Donald Pettit 博士によって国際宇宙ステーション内で行われた円形液膜におけるマランゴニ対流に関する実験が行われた。この一部の実験において振動するマランゴニ対流が発生したが、その発生メカニズムや発生条件等は詳細には分かっておらず、その振動流の発生条
件とメカニズムを解明することとした。数値シミュレーションによってそのメカニズム解明を試み、その振動流に対する液膜変形度の影響が非常に大きいことが判った。さらに、その流動メカ
ニズムに対しても考察を行った。

1 緒言

2 数値解析手法

マランゴニ対流は表面張力勾配によって発生す数値シミュレーションに用いた解析領域の概略

る流れであり、バルク半導体結晶生成プロセスや図を Fig. 1 に示す。
計算領域は円形の液膜であり、製膜プロセスにおいてよく見られ、最終製品の品

流体領域のみを数値的に解いた。液膜は水である質に影響する[1]。このように、数多くの応用面でとし、Prandtl 数は 6.96 で詳細な物性値は既報[7]も見られるが、一般には浮力対流と共存するための通りである。

マランゴニ対流の影響だけを実験的に調査するの支配方程式はNavier-Stokes 式、
連続式、
energyは非常に難しい。このため、数多くの宇宙環境を式である。気液界面形状は動的に変形しないと仮利用した研究が行われてきた[2,3]。

定し、振動流に対する静的な界面形状の影響を調

宇宙環境を利用したマランゴニ対流実験の一つ査した。境界条件として、気液界面に沿ってに宇宙飛行士の Donald Pettit 博士の実験[4,5]がMarangoni 対流条件、断熱条件を課した。また、ある。これらの実験では円形リング内に液膜を作

半田ごての加熱を側壁での温度境界条件としてモ製し、この液膜に加熱した半田ごてを近づけるこデル化した条件[7]を用い、非定常伝熱条件を課しとで、液膜内にマランゴニ対流を発生させる。こた。この計算についてはオープンソースであるのような簡易な実験であったものの複数の特異なOpenFOAM[13]を活用し、これらの条件をすべて

流れや振動流れが発生している。この特異な流れ導入して計算した。離散化手法として時間進行にに対しては、数多くの研究が行われてきている

２次精度 backward 法、energy 式の非線形項には[6-10]。一方で液膜内振動流に対しての研究

NVD スキームである gamma スキーム[14]を用[11,12]は少なく、Pettit 宇宙飛行士の実験によっい、その他空間スキームには２次精度線形補間法て確認されたような振動流の発生条件や発生メカを用いた、速度圧力連成には PISO 法[15]を用い

ニズムについては全く判っていない。

た。計算格子の解像度については既報[16]で報告

そこで、今回の研究では Pettit 宇宙飛行士の実の通りである。

験で見られた振動流の発生条件、発生メカニズムを解明することを最終的な目的とし、液膜内マラ3 結果と考察

ンゴニ対流に対する液膜形状と液膜に加える熱量の影響を調査した。

自由表面に沿った温度分布の時間変化を Fig. 2

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

に示す。加熱部から液膜中心に流れが発生するこカニズムの解明のための数値シミュレーションをとで、中央へ向かって高温部分が広がっているこ行った。液膜が大きく凹んでいる場合には左右にとが容易に分かる。さらに、液面形状の変形が大大きく振動する振動流が発生し、液膜形状に関わきい(a)場合、大きく左右に振動する流れが発生し、らずに液膜端に沿って Hydrothermal Wave が発

変形が小さい(b)場合には振動流は発生していな

生することが判った。Pettit 博士が行った宇宙空

いように見える。このことから、液面形状が大き間での液膜の実験では、液膜が大きく凹んでいたく変形している場合に、宇宙実験で見られたようために、大きく左右に振動する振動流が発生したな振動流が発生すると考えられる。

ものだと推察できる。

さらに、加熱部付近の温度分布を詳しくみると、
特に高温となっている加熱壁付近の等温線が波状

引用文献に曲がっている。これは、Hydrothermal wave
によってこのように、温度分布が空間的に振動し[1] A. Hirata, FDMP, 3, 203-230 (2007).

たものであると考えられる。このように、液膜内[2] M. Lappa, Fluids, Materials and Microgravity:

では局所的に時空間的に振動する Hydrothermal

Numerical Techniques and Insights into Physics, 1stwave と左右に大きく振動する流れの２種類の振edition (Elsevier Science, Oxford, U. K., 2005).

動流が混在していることが判った。

[3] H. Kawamura et al., Int. J. Heat Transfer, 134,

加える熱量を変化させた場合の計算も同時に行

031005 (2012).

った。この場合の計算結果についても、加熱壁付[4] D. Pettit, Saturday Morning Science Videos, 2003

近で Hydrothermal Wave が発生し、液面形状がsee

大きく変形している場合にのみ左右に大きく振動

http://mix.msfc.nasa.gov/IMAGES/QTVR/0601211.m

する流れが見られた。本解析結果では、加える熱ov.

量が振動流に対して大きな影響を与えるというこ[5] D. Pettit, Thin Film Physics, 2011 seeとはなかった。このことから、左右に大きく振動

http://physicscentral.com/explore/sots/episode3.cfm.

する流れは液面形状が支配的な因子であることが[6] I. Ueno et al., Proc. of the second Conference on

わかる。

Microfluidics, Toulouse, 10-203 (2010).

液面形状を変化させた場合には、流動の強さが[7] T. Yamamoto et al., Phys. Fluids, 25, 082108

大きく変化する[7]。これは、液面形状が大きく変(2013).

形した場合に、流路が細くなり、質量保存則を満[8] T. Yamamoto et al., Trans. JSASS Aerospace Tech.

たすために流体の速度が上昇する。流体の速度がJapan, 12, Pe_19-Pe_27 (2014).

上昇した場合にはマランゴニ対流によるプリュー[9] H. C. Kuhlmann, Microgravity Sci. Tech., 26,

ム先端において温度勾配がより急峻となり、さら

397-400 (2014).

にマランゴニ駆動力が大きくなり、速度が上昇す[10] L. Fei et al., Microgravity Sci. Tech., (in Press)

る。このため、液面形状を変化させることで液膜

DOI: 10.1007/s12217-016-9523-8.

内に発生する液膜中心に向かう流れが速くなり、[11] I. Ueno and T. Watanabe, I. J. Trans. Phenomena,

液膜中心付近での自由境界層が薄くなったために、

12, 301-306 (2011).

このような左右に大きく振動する流れが発生した[12] D. Limsukhawat et al., Prog. Comput. Fluid Dyn.,のではないかと推測できる。

13, 132-140 (2013).

さらに詳細な結果やメカニズムに対する議論に[13] OpenFOAM, (see http://www.openfoam.com/)

ついては、すでに誌面上[16]で報告済みであり、[14] H. Jasak, Int. J. Numer. Meth. Fluids, 31, 431-449

そちらをご覧頂きたい。

(1999).

4 結言[15] R. I. Issa, J. Comput. Phys., 62, 40-65 (1985).
[16] T. Yamamoto et al., Phys. Fluids, 28, 032106(2016).

液膜内振動マランゴニ対流の発生条件、発生メ

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

Fig. 1 Schematics of a circular liquid film used in thiscalculation.

Fig. 2 Time development of temperature contours
along a free surface: volume ratio is (a) 0.55 and (b)0.95, RED color indicates high temperature and BLUE
color indicates low temperature.

謝辞
本共同研究制度(若手研究者奨励枠)を活用させて頂いたことを、この場を借りて厚く御礼申し上げます。

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 5%



#

$

¹í¿g
£ř£ª£ªśěºªęġĝíęġ%āq~Ë$ĽĆ 1¢¨ņîæ$ĻōP\DB#đ143ė{şĎŅ6~Ë3 îæ$Ñņ6 3P\DB$ŖĔ6Ėßė{ş6ohā$ 3 úvåtľñ6ĸīò 1hā$Èţ#0ĵ»Īİöō$Ģŝč% 88% 1 95%( ċ#hāf$£B@_Z$Ą¸ 1´B@_Z$Ąħ($ŒĞ#oĵ»Īİöō {ĺ3

Ĝ4

1 "

2 !

ũą¬w#Ă23¢ėāq%ÐŀĦ³Ņ ĩ%ĪÔ AB

1"3šį

x"!#0ĉªĕ#c«­"ČÓ "3Ñ$iÔ{Ăāq2|ìČÓ ¬w

Ô{$ė( {ş3$ė{şĎŅ%BO#Ă4ČÓ(Ô{üčL_FZ{ľ &4ċ#ė{ş$Ş#KG

iÔ{$ēŢ6"1 #ĹĎ3-I]_>Č$ĵ»Īİöō ÇÔ43

eÔ{$±À ρAρB

Ę14 2`è$Ô{ť6Ŕĥ$yď

E6­ı3

50%-50%)6Ò­

1mb$:_F_MXU_

2Ŝ  ¢¨ņîæ$ĻōP\DB#đ143¢¨ņîæ%£ "ĹŢĠ/ũùôÅÀ"!$ùĳ6îæ#ĔĎ¢¨ņ<XB/8:^káĴĿ¦Ûqĥ#đ14ψ=

ρ A − ρB
ρ03ė{ş6}đ¢¨ņîærĻőğ# $Ñņ6~Ë 4&ũÎÑ/čĕρ0 %ĂČÓ$¾±À3$"ĊņĞ:_F_MXU_E6đė{şà 6ĹĎìÊ 3®Ş$îæ6rĻ3

P\DB$ā ČÓ%ĂÑĳ63-#hā3¢ hā ė{ş#d3È

ţ6ń)3

Íļ#"3++#ė{ş#03ĵ»Īİöō$ÇÔ hā$©#0Ù~43®ŧĕ•åtĕ#Ę14

3êŗĔ·èğÃ%mb$ Cahn-Hilliard Ã

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

3 [1, 2]hāf$ċ­$Ąöō ė{ş#!
$0#°d3 %5 "í"3t %êxj %Áøui %ŌÀęġ%hāĄ$ĵ»ĪİöōÇÔőğ(Γ %VNYH7_MXU_Eµ %ªSH^A*ò#ń)3-aúvåtŀĦ6đ

WZ3ªSH^AWZ%iŕjÖ "õŁ6ĸ

3Landau-Ginzburg $ĵĒ9KZ=_ 1
ý-

6 − Academic Center for Computing and Media Studies, Kyoto Universityāq$äœèğÃ%mb$ŎĬÃ
Navier-Stokes Ã

3 
hā$ò ¢¨ņöō#*Èţ6ń)3
-Ĉúvêŗ t = 3 # 3ė{şť$ĥ

∂ui
= 0,
∂xi

tŢ (ψ = 0.6)6 Fig. 1 #Ĝā4 "ė{ş#03¢¨ņöō$ÇÔ ĚŃ4
2Ģŝč (ψ∂p∂ui ∂(ui u j )
+
=−
∂xi
∂x j
∂t

O_>tmb3ť$-3)% 88% "`èhāČÓ$ā4 3¢¨ņöōŨó ăµĢ

1 ∂2 ui
∂µ
+
−ψ
+ fi
Re ∂x j∂x j∂xi

ŝč% 95%( $öō %hāf$Ą#03Ú2,/Ä Â&6ŋĂsŏ4-3 Ĳ143Õ$é
É$ęġ[3]# hāfĽ²43£B@3p %Re %[8LZCå_Z 1´B@_Z($9KZ=_;B@_J 

fi %¡Ť3Navier-Stokes Ãf$ŊĤŬ

ė{ş$ĔĐ#0üŉĕ£B@_Z(k<101)Ť%ė{ş°dŤ2ªSH^AWZ$

1ķ ăµ"5ė{ş#09K

œ6Ŧ  āq 36ĹĎ3

Z=_;B@_J Ś¯43{ ė{ş āq6;GPY^?3íę3hāf$Ąöō#%üŉĕB@_Z £ ġ%Ì´B@_Z$ĄŐ 6ľu3Íļ

A_IČ$Ą¸ ë´B@_Z#ėÆ3Ąħ$3-ĢŗÌ{$ŀĦ#%µ"ó§Ćå© Ę14 2[4]ċ­B@_Z$Ąöō

.üŉĕľuÀ$ũØÿ3BQ>

ĵ»ĪİöōÇÔ#*Èţ$ďľ%ĵ»IZÿ6đêŗĠ{#%ĮÇŤ# 2 úĨĪİöō$ĳ ĕ~ËP\DB$Ěģ#£À$ Crank-Nicolson ÿ6šĮÇŤ# 3 úĨÀ°d3 Ï543íęġ%ĵ»Īİ$ Runge-Kutta ÿ6đöō$ #*B@_Z$hāĄ$°díęġ%`#âÜ4hāf¶×ĕ6ń)3-Fig. 2(a)#t = 3 $ x-y çŢ##Ò­430"`÷ĥèhā6³Ņ 3 ψŪ?[_B@_Zūüŉĕ£B@_Z

2$ģèq61283 $ó§Ć$Ą3Ą¸(ŠĶ$ĥtĮ)´B@_Z$Ą{ēïn%è#ìēïn6ł­3Ąħ(ŇĶ$ĥtĮ)${¼6Ĝ+hā%pĀåť#¡63

Fig. 2(b)Fig. 2(c)#%44Ą¸ äœĕ"ť`÷ĥèĕ"hāČÓ6ĭÞ­½ČÓ

Ąħ äœĕ"ť$Ý£6ĜFigure

$H8X_T8>\B@_Z# [8LZ2

Cå% 120":_F_MXU_E$|ì"Çĵ»Īİöō ÇÔ43

ïn%ė{ş$Ŗ¥#Íļ"R]8IL8C6{ 3Ą¸%hāfā4 ' 2£

Cahn-Hilliard Ã#+43MXU_E%

"û, Đ3ť6Ĺ 2$û,

Berti 1#0ĸ54 2 úvŀĦ :_F_#0ĵ»Īİöō Ä Â&4 Ĳ

ĥ "30"t6đ[2]ľñ%¡

143
+
Ąħ äœĕ"ť$2%

ŀĦť%`Ŋ

6ė{ş6ohām¡#|ìŌÀĈ ¡."ė{ş$,$#
.ĸ

1Ą¸ äœĕ"ť%Ą¸#þ0üŉĕ´"ĵ»Īİöō {¼3
{ 34%Ą¸ 1Ąħ( ŒĞ3Ş
#Ą¸#0Ä Â&4ĵ»Īİöō Ąħ$ňŐ #0{ĺ14
$Ş#Ą 9KZ=_6¤ 9KZ=_;B@_J Ś¯43 Ĳ143

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

(a)

4 "

1.4

íęġ%`÷ĥèhāf$iÔ{āq$ė{
şĎŅ6³Ņ  3 úvåtľñ6ĸhāyċ#hāf$Ąöō ĵ»Īİöō$ÇÔ$ŘŎÑ#ń)íęġ%hā$Èţ#0ĵ»Īİöō$Ģŝč ¢¨ņîæ$Ñņ6 143

0.6
0.6

5 

x

1.4

x

1.4

+Ąöō ĵ»Īİöō$Ľ² 1£B@_ZĄ3Ą¸#0ĵ»Īİöō Ä(b)Â&4Ą¸ 1ŒĞĄħ$ňŐ #1.4

0{ĺ143 Ĳ14hāÅÀ$/ç63"!$lćĕ"ãr#0yhāf$Ąöō6 3 ĵ»Īİ
öō6ĳ ĕ# 143

ìÊzð3



0.6
0.6

1. Pine, D. J., Eswar, N., Maer, J. V. andFig. 1 The self-organized structure at t = 3

Goldburg, W. I., “Turbulent(isosurface of ψ = 0.6): (a) without flow,

suppression of spinodal decomposition,”(b)Phys. Rev. A, 29 (1984), pp. 308 - 313.

(a)

with

turbulent

flow.

1.4

2. Berti, S., Boffetta, G., Cencini, M. andVulpiani, A., “Turbulence and coarseningyin active and passive binary mixtures,”Phys. Rev. Lett., 95 (2005), 224501.

0.6
0.6

3. Okamoto, S., Takagi, Y. and Okano, Y.,

x

1.4

“Interaction between phase separationand turbulence in binary mixture,” 2015

(b)

(c)

Asian Symposium on Computational
Heat

Transfer

and

Fluid

Flow(ASCHT2015), ASCHT15-Wed06-005, November 22-25 (2015).
Fig. 2 The distribution of vortex sheets (blue4. Horiuti, K. and Takagi, Y.,line), vortex tubes (red line) and ψ (gray“Identification method for vortex sheet

scale) at t = 3: (b) vortex sheets are

s t r u c t u r e s i n t u r b u l e n t f l o w s, ”

dominant,(c) vortex tubes are dominant.

Phys. Fluids, 17 (2005), 121703.

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

有機系太陽電池への応用を指向した光機能性有機分子の構造と電子構造の相関解明
梅山有和
京都大学大学院工学研究科分子工学専攻

1 緒言と色素増感太陽電池への応用

色素増感太陽電池における増感色素の多くはカルボン酸部位で酸化チタン表面に吸着しているが、その吸着力が十分でなく、酸化チタンから脱着しやすいことが素子の耐久性を下げる一因となっている 1-3。本研究では、ポルフィリンの meso 位の 2

YD2-o-C8T および YD2-o-C8 の THF 溶液中の吸収スペクトルを Figure 2 に示す。YD2-o-C8T の吸収スペクトルは Soret 帯、Q 帯ともに YD2-o-C8 よりもブロード化、長波長シフトしており、光捕集能が改善されていることがわかった。

カ所に 2,6-ビス（オクチルオキシ）フェニル基を、残りの 2 カ所にそれぞれジアリールアミノ基とトロポロニルエチニル基を導入した新規ポルフィリ
ン YD2-o-C8T を設計した（Figure 1）。トロポロンはチタン原子に対してカルボン酸よりも強く配位することが知られており、酸化チタン表面への吸着基として用いることで耐久性の向上が期待される。そして、吸着基がカルボン酸のポルフィリン
YD2-o-C8 を参照化合物として、電池性能及び耐久性の比較を行った 4。

Figure 2. UV-vis absorption spectra of YD2-o-C8T(solid line) and YD2-o-C8 (dashed line) in THF.

浸漬時間、共吸着剤 CDCA の添加量、セル組み
立て後の暗所経時効果により最適化を行ったところ、YD2-o-C8T において JSC = 15.56 mA cm–2, VOC =YD2-o-C8 においてJSC =0.73 V, ff = 0.68, η = 7.7 %、17.32 mA cm–2, VOC = 0.78 V, ff = 0.65, η = 8.8 %という結果が得られた。光捕集能の改善により、
YD2-o-C8T では 700 nm 以降の波長においても光電流発生が確認されたが、より短波長の領域で低Figure 1. Molecular structures of porphyrins.

い IPCE 値（外部量子収率）となっているために
YD2-o-C8 よりも小さい JSC 値になったといえる。

2 結果と考察

しかしながら、連続白色光照射による素子の耐久

2.1 新規ポルフィリン色素の光吸収特性性を調べたところ、500 時間後に YD2-o-C8 では33％低下したのに対し、YD2-o-C8T では 13％しか京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 9

低下しなかった。トロポロンを吸着基として用いよりも狭い HOMO-LUMO ギャップになったと考ることで色素の脱着が大幅に抑制でき、耐久性がえられる。また、吸着部位の LUMO の電子密度向上したと考えられる。

は色素の励起状態と酸化チタンの 3d 軌道の電子
カップリングに影響を与える。YD2-o-C8T では、

2.2 理論計算
ポルフィリンの最安定化構造およびそのフロン
ティア軌道の電子構造について知見を得るために、
密度汎関数法（DFT）による理論計算（B3LYP）を行った。基底関数系に 6-31G(d)を選択し、Gaussian 09 プログラムを用いた。Figure 3 に
YD2-o-C8T および YD2-o-C8 の最高被占軌道（HOMO）と最低空軌道（LUMO）における電子密度分布を示す。YD2-o-C8T および YD2-o-C8 のHOMO-LUMO ギャップを比較すると YD2-o-C8T（2.23 eV）< YD2-o-C8（2.32 eV）となった。HOMOのエネルギーレベルはどちらもほぼ同じであったが、YD2-o-C8T ではトロポロンの強い電子求引性によって LUMO が低くなったために、YD2-o-C8

吸着部位の LUMO の電子密度が高いことから酸化チタンへの良好な吸着と電子注入が期待できる。
一方、YD2-o-C8T の HOMO の電子密度は吸着部位から遠いジアリールアミノ基に主に広がっており、
電荷再結合の抑制も期待される。
このように、
本計算結果は実験結果の理論的解釈の一助となった点で意義がある。

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

π-Aromatic Molecules as Potential Sensitizers in

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

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

Porphyrins

andPhthalocyanines for Photoinduced Charge Separation

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

[3] T. Higashino, H. Imahori, “Porphyrins as Excellent
Dyes

for

Dye-Sensitized

Solar

Cells:

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

[4] T. Higashino, Y. Fujimori, K. Sugiura, Y. Tsuji, S.

Ito, and H. Imahori, “Tropolone as High-Performance

Robust Anchoring Group for Dye-Sensitized Solar
Cell”, Angew. Chem. Int. Ed., 54, 9052-9056 (2015).

Figure 3. Selected molecular orbital diagrams for (a)YD2-o-C8T and (b) YD2-o-C8 obtained by DFT
calculations with B3LYP/6-31G(d). To simplify thecalculations, alkyl chains on the diarylamino groups
were replaced with methyl ones, whereas octyloxygroups on the phenyl groups were replaced withmethoxy ones.

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

バイオフィルム成長に関する三次元モデルの構築
西谷 公秀
大阪大学大学院基礎工学研究科

1 緒言置によって成長段階や周囲流動場がどのように調査するために、従来の二次元数値モデルを三次元

バイオフィルムとは、壁面上に付着した微生物モデルに拡張した数値解析を行った。

が流路内の流体運動による栄養基質物質の供給によって成長する膜状の集合体の総称であり、その2 数値解析手法

成長過程を明らかにして制御することが工学的に重要な課題となっている。例えば医療現場におい解析対象は Fig. 1 に示す三次元矩形チャネルでて使用されるカテーテル内では毒性のあるバイオ

ある。支配方程式は基質を含む溶液には対してはフィルムの増殖を妨げるために微生物を死滅させ

以下に示す連続の式及び Navier-Stokes 方程式をる薬剤を添加することが必要であるが、そのよう用いた。

な添加剤の使用量を最小限にとどめるために移動
現象論的アプローチによるバイオフィルムの成長制御が有効となってくる。Sanchez ら[1]は、マイクロデバイスを用いたバイオフィルムの成長実験∇·u = 0ρ

(1)∂(ϕ su)+ ρ∇(ϕ suu) = −ϕ s ∇p + η∇ 2 (ϕ su) + F∂t(2)を行い、実際の成長挙動の観察に成功している。

ここで、ρは溶液の密度、u は溶液の流速、t は時

しかしながら、成長挙動に関する定量的な評価は間、p は圧力、ηは溶液の粘度、F は溶液とバイオ

実験データのばらつきから困難となっている。一

フィルム間の摩擦力、ψs は溶液の体積分率である。

方、高木ら[2]は、供給物質輸送と流体の運動量輸

一方、バイオフィルムの運動を解くためには、成送を同時に解く二次元数値シミュレーションを実長と代謝を考慮した以下の Cahn-Hilliard 方程式

施し、実験と同様なバイオフィルムの成長を数値
ムの初期条件の依存性については二次元モデルで

を用いた。
⎡
⎞
⎪⎫
∂ϕ n
δ f ⎤ ⎪⎧⎛ µc= ∇ ⎢λϕ n− mϕ n ⎬⎥ +Y ⎨⎜ + m ⎟ϕ n∂t
δϕ n ⎦ ⎩⎪⎝ Y⎠ Kc + c⎣⎭⎪あるために限定的な解析のみがなされていた。そ(3)的に表現しているが、成長核となるバイオフィル

こで本研究では、実際のバイオフィルムの初期配

ここで、ψn はバイオフィルムの体積分率、λはモ
ビリティパラメータ、f は自由エネルギー、Y は
収率、µは最大比増殖速度、m は保持定数、Kc は半飽和定数である。これらの支配方程式は有限体
積法を用いた空間離散化を行い、ルンゲクッタ法
を用いて時間積分を行った。計算格子は直交不等間隔格子を用いて作成し、総格子点数は約780,000 点である。解析条件としては、流入平均速度とチャネル幅で定義されるバルクレイノルズ

Fig. 1 Analysis model of biofilm growth in arectangular channel.

数を 0.1 とし、溶液のシュミット数を 435 とした。
初期状態ではバイオフィルムを底面に付着させ、基質を含んだ溶液を流入させてバイオフィルムを

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

(a)

(a)

(b)

Fig. 2 Initial biofilm distributions on the(b)bottom wall: (a) RUN1, (b) RUN2.
成長させた。初期条件としてψn = 0.5 のバイオフ
ィルムを 16 個、底面にランダムに配置し、Fig. 2に示す 2 種類の配置(RUN1、RUN2)について数値解析を行った。

3 結果と考察バイオフィルムは時間経過に伴い成長し、指数関数的に増加した。二ケースそれぞれのバイオ

t = 140 で 2 %弱である。
フィルムの総量の差は、
これは高木ら[2]が行った二次元数値解析における差より小さい。今回の解析では Sc 数が大きい

Fig. 3 The patterns of stream lines aroundbiofilm

and the grownbiofilm(gray

contour): (a) RUN1, (b) RUN2.

4 結言

ため拡散の影響が小さく、栄養基質の輸送に関しては対流が支配的である。そのため、三次元数

より実際のバイオフィルム成長過程に近い数値モ

値

解析の差が二次元数値解析の差より小さくデルを構築するために、二次元モデルを三次元モなった理由を速度場の観点から議論する。バイデルに拡張し、バイオフィルムの配置依存性を検

オフィルム周囲の流線を Fig. 3 に示すと、流体が討した。その結果、数値解析を三次元に拡張した

バイオフィルムの上方だけでなく横を回りこむことで、溶液がバイオフィルムの横を回りこむ

ように通過していることがわかる。一方、二次挙動を表現できた。よって三次元におけるバイオ

元数値解析ではバイオフィルムの立体的な障害

フィルム成長過程では下流に基質が十分に供給により上方を通過する流れのみが発生し、下流され、バイオフィルム初期配置が成長に与えるの底面付近では基質が不足するが、三次元数値

影響は二次元モデルに比べ小さくなることが示

解析では下流のバイオフィルムに基質が十分にされた。

供給されることになる。これを詳しく観察するために底面近くの基質濃度分布を Fig. 4 に示す。

参考文献

基質濃度はバイオフィルムの上方と横方向で同程度である。よって基質濃度が上方と同様にバ[1] Z. Sanchez, A. Tani, N. Suzuki, R. Kariyama,

イオフィルムの横を通って供給される。ゆえに

H. Kumon and K. Kimbara, “Assessment of

どの初期配置でも基質が二次元数値解析より行change in biofilm architecture by nutrient

き渡りやすく、成長の差が小さくなる結果とな

concentration

った。

microdevice flow system”, J. Biosci. Bioeng.,

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

using

a

multichannel

(a)115, pp.326-331 (2013).

(b)

Fig. 4 The concentration distributions of
substrate around biofim on the cross section of
the plane of x = 0.1 mm: (a) RUN1, (b) RUN2.
[2] 高木洋平、大畑貴嗣、岡野泰則、Z. Sanchez、
金原和秀、”流路内におけるバイオフィルムの成長シミュレーション”、第 63 回理論応用力
学講演会、OS05-01-03 (2014).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 13ԝ‫و‬ड़⏼⏠ ␆␛ ⏫ ⏵ ⎯೑Α଑࠵ՃՕஞѽॾॵ෱ƄխງƄ (λֹ੒Ӂ੒Ӂώ࢙Ӂ٥֋ѕ)
yoshi.saiki@r.hit-u.ac.jp೑Α଑ߘ‫⎯ٴ‬Ϊ⎫ ⏌ ࠵ՃՕஞ⎍֩ਣ⎚⏌ ۚਏ⎰‫ٴߘۿ‬

1.2
1ຢӁ‫⎰⎨ذ‬ಒ࣢⎬⎃⏋ ⎶⏍⎧ ⎅⏌ ⎩ ۪⎉ ⏊ ⏍⏌ ゎ ඕ٥֋⎨0.8

೑Α଑࠵ՃՕஞ⎯‫⎬ݘ߹ ⎧ ⎤ ⏉⎬ޔظ‬ҷ௨⎘ ⎠ࣘ⎨ろ ۚਏ

0-0.2Ιч⎨଑զ⎖ ⏍⏌ ␄ ⍬␥ ⏵ ࣘ 3 ߘ‫ࠃٴ‬ਉ (mod 2π)

=

yn+1
zn+1

=
=

0.4
0.2

⎯ఙ‫ ⎬ڤ‬ਣ⎚⏌ ࡵ࠵ՃՕஞ⎬ મෳ⎘ ⎠ ѽॾ⏓ ⎋⎔ ⎫ ⎤ ⎠ ゎ

xn+1

0.6

λ2, λ3λ1, λ2, λ3⎰⎞⎯ࠦ⎯ۚਏ⏓ ⏃ ⎥ଊߘ‫ٴ‬ຢӁ‫ۚ ⏓ذ‬ત⎘ ろ ⎞⎯ߏ߹⏓

0

200

400

600

800

1000

0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08
-0.1
-0.12

0

200

400

600

800

1000

3xn + r(sin xn + sin yn + sin zn ),य 2: ү࠵ՃՕஞ (࠵Ճ 2-7) ⎯੓ 1,2,3 ␦ ⏞ ␓ ␊ ␑ ‫ॅޱ‬yn − r sin yn + r(1 − cos xn ),(ܳ). ੓ 2,3 ␦ ⏞ ␓␊ ␑ ‫( ॅޱ‬ϓ). ੓ 1 ␦ ⏞ ␓␊ ␑ ‫⎰ॅޱ‬zn − r sin zn + r(1 − sin xn ) + r(1 − sin yn ). ࣢⎬।⎨⎃⏌ ⎍ , ੓ 2, ੓ 3 ⎰।ろ ೦⎯⏃ ⎯⎍֩ਣ⎘ ⎧ ⎅⏌ .
(1)

⏓ ࠟ⏋ Ά⎇ (r = 0.1). ⎔ ⎯ຢӁ‫ؠ଺⎰ذ‬ଯ⎬ ⎰⏧ ⏦⏵ ଯ
⎫ ⎶⏌ ⎿ ⎅⏓ ߟ⎚ (य 1)ゎ ߹‫ ݘ‬3 ⎥⎯␦ ⏞ ␓␊ ␑ ‫⎰ॅޱ‬ろ2π

λ1 = 1.09833, λ2 = −0.00634, λ3 = −0.01383 ⎨⎃ ⏌ ゎ ⎠

z

2π

π

z

0
2π

2π

π
0

y π

0 0

π x

2π

2π

y π

0 3.0

π

x

3.3

2π
π

y

z

π

0
y π

0 0

π x

य 3: ຢӁ‫⎯ذ‬೑ബ෣ (ܳ) ⎩ ⎞ ⎯ұ੒य (ϓ).

2π
0

2π

2π

π

π

0

z

z

2π

0

π
x

2π

0

0

2π

π
x

xn
→ π Ƅ (n → ∞)3
⎫ ⏌ xn ⎯࠴ਖগ⏓ ֆ⏂⏍⎱⏉ ⎅ . ⎔ ⎔ ⎨ろ ߹‫⎨ޔظ⎯ݘ‬⎘ ⎠ (य 3)ゎ ү (y, z) ⎬਻⎘ ⎧ , limn→∞⎰ ∞ ⎯੐⏐⏋ ⎬ 50 ⏓ ๆ⎅⎠ . ⎔ ⎯೑ബ෣ࣘ⎨⎯Օஞ⎯ࠄϳय⎰य 4 ⎯⏉ ⎇ ⎬⎫⏋ ろ ⎞⎯␦ ⏞␓␊ ␑ ‫⎰ॅޱ‬।⎍ 1 ‫ڃ‬⎨ 0 ⎍ 2 ‫( ⏌⎅ ⎧ ⎤⎫ ⎩ڃ‬य 5)ゎ ⎘ ⎠⎍⎤ ⎧ ೑ബ෣ࣘ⎯Օ
0

2π

π
y

ஞ⎍ 2 ߘ‫࠵ࡵ⎯ٴ‬ՃՕஞ⎨⎃ ⏌ ⎔ ⎩ ⎍ߟܱ⎖ ⏍⏌ ゎ ߹‫ݘ‬ろय 6 ⎬ ⏉ ⎤ ⎧ ‫ ⎬ ⎐פ ⎬⎅ڞ‬ਣ‫ ⏌⎚ݚ‬೑ബ෣ࣘ⎯ࡾՃ୅⎰ろ300, 600, 900, 1200 Ѿ౞೹⏓ ‫ ⎬ ⎐פ ⎬⎅ڞ ⏃ ⎧ر‬ਣ‫⎘ݚ‬य 1: ⏧ ⏦⏵ Օஞ⎯ࠄϳय (xyz,xy,xz,yz).

ਞ⎒ ⏌ ⎔ ⎩ ⎍ҷ௨⎖ ⏍ろ ೑ബ෣ࣘ⎬ ⎰೑Α଑ढ़ഇ⎰⎫ ⎅⎡ろ ૳࣢⎯ଊߘ‫ٴ‬ຢӁ‫ ⎩ذ‬Ϊ⎫⏌ ୅⎰೑Α଑ߘ‫⎯ٴ‬Ϊ⎫⏌⎔ ⎩ ⎍⏐⎌⏌ ゎۚਏ⎍֩ਣ⎘ ⎧ ⎅⏌ ୅⎨⎃⏌ ゎ ⎔ ⎯ຢӁ‫⎬ذ‬ඝ⏂‫⎠⏍ ⎿ܛ‬ࡵ࠵ՃՕஞ⎯Ѿୃॅ⏓ ૛⎺⏌ ⎠ ⏂⎬ ࡚⏀೒⎎ ␌⍬⏯ ␑࠵ՃՕஞ⏓ ٠ࡥ⎘ ⎧ ⎞⏍⏊ ⎯␦ ⏞␓␊ ␑ ‫ ⏓ॅޱ‬૛⎺⎠⎩ ⎔ച‫ ה‬WBN [1]⏎ ろ 1,2,3 ⎯೑Α଑ߘ‫࠵⎥ ⏃ ⏓ٴ‬ՃՕஞ⎍֩ਣ⎘ ⎧⎅⏌ ⎔ ⎩⎍ҷ௨⎖ ⏍⎠ (य 2)ゎ Ϊ⎫⏌ ೑Α଑ߘ‫ۚ⎥ ⏃ ⏓ٴ‬ਏ⎯֩ਣ⎯ఙ‫ ⏌⎚⎬ ⎩ ⎔ ⏌੾ ⏓␜⏶␇⏧ ␝ ⏌⎃⎬ڤ‬ゎ

WBN
p (f )(x) :=mod ⏓ ෌߄n=0

⎘ ⎧ ۪⎉ ⏌ ⎩ ⎔ ⎯ࠃਉ⎨ x ൘‫ ⎰ ⎧⎅⎥⎬ڽ‬1 Ѿ⎯౞೹⎨Ҟ
⎮ 3 ఝ⎖ ⏍⏌ ⎔ ⎩ ⎬મෳ⎘ ろ x ߫⎬Юઐଯ⎫ ೑ബ෣⏓ ஔ଑

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

N−1ŵn,N f (xn )(2)⏓ ๆ⎅⏌ ゎ ⎔ ⎔ ⎨2πŵn,N =

2π
z

π

y

π
0

2π

0

π

x π

0
3.0

2π

3.1 π

3.2

x

3.3

wp (t) :=

2π


�
�
exp p −1 p ,t (1−t)

0,
z

π

z

(3)

⎨⎃ ⏋ ろy

2π 0

wp (n/N )N −1
Σj=0 w(j/N )π

for t ∈ (0, 1)(4)

for t ∈/ (0, 1).

⎨⎃ ⏌ (Ιч⎯‫ ⎰⎨ޔظ‬p = 1 ⎩ ⎘ ⎠ ). ೑ബ෣ࣘ⎯ࡵ࠵ՃՕஞ⎯ y ढ़ഇ⎌⏊ ֆ⏂⎠ Ѿୃॅ ρy ⎰ろ

0
3.0

य 4:

3.1 π

3.2

x

0

3.3

0

2π

π
y

ρy = 0.1731838339716555313719

x ߫ ⎬ Ю ઐ ଯ ⎫ ೑ ബ ෣ ࣘ ⎨ ⎯ Օ ஞ ⎯ ࠄ ϳ य z ढ़ഇ⎌⏊ ֆ⏂⎠ Ѿୃॅ ρz ⎰ろ(xyz,xy,xz,yz). ೑ബ෣ࣘ⎯Օஞ⎰೑ബ෣⎌⏊ Юઐ൘‫⎬ڽ‬೑Α଑⎨⎃⏌ ⎠⏂ろ ೑ബ෣ࣘ⎯୅⎌⏊ ౞೹⎖ ⎜⏌ ඡ⎬ॅખ‫ ⎬ܴڬ‬৷୻⎚⏌ ഇ⎡⎒ ೑ബ෣⎌⏊ ๸⏍⎧ ⎘ ⎿ ⎇ ゎ ⎞ ⎔ ⎨ろࠃਉඡ⎬ ઩λ೑ബ෣⎬ ω⎎ ෸⎘ ⎧ Օஞ⏓ ֆ⏂⎠ ゎ
1.065
1.06491.06481.06471.0646
1.06451.06441.06431.0642
1.06411.064

ρz = 0.1600374953195806712604(6)⎩ ⎫⎤ ⎠ゎ ⎔ ⎯ࡵ࠵ՃՕஞ⎯ਣ‫⎍ݚ‬Ϊ⎫⏌ ೑Α଑ߘ‫⏃ ⏓ٴ‬⎥࠵ՃՕஞ⎯ᘥවड़⎫ ⏊ ⎴⎬ ԝ‫و‬ड़⏼⏠ ␆␛ ⏫ ⏵ ⎯Ք‫ٺ‬⎬ ⎫ ⎤ ⎧ ⎅⏌ ⎩ า২⎖ ⏍⏌ ゎ

0.01
0.005

‫۪ދ‬ഓ٤

0

λ2

λ1

(5)

-0.005
0

25000

50000
t

75000 100000

-0.01

0

25000

50000
t

75000 100000

0.01[1] S. Das, C. B. Dock, Y. Saiki, M. Salgado-Flores,
E. Sander, J. Wu and J. A. Yorke. Measuring
quasiperiodicity, Europhysics Letters 114, 40005,
2016.

0.005
λ3

0
-0.005
-0.01

य 5:

0

25000

50000
t

75000 100000

೑ ബ ෣ ࣘ ⎨ ⎯ Օ ஞ ⎯ ␦ ⏞ ␓ ␊ ␑ ‫( ॅ ޱ‬λ1=

1.0646, λ2 = 0.0000, λ3 = 0.0000).

2π
z

π
0
2π

y π

0 3.0

π

x

3.3

य 6: ೑ബ෣ࣘ⎬ ‫ ⎠ ⎤ ⎩ ⎬ ⎐פ ⎬⎅ڞ‬13 ⎯ࡾՃ୅ࡍ‫⏓܄‬300, 600, 900, 1200 Ѿ౞೹⎖ ⎜⎠ ࡍ‫܄‬. ␦ ⏞ ␓␊ ␑ ‫ॅޱ‬(य 5) ⎩ ‫ ⎩ ⏌⎜⏐܄‬೑ബ෣ࣘ⎯⎶⏌ ⎿ ⎅⎰ࡵ࠵Ճଯ⎨⎃⏌⎔ ⎩ ⎍ߟܱ⎖ ⏍⏌ ゎ

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

運動量輸送と熱輸送の非相似性の最大化ôó Õ
­ť­°­°Ŧ¦ġ¾°ĠĦģýľ×¶ç

1 緒言|qĚŋīÙű!Ż'jć Ğţ§ŻŘĵ Æ3ċŐ/ŖÅ§3òâ²

8N\=a ñĆĕ ŉČ-ŻčkãfŻăŤé ş3ĕcíìŖÅ§Óľc ü¾°ýĂ)-0

3ïë/ ú­Čœl Ż Ì/żčkã Óľ3c/)!Żr

 GMRES ą3.z0MYaK_ą3ĕćq ćÜÛ3pĈŻ/!©'õżŖÅ+ĊÅĪ «é!ć0ì'Żnč3tř/Ñņ /żŻ

#EP_ìRa[8įéŻª¤ĜìH³Ē!bŀūōų/0/żs

7QC7R¬ŰÆ ķÉĳŅ0/ªŮ¨Ė"Żjć!»ćā$ŴłŷčŒŕľ 3

öm3ċŐ/¦ÄºŢŻ;Za<_ą3ĕñ/ ŻbìŻśŠ Œŕ tř*'ÇìĥÆ3Ūèżì 86[5Eŏ ŻĳùŻåæä® ©­3ÞŌ¿!ŧżìĥÆĭ Ā¿!ŋīý iě'ż

D]_V!śŠŒŕčŒŕţŌ¿ O(10-15)ĥÅ/żévŋīAaL!

ĝoÓ -ĲŶĚ, ğ-0/żFortran90 r×ŻIntel A_P6Z3ĕAôĠĦ ěĚ!Ż

ĝoÓ3½ŻŭĝoÓ3_ P 6 \ 3 Ń   ż ŋ ī A a L ! MPI tř/Ź|qĚ!ŻªŮĊÅ b² ²

OpenMP 3ı(2hżĴŋ}ē+0ÃŃÃøţć0ŻªŮåæä®3ňĪ UEK}ē! IDL # Paraview 3¸ u -ŷªŮčćõ3Î/źŖÅ§ĕż

3Ň / /ż ,ðŝŖÅ§3Ă)/  ~÷0"Żnčtř3ěĚ

3 結果と考察ćqÏ àšÊħòË0/żŽ!]6O\FéŹÃø ŖÅ¿ 

2 計算手法ćőŷ ¦ źÎ-0ðŝďÔ E:ZaèŚ8N\=aèŚ3Ģ

cdb² ĊÅ¿3lſþyÃŃÃø

/ż³ķ!»ćďÔ3Ģ.ŻE:ZaèŚ

ţ>8IKć3µŎŻE:ZaŹĊÅź Ķè

#8N\=aèŚ!Īv3/żbìjŚđŹªŮčćõĝÈ/ź8N\=a ĶćďÔ!Żgļ!±{!bĿ*èŚđ ¿3ěĚăŤé/żáİ µŎŻŔĝoÓ )Ż%&Īv3

/ŖÅ§!²Â/* ŻĮömĢ/  /żµŻôĠĦÎ-Żŭ¡Ĺćq /Řĵ ÆĊÅ0ðŝďÔ!jćďÔ,.*ŴłŷE:ZĤćßèìĥÆ3ĕż«ą3ĕŻ¨ĖaèŚ p8N\=aèŚŜ×0

öm3ĻÖ/ ŻăŤé ĩb« Ŭ/  2/żž!p]6O\Fé //wėömŹŭķÉ 96Za`Z?Z_DY

E:ZaèŚ3Ģ.Żć0ìbüìĥÆź3·~Ż 03évĚŊ,Į3ōžþy ŖÅ§Î-0ðŝ.ěĚăŤé3ð­/ðŝŖÅ§3ÎżďÔ*S^IK/żú)p]6O\F

16 − Academic Center for Computing and Media Studies, Kyoto Universityé!»ćŊ ðŝďÔÎ-0Ż]6O\
Fé ©/ć0ìbüžþyŊ Ē: Scalar dissipation [Optimal]: Energy dissipation [Optimal]

10-1: Scalar dissipation [Turbulent]: Energy dissipation [Turbulent]0/żſ!žþyŊ ĊÅ#ŖÅT>K\Á3Ģ/ż

Ŋ!ŻĪĊÅ¿3ñ/

10-2

ÃŃÃøţ /Żžþyŭ¡Ĺ²ÂŖÅ§Ż8N\=aèŚb² öm dª¤Ĝì čćõ3ð­/ðŝďÔ[1]ŵo.Ż­E@a\ ÐĔ^a\,û×0/żbìôĠĦ!Ż-]6O\Fé3©/ ,ſþyûŗ3ñ/Ŋ ðŝďÔÎ-0żŷ]6O\F

10-3

10-4
10

100

1000

10000

Fig1. Scalar and energy dissipation for optimal and turbulent states.

éŻžþyŊ!ůČ/żƀ!¼
ĜÌ /ć0 ûŗ3Ģ/żĉûŗ0.020

ň!ŖÅşJ_G\ ĩže«Š𝑄𝑄 > 0

0.015

3ĕ.Ż ĪvŮ!EP_ìĉÅ «×𝜔𝜔!! = 𝜔𝜔! −𝜔𝜔! !" 3ĕŁl -0

1.456

0.010

/żÝ~0Ĭďĉûŗ!ć0ì3.Ż0- ­Ş ÿ 𝜔𝜔!! 3Ģ/ 
2/żÍŻ

[3D]

1.450

0.005

247.5
1.446
245 246 247 248 249 250

ĸĉ!Ã£ć ĉÅµ0.000Ĩ EP_ìĉÅ3*,EP_ìx//
Ź5_HB6>^_ĉź
ż
'Ż

0

aèŚ žþyŊ-Ŵł©Ų¥3Ģ/żŸ ĪvŮ!ĉĬ Ĝd (ŉ´
0/ Żµę ĪvŮ!ĉĬĄ)-0/ ŏ/Ļ-0/żôĠĦ,Î-0ðŝďÔ!ú)ŷnčĎÓ3ñ/ ŻnčtřÏ ŉČ-Ļ/Ż ŖÅ§ * 3ćqÏ,
³ũŜ×/ !ū/1ż -ŻăŤé ðŝų3ŊyÁ/  2/ż 0!Ż5_ţŨ ŻĽî¯¢/4êć,,.Đ

100

200

300

400

500

Fig2. Scalar dissipation for 𝑅𝑅𝑅𝑅 = 0 − 500.

ęŸ ĪvŮ!¹Ø E:ZaèŚ#8N\=

HB6>^_ĉÀ l -0/ĊÅ ĪĊķ[2D](×10 -2 )

1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

0.0

0.5

1.0

1.5
z

2.0

2.5

3.0

-1.0

Fig 3. Temperature and velocity field of the optimal state at 𝑅𝑅𝑅𝑅 = 200.

,Î-0ğŇ Ž/5_HB6>^_ĉ,/nčtř X:MFW3Òĕ/ ,ŻënčtřÚń ŢĘĺ /  òË /ż

4 参考文献[1] Hassanzadeh, P., Chini, G. P. and Doering, C. R.
“Wall to wall optimal transport”, J. Fluid Mech., 751,

Fig4. Vortex tubes and regions in which increase in scalar (white) andenergy (black) dissipation is observed.

627–662, 2014.

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

高分子特性を考慮した界面活性剤添加溶液の抵抗低減効果ŕØ é°
ň¢¢ŊĄ¬¢ăĊĆ

1 緒言

2 数値解析手法

êj jŐ$ņ#%Ĕ½ËÌÄÃ ùÚăĊ%)%+#ŏě½YáďÐa

$ÄÃ#. "6HS9XÊ&}

ê$ DNS 4ğÍĿÐć²%f[#ąĹÊ ³ ī21Ķ¯$ê2%ķ± ķ

ę$²

aê÷¾1 $aêËÌÄÃ4iî1Ðç$Y#êj^#®Ğ"ŕy òë4í~1 .ËÌÄÃ î©1E

N>Ü 1EN>Ü%ºł"í~|4~1 #.äİþ¥Ó#ÄÃiîÜ4¹/21/ĂåK5MQ5Vıĵ-õdÉ

Navier-Stokes Ðć²1

∂ui
=0
∂xi

(1)

2
∂ui ∂(ui u j )∂p 1 ∂ ui+=−+∂xi Re ∂x 2j∂t∂x j"!))"yŁ$¼ú ×¸21

1uΩRo j k+ εijkEN>Ü#.1ÄÃiî$Ņ%_#ŕy(2)

ó±#m¡1 ĝ/21 ħĖ"âàui %ķ±t %Ôņp %}Ωi %ģ

%y#ĨÒ "2"ÈÞ21

ķ±0Ï%ôãs21T5ÄÃiîO8G>N%ŕy $ĳq$Ōĉ#Į$´4ąIS>Ï Re % 3000 $* .1ķ±û« *$~ êj6HS9X$

U=LXÏ Ro #Ň%1/Ro = 0  0.5

ŕy µ½6HS9X'$É#.1Ü$cĈĮ´±$m¡½4ßĦĮ% z įÐ$

œ#y/21*#]ĤÛŒ%ďÐþ"Yáaê4§Ī

EN>Ü$ýø#%ŕy $T7U<Xö½

1+YĲ

2π$ČÐjŒ ì Ň31 ĝ/21ÄÃiîĥđÝ % 2563!$."Ô#ýø1 4bñ1+#%

úōÎ Runge-Kutta ç&ŕy $ĚÔņ aê^$ïĻ$Ôņ=;XÍĿÐć²%Öŉ­yç4

Crank-Nicholson ç#.Ôņĉy»ġ1 21

ã#ùÀ2aê#§ŕy 4Ĵ~[1] ¤Ŕþ#2/$ö½ł4Ôċņþ#ĩ(aê^$ŕy $Ļ4Ĵĭ¦1 %Ŏ y"ßĦ%"2""ŕy í~òë$aêĤÛ%ŕy ¼}$ýÚăĊ%Õ,ĕ"8IG8Saê$YªÐć²4 Navier-Stokes Ðć² 8BMRV:

1YáďÐaê4§Ī YįÐ#Į4

1 ÚăĊ%aê /ŕy '$·ő

~    Į a ê $ Ā Ç Ï p ĥ đ (Direct*$4ĝ¿1 one-way 8BMRV: ŕ

Numerical Simulation, DNS)4¤ÑùÀ2

y Ń$PDS#%Finitely Extensibleaê#ŕy 4ÂtĴĭ1 #.0Nonlinear Elastic model (FENE)[2]4úŕy $µ½6HS9X aêïàĸ $āŇ4

$ FENE PDS4úĩ(S9XEl %f[$.#ĥđ21

S4Ňn1

18 − Academic Center for Computing and Media Studies, Kyoto Universityŕy $µ½6H

aê^$ïàĸ4Åx1+#ã²]/nkTb(N ) N −1 ⎛ [r (m) ]2 ⎞El = −∑ ln ⎜1− b(N ) ⎟
2
⎝
⎠
m=1

21ķ±ĿCV@S$Ďc\ł Q 4đx

(3)n %ŕy $Ï¦±k % Boltzman£ÏT %Ę§ð±N %LX>Ï(e$ĤÛ% N = 2)m %LX>$č /$üb(N)&r(m)%22°Ġ÷¾$ńôãs2
JH$Õhńń LX>ņĬō1ŕy $ÕÙčņĬō$c`4 b = 3000

1
S S + Ωij Ω ji ,2 ij ji
1 ⎛ ∂u ∂u ⎞Sij = ⎜⎜ i + j ⎟⎟ ,2 ⎝ ∂x j ∂xi ⎠1 ⎛ ∂u ∂u ⎞
Ωij = ⎜⎜ i − j ⎟⎟2 ⎝ ∂x j ∂xi ⎠

Q=

ŕ

(

)

(4)

Kolmogorov $Ôņ=;$Ďc\ł Q 4úaê^$ïàĸ4

XSτ (= 0.77)£Ĝ21 Weisenberg Ï(We =ĢėÜ4 Fig. 2 #ąÚăĊ$ DNS

λ / τ )410

DXA% Re Ï i+01 Q $p 

y $ÕĚÔņληη

/90$pĥđ4ğ" Fig. 2(a)#ąĮô$$ė

3 結果と考察Ü#èÿ1 ?U#ĳ Q = 1 $ďpŐ%«÷$àĸ Åx2  Q = 1.5 $ďpŐWe Ï

%Ē÷$àĸ Åx2 0ï«#æ$āŇ4ąWe Ï ~1 ŕy $µ½ïĒ ¶À21 [3] 3 1YÐ

6HS9X,äl~1

Į4~ Fig. 2(b)$ėÜ%ïàĸ ujþ

Figure 1 #ŕy $µ½6HS9X

3 1 aê#Į ~31  We Ï# ŕy $µ½6HS9X'$É /#~
1

(a)

3 1Į$Ü% Navier-Stokes²# } ]/2"3 "=;XS$àĸ#r 6HS9Xèt ğ321aê^%6HS9X8=;XF#.èt26HS9X ºĖ"aêï$=;XS)gĽ2$ėÜ=;XS$¨ŕy# ,êj /ŕy '$6HS9XgĽĮ$Ü#.oĺ21 ĝ/21
(b)Fig. 2 The isosurfaces of Q at t = 40 (outerFig. 1 The rotation strength dependency of
Weisenberg number on elastic energy.

surface with Q = 1 and inner ones with Q =1.5 : (a) 1/Ro = 0 (without rotation), (b) 1/Ro
= 0.5.

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

#ē"3z×$ àĸ ºĖ2"

3 0Fig. 1 $ėÜ kĝ1 ŕy %aê^$ àĸ1
ï«$Ĳ#Ō)06HS9X8=;XF$ļć# aêĻ$6HS9X ŕy $µ½
6HS9X#,É21 ĝ/21

4 結言ÚăĊ%ŕy í~òë^# ÄÃiîÜ#" 1ŕy $Æ4ĩ(1+#Į4~YáďÐaê$ĀÇÏpĥđ4ğùÀ2aê#ŕy 4Ĵ~ŋ$ŕyÆ ŕy -êjņ$6HS9XgĽ#ĤÛ4ğĮ4~1 #.aê^$ïàĸ$ºĖ Á{21  =;X
S$ïàĸ1«÷$ï«Z#ŕy Ń Ōĉaê$6HS9X$Yľ ŕy $µ½6H
S9X' É21

Ò/ #"e§Ī ê2%YáďÐaê w Ē v ê 2 -  Į c ŀ w Đ v ê 2(Taylor-Couette ê2)# ,á"ĤÛ4ğŕy $ĿWŌĉ ÄÃiîÜ$Ňn4
ĩ(1b£1

参考文献
[1] P. K. Prasinski, B. J. Boersma, F. T. M.
Nieuwstadt, M. A. Hulsen, B. H. A. A. Van
Den Brule and J. C. R. Hunt, “Turbulent
channel

flow

near

maximumdragreduction: simulations, experiments andmechanisms”, J. Fluid Mech., 490, pp.
251-291 (2003).
[2] R. B. Bird, C. F. Curtiss, R. C. Armstrong
and O. Hassager, Dynamics of Polymeric
Liquids, 2nd Edition, Volume 2, Willey(1989).
[3] K. Horiuti and Y. Takagi, “Identification
method for vortex sheet structures inturbulent flows”, Phys. Fluids, 17, 121703(2005).

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

周囲気体を考慮した HZ 液柱内温度差マランゴニ対流場の数値解析
中村 翔一
東京理科大学大学院理工学研究科機械工学専攻

1 緒言

2 数値解析手法

自由表面上に温度差を与えると，それに起因した表面張力差を駆動力とする対流が発生する．こ

数値解析で用いた計算対象を図１に示す．これの流れは，マランゴニ対流と呼ばれており，特には，液柱(Pr=4)の周りにシールドを取り付けその微小重力下やマイクロスケールにおいて支配的に中にアルゴンガス(Pr=0．67)を充満させている．

なることが知られている．

支配方程式は円筒座標系非圧縮非定常の３次元ナSmith と Davis [1]は，薄液膜において自由表面に

ビエストークス方程式，エネルギー方程式，連続

対して水平方向の温度勾配が発生している場合，の式であり不等分割スタッガード格子を用いて解

マランゴニ対流により hydrothermal wave(以下いた．また，有限差分法を用い 2 次精度中心差分

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

により離散化した．連続の式，ナビエストークス

用いて発見した．現在，図 1 のような Half-Zone

方程式のカップリングにはフラクショナルステッ

液柱(以下 HZ 液柱)はマランゴニ対流の不安定性

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

を研究するための代表的な系となっている．HZ

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

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

進行方法は，液柱には 2 次精度の数としてマランゴニ数(以下 Ma)がよく用いられAdamus-Bashforth 法，周囲気体は，軸方向の粘ている．Kamotani ら[2]は HZ 液柱において高プ

性項に 2 次精度の Cranck-Nicolson 法その他の項

ラントル数流体(Pr=24～49)を用いて実験的に自は Euler 法 を 用 い た ． 粒 子 の 計 算 方 法 は

由表面での熱輸送を示すパラメータであるビオ数

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

(以下 Bi)が臨界値に与える影響を示した．2014

この計算コードは，Open MP と MPI を用いたハ

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

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

Kamotani ら[2]の実験を元に Bi に着目している．
さらに，
2018 年から実施予定の日欧共同宇宙においても，Bi と可視化用トレーサ粒子が一本の線内3 結果と考察に集合する現象である粒子集合現象に焦点を当て

今までの先行研究では、Bi を考慮し液柱内流れた研究が行われる予定である．しかし，実験装置の不安定性を調べる際計算コストが増加するため、の制約上得られない情報も多く，数値解析による周囲気体を無視し簡略化した境界条件を定義し，アプローチが必要である．また，液柱だけでなく周囲気体の数値解析はほとんど行われこなかった．

周囲気体も着目しているため，計算コストが非常

まず，二次元定常流から三次元振動流へ遷移すに高くなりスーパーコンピュータの利用が不可欠る臨界マランゴニ数 Macr と低温壁温度 Tc の関係である．本研究では，実験に近い条件で従来よりを示す．臨界マランゴニ数 Macr の算出方法はも広範囲な領域での数値解析を実施し，周囲気体

Imaishi ら [4]の方法を用いた．局所的なビオ数

が HZ 液柱内流れの不安定性に与える影響を調べを軸方向に平均したビオ数 Biave は Macr 上での値ることを目的としている．

である．今回の数値シミュレーション結果（図 2）

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

Macr と Tc の関係を見ると Tc が高くなるにつれてMacr が小さくなり振動流へ遷移しやすくなって25000

いる．これは地上実験でも見られている傾向であ

20000
Macr

り[2]，表面での熱損失により表面流速が加速され


2
1
0.5
closed dot ; constant viscosity
open dot ; variable viscosity

るためである．次に，アスペクト比 Γ（液柱高さ/

15000半径）の影響を見ると Γ が大きくなるにつれ Tc

10000

依存性が高くなる傾向が見られた．これは，低温
壁温度 Tc との臨界高温壁温度 Thcr の関係（図 3）2030により説明できる．これを見ると，アスペクト比が大きくなるにつれ，Thcr が小さくなっているこ

40

50

Tc[℃]

60

図 2．低温壁温度と臨界値の関係
横軸: 低温壁温度
縦軸: 臨界マランゴニ数（対流場の遷移点）とがわかる．そのため，Γ が小さくなるにつれ液柱近くの気体温度が上昇し，液柱表面での Bi が低下する．つまり，Tc を同じよう変化させていっても Γ により Bi が異なるためこのような傾向となっていると考えられる．さらに、動粘度温度依

60

存性の影響を見た．これは，数値シミュレーショ
Thcr[℃]

ンで調査できないパラメータであり，図 1、2 の中塗りの点が動粘度温度依存性なし、中抜きの点が動粘度温度依存性ありの結果である．これによ

40


2
1
0.5
closed dot ; constant viscosity
open dot ; variable viscosity

20

り動粘度温度依存性を考慮したことにより Macrと Thcr が小さくなることを示した．

20

4 結言

30

40

Tc[℃]

50

図 3．低温壁温度と臨界高温壁温度の関係
横軸: 低温壁温度
縦軸: 臨界高温壁温度（対流場の遷移点）

より実験に近い条件である周囲気体を考慮した
HZ 液柱を用いたい数値シミュレーションを行い，地上実験と同じような傾向を得ることができた．
さらに，Γ による影響を調べ今後の実験結果に繋げることができた．しかし，この周囲気体を考慮

参考文献

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

[1]

M．K．Smith and S．H．Davis， J Fluid
Mech， 132， 119-144 (1983)[2]

Y． Kamotani_et_al， Int． J． Heat
Transfer， 46， 3211-3220 (2003)．

[3]

E．Hofmann， H．Kuhlmann， Phys．
Fluids， 23， 072106 (2011)[4]

N．Imaishi_et_al， Crystal Growth 230，
164-171 (2001)

図 1 計算対象

最後に，本研究制度（若手奨励枠）を活用させて頂きましたことに大変感謝致します．

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

60

高プラントル数を用いた HZ 液柱内温度差マランゴニ対流の数値解析
カオス・乱流化過程の解明に向けて
茂木 孝介
東京理科大学理工学研究科機械工学専攻

1. 緒言
自由界面を有する系において，気液界面に温度差や濃度差をつけることに起因して，表面張力差を駆動力とする対流が発生する．この対流はマラ
ンゴニ対流と呼ばれ，微小重力環境下におけるマ
イクロスケール流体の駆動力の一つとして注目さFig. 1 The schematic image of a liquid bridge.

れている．一方で，floating-zone 法のような高品
質な材料生成過程において，マランゴニ対流の不2. 計算条件

安定性によって生じる振動流が材料生成に悪影響を及ぼすことが知られている 1 ．故に，マランゴ

液柱の高さを H，半径を R とし，アスペクト比Г = H/R と定義する．今回の計算では，Г = 2.0 とニ対流の不安定性に関する理解，制御の点で研究が行なわれてきた．

し，液柱は表面の変形を考慮せず，常にストレー

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

トな形状を保持しているものとする．また，計算系全体は無重力状態であることを想定した．
流体は非圧縮性ニュートン流体であるとし，Pr

物理実験に共同研究グループとして参画している． = 16 とした．無次元化された流体の運動と熱輸送
今回，宇宙実験のテーマの一つである，「マランゴを記述する支配方程式は，Navier-Stokes 方程式，
ニ対流におけるカオス・乱流とその遷移過程」の解明に向けて，直接数値シミュレーションを行な

連続の式，エネルギー方程式となる．

った．

熱の授受を示す無次元数，ビオ数(Bi)で表した．マラ

本研究では，計算対象として floating-zone 法のまた，液柱自由表面上での熱の授受を自由界面でのンゴニ対流の強さを示す無次元数，マランゴニ数

代表的なモデルである half-zone モデル（ HZ モデル，図 1 ）を計算対象とした．上下ロッド間に液(Ma)をとし，ビオ数とマランゴニ数をパラメトリッ

柱を形成させ，上部壁を加熱，下部壁を冷却することで液柱自由表面上に温度勾配を付与し，マラ

3. 結果

ンゴニ対流を発生させるモデル系である．宇宙実
験の結果との比較を行なうために，高いプラント
ル数( Prandtl 数，Pr = ν/κ，ν：流体の動粘度，κ：
流体の熱拡散率)で，かつ，地上実験と比べて長い

クに変化させて計算を行なった．

図 2 に Bi を変化させたときの臨界マランゴニ数(Mac)と各対流場の遷移過程を示す．Mac は各 Bi 数における臨界マランゴニ数を示す．SW は standing

waves である．
TW1(a)，
1(b)は traveling waves である．

液柱を対象とした．また，周囲気体への熱損失を考慮し，対流場が定常な流れから非定常な流れへ

臨界値の傾向は先行研究で実施された、Pr = 16，Г =と遷移する臨界値を調査した．また，遷移後の十分発達した対流場の構造を調査した．

ている 2．Bi が 1 以上の領域において，対流場は SW

1.8 の線形安定性解析の結果と類似した傾向を示しへと安定的に遷移し，維持されていることが確認さ京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 23

れた．また，臨界値からやや離れた Ma においては，
臨界値直上の SW とは傾向の異なる別種の SW が確認された．図 3 は臨界直後に見られた液柱自由表面
上の温度変動分布の 1 周期分の時間変化を観察した(a)

ものである．臨界マランゴニ数に近いほうに着目する( 図 3( a ) )．上部ロッド側から cold spot が現れ，
下部ロッドに向かい，それと入れ替わる形で，今度は hot spot が現れる．hot spot も cold spot 同様に上部から下部へと流れていき，再び cold spot が現れるという周期的な流れとなっている．
次に臨界マランゴニ数から離れたマランゴニ

数の対流場に着目する( 図 3( b ) )．図 3(a)と大きく異なるのは cold spot の形状である．cold spot の出現と同時に一対の斜交した cold line が現れ，
cold

(b)

Fig. 3 The snapshot of temperature deviation
distribution on the free surface,( a ) Ma = 19000, Bi = 2.0 ( b ) Ma=36000, Bi = 2.0.

spot が現れたところで交差した状態を維持したまま下部ロッドに向かっている様子が確認できる．
これは Matsugase et al.によって報告された，‘X’状構造に非常に似通った構造である 3．本計算結果
は実験で確認された‘X’状構造を，初めて数値シミュレーションによって再現することができた．
この温度変動分布の構造を明らかにするために，周方向波数: m ≥ 2 による空間フーリエ分解を施した(図 4)．これを見てみると，オリジナルのデータでもっていた，cold line の斜交している構造

Fig. 4 Space Fourier decomposition for m ≥ 2(Ma = 36000, Bi = 2.0).

4. 結言を維持していることがよくわかる．これより，こPr = 16 の流体で，無重力状態を想定したスト
レートな長液柱を対象とし，周囲気体との熱損失の‘X’状構造は基本波によって構築されたものでを変化させた際の数値シミュレーションを行なっはなく，より高い周方向モード数によって構築されていることがわかった．これは対流場の高調波

た．ビオ数が 1 以上において，対流場が振動流へ成分の影響が増大したことによるものであることを示唆する結果である．

と遷移すると，対流場の形態は SW へと遷移し，維持されることを明らかにした．SW の対流場を強くしていくと，‘X’状構造を持った対流場を数値解析で初めて再現することができた．また，‘X’
状構造は対流場の高調波成分の影響が増大したことによるものであることを示唆した．
今後の予定は，周囲気体から熱を受け取った際の対流場の遷移過程を調査し，今後の宇宙実験への一つの知見として貢献していく予定である．

5. 参考文献
1 ) A. Cröll et al., J. Cryst. Growth, vol. 191, pp. 365–376,
Fig. 2 Critical stability curves as a function Bi
for Pr = 16, Г = 2.0.

1998.

2 ) B. Xun et al., Int. J. Heat Mass Transf., vol. 54, no. 9–
10, pp. 1698–1705, 2011.

3 ) T. Matsugase et al, Int. J. Heat Mass Transf., vol. 89,
pp. 903–912, 2015.

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

パターン学習を用いた経験的セル・オートマトン構成法とそれによる偏微分方程式モデルの模倣
中野直人 ∗∗

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

本稿は川原田茜講師（京都教育大学数学科）と宮路智行特任講師（明治大学研究・知財戦略機構）との共同研究に基づく．

1

接セルの状態から決定するというモデルである．例
えば，空間 1 次元の M 状態 3 近傍 CA では，(a, b, c) → d

はじめに(a, b, c, d = 0, 1, . . . , M − 1)という局所規則という写像を用いる．これは，現時

複雑な挙動を見せる現象に対する予測モデルの構
築は自然科学のみならず様々な分野における強い関心事である．いわゆる非線形性が強い，もしくは自由度の大きな力学系ではその予測は一般に難しく，数理科学の重要な問題として長く研究されてきた．また，計測技術の向上によって大量のデータが取得できるようになり，そのデータから対象とする現象の法則を経験的にモデル化する試みがなされてきている．この手法では，所与のデータセットに対象とする現象が表現するあらゆる振る舞いが含まれているとみなし，機械学習的なアプローチでモデル方程式を構成する．これには大規模なデータ処理とモデル
パラメータ決定を必要とするが，大規模計算機の利用など，昨今の計算機能力の向上が経験的モデリン
グの実行を可能としている．
しかし，複雑な現象であればあるほど予測モデルを立てるのは難しく，その表現のためにそのようなモデルにはしばしば多くのパラメータが含まれる．
それらのパラメータを実際のデータにフィットさせることは困難を伴う．そのため，現象の遷移の「骨子」を表すようなセル・オートマトン（CA）によるモデリングは，ノンパラメトリック解析の一つとして，方法論の確立が必要とされている．経験的に構成された局所規則を調べることによって対象となる現象のメカニズムの解明につなげるのである．

CA モデルとは，空間を有限サイズのセルで区切り，次の時間ステップでの状態を，現時点における隣点での状態が b であり，かつ左隣と右隣の状態がそれぞれ a と c となる場合，次の時点での状態を d で与える．このように，CA は局所規則によって各セルの状態を更新する非常にシンプルな数理モデルではあるが，game of life や von Neumann の自己複製機械などの複雑なパターンを表現するものもあり，豊富なパターン生成を為し得るという意味では現象の模倣に用いるに足るポテンシャルを有している．
所与のデータから経験的に CA をモデルする手法
は，最初に Kawaharada and Iima [1] によって提案された．彼らはその手法を生物対流のパターン画像に対して適用し [2]，実験データと定性的に一致するようなシミュレーション結果を出力する CA モデルを構築している．また，その方法論をある種の偏微分方程式（PDE）モデルに適用し，その数値解のデータセットを用いて，その解の振る舞いと似たパ
ターンを表す簡易 CA モデルを構成しようと試みた[3]．しかし，PDE の数値解を用いる場合は，得られる局所規則はわずかな時間ステップの間に定常状態に陥るパターンを与えることが多く（すなわち CAのパターンが時間とともに変化せず固定される），現象の模倣には程遠いのが現状であった．彼らは確率的 CA を導入することでその問題点を突破したが，しかし決定論的 CA モデリングに関しては依然としてパターンの固定に対する解決策が与えられていなかった．
それに対して，我々は時空間スケールを CA のパターン伝播に合わせることでその問題を解決した [5]．

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

1
0.8
0.6
0.4
0.2

0

50

100

150

図 1: データセット U

(j)

200

250

300

350

5

200

4

300

3

400

2

500

1

400

方程式の初期値の例.
1
0.8

0.4

100

150

200

250

300

350

400

(b) PDE simulation: M= 7

0

6
5

200

4

300

3

400

2

500

1
0

50

100

150

200

250

300

350

400

0

図 3: (a) 本手法によって得られた統計的 CA によるシミュレーション．CA の初期値は (5) を状態数

0.2

M = 7 で (3) によって離散化した． (b) 同じ初期値
0

50

100

150

200

250

300

350

に対する Burgers 方程式の数値解を状態数 M = 7 で400

図 2: CA シミュレーションのための初期波形: 赤線
は (5), 緑線は (6), 青線は (7) で与えられる.
そこでは Burgers 方程式の数値解を対象として，その数値解から改良された統計的 CA 構成法を用いてCA モデルを構成している．その得られた CA モデルは元の Burgers 方程式の持つ粘性衝撃波面の伝播などの諸性質を良く再現できている．さらにこの経験的 CA モデル構成法は数値解析的なアプローチも持ち，対象とする PDE との関連性も見出されており [4,6]，今後の発展が非常に期待される研究である．
本稿では [5] での研究について概略を説明する．

2

50

100

600

0.6

0

0

0

を構成するための Burgers

6

100

600

0

(a) CA simulation: M= 70

u(t, 0) = 1, u(t, L) = 0 (0 < t < T ),
 u(0, x) = u (x) (0 < x < L)

与えると，x の正の方向に波が伝播し，さらにその波面が衝撃波面を形成することが良く知られている．
波面の伝播速度はその面の両側の u の値によって定まり，速い波面が遅い波面を吸収する様子が観察される．
統計的 CA 構成法のため，様々な初期値に対してBurgers 方程式の数値解を構成する必要がある．こ
こでは適当な離散スキームの下，j 番目の初期値に対する数値解 ui,n = u(j) (i∆t, n∆x), i = 0, 1, . . . , I,(j)n = 1, 2, . . . N を以下のように構成する．上記の波面伝播の性質から，初期値 u0,n の値を，n = 0 のとき

データ
ここで用いるデータは Burgers 方程式

ut + (u2 /2)x = νuxx
(0 < t < T, 0 < x < L),(3) によって離散化したプロット．(a) と (b) 両方とも横軸と縦軸はそれぞれ空間方向と時間方向を表す．

1, n = 1, 2, . . . , N/2 − 1 のとき [0, 1] 上の一様乱数,
n = N/2, N/2 + 1, . . . , N のとき 0 を与え（図 1），
数値解を生成する．ここでは J = 100 セット作成し，
以下のようなデータセットを得る．

(1)0の数値解データである．ここでは，L = 400, T = 600,

Re = 10 としている．u0 ≥ 0 を満たす初期条件を

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

U (j) := {ui,n | i = 0, 1, . . . , I − 1,

U :=

J
∪

j=1

n = 1, 2, . . . , N − 1},
U (j) .

(2)

3(a) CA simulation: M= 7

0

6

100

5

200

4

300

3

400

2

500

1

600

0

50

100

150

200

250

300

350

400

(b) PDE simulation: M= 7

0

0

6

100

5

200

4

300

3

400

2

500

1

600

0

0

50

100

150

200

250

300

350

400

手法
統計的 CA 構成法を以下のようにしてデータセッ

ト (2) に適用する．この手法は，所与のデータセッ
トについて状態変化の最大頻度を与えるパターンを
局所ルールとして採用する方法である．

CA の状態数を M ∈ N, M ≥ 2 とするとき，
�
{(
)
�
(j)
(j)
(j)
Ua,b,c =[ui,n−1 ], [ui,n ], [ui,n+1 ] = (a, b, c)�
(j)ui,n ∈ U}

とおく．ただし，a, b, c = 0, 1, . . . , M − 1 である．ここで，[x] は x の状態数であり，[x] = floor(xM )(3)で与えられる．この Ua,b,c は，隣接 3 セルの状態数が (a, b, c) と並んでいる三つ組みのデータ全体からなる集合である．さらに，その三つ組みのデータの

図 4: (a) 本手法によって得られた統計的 CA によるシミュレーション．(b) Burgers 方程式の数値解．
用いた初期値は (6) である．離散化と軸はいずれも
図 3 と同様である．

次の時間ステップでの行き先を
�
{
� (j)
(j)
Ua,b,c (d) = ui+1,n � [ui+1,n ] = d,
(
)
}
(j)
(j)
(j)
ui,n−1 , ui,n , ui,n+1 ∈ Ua,b,cと定める．すると，隣接 3 セルが (a, b, c) となるよ6

うな三つ組みのデータの次時間ステップでの状態の

100

5

最頻値は

200

4

300

3

によって与えられる．ここで，# は集合の要素数を4002示す．この計算をあらゆる (a, b, c) に対して適用す5001ることで，M 状態 CA の全ての局所規則を決定する

600

0

ことができる．

(a) CA simulation: M= 7

0

0

50

100

150

200

250

300

350

400

(b) PDE simulation: M= 7

0

6

100

5

200

4

300

3

400

2

500

1

600

0

0

50

100

150

200

250

300

350

400

d = arg max#Ũa,b,c (d′ )
′
d

(4)

尚，この手法ではあらゆる数値解の変動パターンを記録するために非常に多くのデータを用いている．
また，状態数 M における決定するべき局所規則の数が M 3 となるため，M が増大すると計算負荷が高い．しかし，局所規則の計算アルゴリズムは並列化
計算と相性が良く，スーパーコンピュータの利用意
義の高い研究となっている．
さらに，上記の手法によって構成した統計的 CAの局所規則 (4) を用いてシミュレーションをおこな

図 5: (a) 本手法によって得られた統計的 CA によるシミュレーション．(b) Burgers 方程式の数値解．
用いた初期値は (7) である．離散化と軸はいずれも
図 3 と同様である．

い，得られた CA モデルの再現性を確認する．ここでは 3 つのケースに対して確認する．

1 つ目は，シミュレーションの初期値をシフトされたシグモイド函数を用いる：

u0 (x) =tanh((80 − x)/4) + 1(0 ≤ x ≤ L).
2(5)

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

この初期波面の傾きは空間スケールに対して相対的た CA モデルは良く Burgers 方程式の持つ性質を再にきつくなっている（図 2 赤線）．

現できていると言える．

2 つ目は，階段函数
(0 ≤ x ≤ L/4),
 1
u0 (x) =1/2 (L/4 < x ≤ L/2),

0(L/2 < x ≤ L),

初期値を (7) とした場合は, Burgers 方程式の数値
解は値の傾斜が徐々に切り立ち，粘性衝撃波面を構(6)である（図 2 緑線）．これは不連続点を持つため衝撃波面が初期から存在するケースであるが，方程式に粘性があるために数値シミュレーションは破綻せず進むことができる．

3 つ目は，区分的に直線の函数である（図 2 青線）．
{
1 − 2x/L (0 ≤ x ≤ L/2),
u0 (x) =(7)
0(L/2 < x ≤ L).
(5) と比べると傾きは緩くなっている．
また，統計的 CA 構成法によって構成された CAモデルの再現性の比較のため，ここでは上の (5), (6),(7) と同じ初期値に対して Burgers 方程式の数値解も計算する．

4

結果
ここでは M = 7 のときの結果について述べる．

J = 100 セットのデータを用いると，全ての状態の三つ組み (a, b, c), a, b, c ∈ {0, 1, . . . , M − 1} に対して，局所規則 (4) が決定可能であった．

図 3 に本手法で得られた統計的 CA モデルによるシミュレーションと Burgers 方程式の数値解をプロットした (5). どちらの図も比較のために (3) の方法
で離散化している（図 4, 5 でも同様である）．図が示すように，両方とも波面（白と黒の境界面）の伝播速度は定量的な一致をみせている．また，統計的CA の結果は，状態数 6 の領域が散逸によって消えたり，波面が停止したりすることはない．これは既存の手法では達成できなかったことである．
初期値 (6) の場合の結果を図 4 に表す．いずれも伝播の速い波面（状態数 6 と 3 の境界面）が遅い波面（3 と 0 の境界面）に到達することが見て取れる．
さらに，2 つの波面が合体した後は，波面の速度が早まっていることも一致している．これによって，統計的 CA シミュレーションでは 3 つの異なる波面の伝播速度を表現できるモデルとなっていることがわ
かる．そしていずれの速度も元の Burgers 方程式の伝播速度とほとんど一致しており，統計的に構成し28 − Academic Center for Computing and Media Studies, Kyoto University

成しようとする（図 5b）．統計的 CA シミュレーショ
ンの場合は，時間ステップが経っても状態数 0 と 1（白と黄）の境界面と 1 と 2（黄と橙）の境界面が動かないが，状態数 3, 4 と 5 の境界面はそれぞれ異なる速度で左方向へ動き，やがて 1 つの波面へと合体する．波面の固定が起きているために CA モデルはBurgers 方程式の解の挙動を完全に模倣できているとは言えないが，しかし依然としていくつもの波面の伝播速度を表現できることによって，Burgers 方
程式の解の漸近状態と良く似た挙動を最終的には取ることがわかり，良い一致を与えていると言える．
参考文献[1] A. Kawaharada and M. Iima. Constructing cellular automaton models from observation data. In
2013 First International Symposium on Computing
and Networking, pp 559–562 (2013).
[2] A. Kawaharada and M. Iima. An applicationof data-based construction method of cellular automata to physical phenomena. The special issueof Journal of Cellular Automata, to appear.
[3] A. Kawaharada, E. Shoji, H. Nishimori, A. Awazu,
S. Izumi and M. Iima. Cellular automata automatically constructed from a bioconvection pattern. Recent Advances in Natural Computing, Mathematics
for Industry, 14 pp. 15–25 (2015).
[4] A. Kawaharada, T. Miyaji, and N. Nakano. An
analyzable method for constructing a cellular automaton from a continuous system. In 2015 Third
International Symposium on Computing and Networking, pp 418–423 (2015).
[5] A. Kawaharada, T. Miyaji, and N. Nakano. Properchoice of spatio-temporal scale and dataset subsampling for empirical CA construction. In 2015 Third
International Symposium on Computing and Networking, pp 424–429 (2015).
[6] A. Kawaharada, T. Miyaji, and N. Nakano. Analysis of a method for constructing a cellular automaton from a continuous system. Int. J. Networking
and Computing 6(2), pp 230–242 (2016).

共有結合性有機骨格構造体(COF)の電荷移動シミュレーション
鬼頭(西岡)宏任
名古屋大学大学院理学研究科 (現：筑波大学計算科学研究センター)

1 INTRODUCTION

考えられる。そこでまず、TP-COF の３×３×４
スパーセルに対し、周期境界条件と室温の NVT

共有結合性有機骨格構造体(COF)と呼ばれる多アンサンブルを課した古典 MD シミュレーショ

孔性高分子は、その構造特性(自己集合、熱安定性、ンを実行した。次に、MD シミュレーションのス

軽量)と機能特性(ガス吸着貯蔵、不均一性触媒、ナップショットに対して量子化学計算を実行し、光伝導)から大きく注目されている[1]。ポルフィ

ピレンユニット間のホール移動積分を求めた。

リンやフタロシアニンなどの芳香族分子が連結して作られた２次元シート形状の COF は、シート間のπ-π相互作用によってグラファイトのように積層する。半導体性特性を示す COF では、その骨格に使われた芳香族分子のπ軌道を使って、二次元面に垂直なカラム方向に効率よく電荷が移動していると考えられているが、その詳しい実態とメカニズムは良く分かっていない。
本稿では、古典分子動力学(MD)シミュレーショ
ン、量子化学計算、キャリア移動シミュレーショ
ンを組み合わせたマルチスケールシミュレーショ
ンを用いて、COF 中の電荷輸送の理論解析を行った研究成果[2]について簡単に紹介する。

図 1(a)TP-COF (b)移動積分計算結果の分布(c)移動積分データのパワースペクトル2 移動積分の揺らぎ

移動積分の計算には、フラグメント分子軌道
本研究[2]では、p 型半導体特性を示すことが実験的に確認されている TP-COF[3]を対象とした。
図 1(a)に示すように、TP-COF はトリフェニレンとピレンが連結した COF で、イオン化ポテンシ(FMO)法から電荷移動パラメータを計算する手法を用いた[4]。今回は、各フラグメントの電子状
態計算を密度汎関数(DFT)法から DFT の精度を再現する半経験的手法である密度汎関数強束縛ャルエネルギーの低いピレンユニットの方を通っ(DFTB)法に置き換えた FMO 法[5]を用いることて、ホールキャリアがカラム方向に移動しているで、更なる計算コストの削減が実現されている。

と考えられる。このホール移動のシミュレーショ

図 1(b)には、60 フェムト秒毎に約 492 ピコ秒

ンを実行するためには、カラム方向に隣接したピ

間、移動積分を計算した結果の分布をプロットしレンユニット間のホール移動積分値をパラメータている。移動積分の平均値と揺らぎの標準偏差のとして見積もる必要がある。その際、シート間は

結果は、38meV と 102 meV であり、構造の熱揺

π-π相互作用によってスタックしているため、移らぎによって移動積分が大きく揺らいでいること

動積分は構造の熱揺らぎの影響を大きく受けるとが分かる。またガウス関数(図 1(b)の実線)でよく京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 29

フィットすることができるため、移動積分の揺らョン[6]を実行する。図２にはキャリアの平均二乗

ぎは正規分布に従っていると仮定できる。

変位の時間依存性をプロットしている。この結果

図 1(c)には、移動積分時系列データに対してパは直線で良くフィットすることができ、その傾き

ワースペクトルを求めた結果である。スペクトルが拡散係数となる。このシミュレーションから、の強度から、65cm−1 以下の遅いモードのみが移動

1.93cm2V−1s−1 という非常に高いキャリア移動度

積分の揺らぎに寄与していることが分かる。
また、が TP-COF に対して見積もられた。TP-COF に

25 cm−1 辺りに顕著なピークが現れる。もう一つ対するキャリア移動度の実験値は報告されていなの顕著なピーク(5 cm−1 辺り)は非常に遅いモード

いが、我々の計算値は、別の p 型半導体 COF のであることと、移動積分が正規分布に従うことか

実験値(>1.3 cm2 V−1s−1)[1]と良く一致している。

ら、移動積分の揺らぎの性質は 25 cm−1 の振動数

今後、今回の手法を、より良い半導体特性を示すをもつ調和振動子を使って再現できる。

COF の理論設計へ応用することを目指していく。

3 電荷移動シミュレーション
TP-COF におけるホールキャリア電荷移動シ
ミュレーションは、上で述べた移動積分の揺らぎの効果を取り入れる必要がある。そこで以下の一次元強束縛ハミルトニアン[6]を用いて、ピレンを通ってカラム方向にホールキャリアが移動していくことをモデル化する。
H = ∑α x j,1 j jj+∑⎡⎣−τ + β ( x j+1,2 − x j,2 )⎤⎦( j j +1 + j +1 j )
j

⎛1
2⎞
1
+∑⎜ m1v 2j,1 + m1 (ω1 x j,1 ) ⎟
⎝
⎠
2
2
j

⎛1
2⎞
1
+∑⎜ m2 v 2j,2 + m2 (ω 2 x j,2 ) ⎟
⎝2
⎠
2
j

図 2 キャリアの拡散シミュレーションから得られた平均二乗変位の時間発展

4 ACKNOWLDGEMENTS
本研究では、京都大学学術情報メディアセンタ

ここで、サイト j がカラム方向のある一つのピレンを表す。右辺第二項では、ω2=25 cm−1 の固有振動数を持つ調和振動子 xj,2(第 4 項)とキャリアが線形に相互作用することで、図 1(b)(c)に示した移動積分の揺らぎが引き起こされると仮定している。
(その他のパラメータについては、文献[2][6]で説明されている)。
キャリア移動度は、キャリアの拡散係数 D から
アインシュタインの関係式𝜇𝜇 = 𝑒𝑒𝑒𝑒/(𝑘𝑘! 𝑇𝑇)を用い

ーのスーパーコンピュータ共同研究制度(若手奨励枠)を活用させて頂きました。
[1]X. Feng, et al., Chem. Soc. Rev. 41, 6010
(2012)[2]H. Kitoh-Nishioka, K. Welke, Y. Nishimoto,
D. G. Fedorov, and S. Irle (in preparation)[3] S. Wan, et al., Angew. Chem. Int. Ed. 47,
8826 (2008).
[4] H. Kitoh-Nishioka, K. Ando, Chem. Phys.

て見積もることができる。そこで、一次元系の中Lett. 621, 96 (2015); 鬼頭(西岡) 宏任, 安藤耕司,

心サイトにホールの波動関数が局在していると初京都大学 学術情報メディアセンター 全国共同利期条件を設定し、波動関数がどのように拡散する用版［広報］Vol.14 No.2 pp. 8-9. (2016)

かをシミュレーションする。その際、キャリア波[5] Y. Nishimoto, et al. J. Chem. Theory

動関数の時間発展は時間依存シュレディンガー方Comput. 10, 4801 (2014); Y. Nishimoto, et al. J.

程式を、調和振動子の運動の時間発展は古典ニュ

Chem. Phys. Lett. 6, 5034 (2015).

ートン方程式を解く、量子古典混合シミュレーシ

30 − Academic Center for Computing and Media Studies, Kyoto University[6] A. Troisi, Adv. Mater. 19, 2000 (2007)

理想化モデルを用いた, 竜巻を模した渦の力学に関する研究
松嶋 俊樹 ∗∗

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

はじめに

1

を, C ∞ 級条件, エネルギー有界条件を満たすように発達した竜巻は, 渦の中心に下降流をもつような,

2 セル構造を持つことが知られている. そのような構造は, 不安定性によって, 渦の鉛直軸方向に対して左手回りの螺旋構造をもったモードや, 鉛直方向にほぼ一様なモードが発達し [1], 竜巻の親渦の中に複数の子渦が含まれる多重渦構造との関連が示唆されている. そのような現象の力学を理解することを目指する. 非圧縮性 (2) により, 速度場はトロイダルポテ
ンシャル ψ とポロイダルポテンシャル χ によって以下のように表される.

u = ∇ × (ψez ) + ∇ × ∇ × (χez ),(3)

半無限領域スペクトル法による数値モデルは, ψ, χ, ϕ =∇2 χ についての方程式を基礎として構築する.

し, 著者らは, 竜巻を模した渦の数値計算のための,平らな平面の上の半無限領域において, 下端で粘着境界条件を課すことが可能なスペクトル法モデルの開発を行ってきた [2, 3]. 本課題では, そのような数値モデルの検証実験として行った渦輪の壁衝突の数
値実験の結果を解析した. 本報告では, はじめに開発した数値モデルについて簡単に説明し, 次に行った解析について説明する.

2

数値モデルの概要

2.1

支配方程式

2.2

離散化

問題領域は半無限であるため, 空間離散化には無限領域のスペクトル法を用いる. まず, 水平方向に対して, 円柱座標系 (r, θ) における無限平面を, 以下のように射影する.
√
1+µr = 2R,1−µr ∈ [0, ∞), µ ∈ [−1, 1].

(4)

式 (4) を介して, 変数 ϕ は以下のように球面調和関数展開する.

ϕ=

三次元非圧縮ナヴィエ・ストークス方程式は, ベク
トル不変形で,

M
∑

N
∑

m
imθ
.
b̃m
n (z, t)Pn (µ)e(5)

m=−M n=|m|

係数 b̃mn (z, t) に対して, さらに以下のような展開を
2

|u|
∂u
+∇
− u × ω = −∇P − ν∇ × ω,
∂t
2

(1)

∇ · u = 0.

(2)

ここで, u は速度場, ω は三次元渦度場で, ω = ∇×u,

P は順圧圧力関数, t は時刻, ν は動粘性係数である.
以上の方程式を, 円柱座標系の領域 r ∈ [0, ∞), θ ∈[0, 2π], z ∈ [0, ∞) において z = 0 での粘着条件を課して解く. また, 基底関数によって展開された速度場考える. 各 m(m = 0 については省略) についての五重対角対称行列 Qm
ij を, その (i, j) 成分が

Qm
ij =

∫
√
1 1
i(i + 1)j(j + 1)
(1 − µ)2 Pim (µ)Pjm (µ) dµ,2 −1(m ̸= 0; i, j = |m|, |m| + 1, 2, . . . , N ),(6)となるように導入する. 行列 Qm
ij の正規直交固有ベmクトルと対応する固有値をそれぞれ vm
n , λn (n =1, 2, . . . , N − 1) と書き, 固有ベクトルの成分を vm
n =京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 31m(vin) (i = 1, 2, . . . , N − 1) と表す. これによって, 係m
数 b̃n (z, t) を以下のように固有ベクトルで展開する.

ϕ̃m
p (z, t) =タに対してブロックサイクリック分割を採用し, 実
空間では各動径方向の分点のインデックス j のデー
タに対してブロック分割を採用した.

N √
∑

m
n(n + 1)b̃m
n vnp

n=|m|

(p = 1, . . . , N − |m| + 1).

(7)以上の展開を用いると, 下端の粘着境界条件の代わりに, 以下の ϕ の積分制約条件が得られる.
( √
)
∫ ∞
m
λ
p
z dz = 0ϕ̃mp (z, t) exp −
4R2
0

3.2

解析と可視化

渦の時間発展を追跡するため, 粒子を用いた解析を行う. 時刻 t = 1 間隔で出力された三次元速度場と渦度場のデータ (ただし, 領域 r ≤ 2, z ≤ 2 に限定)

(8)

鉛直方向の離散化には, 水平方向と類似の以下の射影を用いて, 初期に空間的にランダムに配置した 20000
個の粒子を移流させる. さらに, 移流させた各粒子から出発し, 渦度ベクトルの方向へと辿ることによって, 渦線を構成する. 渦線をこのように構成すると,非粘性の場合, 流体力学に関するヘルムホルツの第

z=H

1+x
,
1−x

z ∈ [0, ∞), x ∈ [−1, 1].

(9)二法則より, 各渦線が時刻に依存せず常に同じ流体粒子群から構成される. 今の場合は粘性流体を考えを行った後, チェビシェフ多項式を用いる. 詳細は省ているため, 渦線を構成する粒子は入れ替わるが, 少略するが, 論文 [3] では, 鉛直方向の展開を与えていなくとも渦線の一点で束縛される. したがって, 時間

る. また, 粘性項に陰的解法を用いて時間積分する際

局所的には, 渦線をトレーサーと見て, 渦の追跡を行の精度を改良するための反復手法を導入し, 時間離うことができる. 粒子の移流計算と渦線の構成には,散化した偏微分方程式が高精度に解くことができる

ParaView (http://www.paraview.org/) をシステム
D の CRAY XC30 向けにクロスコンパイルして用いた.
初期の渦線と y = tan(18◦ )x 面の方位角渦度を図 1に示した. 初期の渦輪はほぼ軸対称な分布であるが,方位角波数 5 の倍数において, 有限振幅の擾乱を入れている. 図は省略するが, 時間が経つと, 渦輪は自己誘導速度で壁に衝突し, 二次渦輪が生成される. その後,二次渦輪は初期に与えた主渦輪の周りを回りこんだ後, θ = 18◦ + 72◦ × i(i = 0, 1, . . . , 4) のまわりで主渦輪の下に潜り込み, θ = 36◦ + 72◦ × i(i = 0, 1, . . . , 4)
では中心軸へ移流されて主渦輪から離れる. 結果として, 図 2 のように花びら状の構造を作りつつ, 流れは複雑になっていく.

ことを示している.

3

結果

3.1

数値実験

数値モデルのパフォーマンスを評価するため, 三次
元的に屈曲した渦輪の水平な壁との衝突の数値実験を行う. 必要な計算資源量を削減するため, 球面調和
関数展開において, 方位角波数の 5 の倍数のモードだ
けを考慮するセクターモデルを用いる. 水平方向に用いた写像における球の半径を R = 0.55, ルジャンド
ル多項式およびルジャンドル陪関数展開の切断次数をN = 745, 分点数を J = 750, 方位角方向の分点数を
2J/5 = 300 としている. 鉛直方向に用いた写像のパラメータを H = 1.2, 切断波数を L = 633, 分点数をK = 640, 時間積分の時間ステップを ∆t = 0.0125,
動粘性係数を ν = 1.45 × 10−5 とする. 粘性項と非線形項の時間離散化にはそれぞれ BDF2 と BX2* を用いた. また, 境界条件は不透過条件と粘着境界 (すなわち積分制約条件) を課した. 数値コードは MPIを用いて並列化した. このとき, 各プロセスでのデータの分割において, 波数空間では各方位角波数のデー

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

4

まとめ開発した半無限領域のスペクトル法モデルを用いた数値実験として渦輪の壁衝突の数値実験を行い, その解析を行った. これらの詳細について記述した論文 [3] については現在印刷中である. 今後は, 開発した半無限領域のスペクトル法モデルを用いて, 竜巻を模した渦の三次元の数値実験を行う予定である.

図 1: 時刻 t = 0 における渦度の大きさを色で示した三次元渦線と, y = tan(18◦ )x 面の方位角渦度.
図の中の黒点は, (r, θ, z) = (2, 2πi/1500, 0) (i =0, 1, . . . , 1499) を示している.

図 2: 時刻 t = 200 であることを除き図 1 と同じ.
謝辞

本課題は, 京都大学学術情報メディアセン

ターの共同研究制度 (若手・女性研究者奨励枠) と,特別研究員奨励費 (No. 15J01380) の助成を受けた
ものです.

参考文献
[1] 松嶋 俊樹・石岡 圭一, 2014: サリバン渦の安定性について, 日本流体力学会誌ながれ 33, 47-60.
[2] 松嶋 俊樹・石岡 圭一, 2016: 半無限領域のスペクトル法による渦輪の壁衝突の高解像度計算, 日本流体力学会誌ながれ 35, 89-92.
[3] Matsushima, T. and Ishioka, K., 2016: A highly
accurate spectral method for the Navier-Stokes
equations in a semi-infinite domain with flexible
boundary conditions, Fluid Dynamics Research, inpress

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

オフロード路面における車輪の走行性能解析
惑星ローバ用車輪におけるラグ周辺部分のミクロな砂の挙動解析
西山健太 ∗∗

1

京都大学農学研究科地域環境科学専攻

緒言

解析条件

3

車輪と土の相互作用問題は，様々な要因が介在す

ドイツ航空宇宙センター (DLR) で実施された，珪る複雑な研究課題であり，これまでは実験的あるい

砂を充填した土槽での火星ローバ用車輪の走行試験

は半経験的手法により解析されてきた. 近年は数値解[2] を参考に解析条件およびパラメータを決定した．

析手法による車輪の走行予測が試みられている．高

供試車輪は，B3 車輪と呼ばれる弾性車輪と，B6 車

精度な走行予測は，エネルギー効率の高い車輪の開輪と呼ばれる剛性車輪の 2 種類とした．車輪は直径

発や，車両による車輪の駆動力制御手法の確立のた
本研究の目的はオフロード路面における車輪の走300 mm，ラグは 12 個のタブルグローサが等間隔に配置されラグ高さ 8 mm，厚さ 1.2 mm，2 つのラグの間隔は 10 mm である．また，DEM の円要素の直行を高精度にシミュレーションすることである．特に径は 3, 2.5, 2 mm (混合個数比 1:2:3)，要素数 47453

ラグを有する惑星ローバ用の金属製弾性車輪のけん

個，時間刻み 1.25×10−6 s とし解析を実施した．なめに必要と考えられる．

引性能について，ラグ間やラグ近傍の応力状態を把お，DEM 要素直径はパソコンによる解析 [1] で用い握してミクロな観点から検討することを目指す． こた直径の半分とした．

れはラグの数・幅・長さ等のラグ形状を決定する際に，車輪全体に発生するマクロな走行性能の観点から議論する既往の研究では不十分なことによる．

2

解析手法
車輪を有限要素法 (FEM)，土を離散要素法 (DEM)

により 2 次元でモデル化し，陽解法による数値積分を用いて解析を行う．車輪にけん引負荷とトルクを加えることで車輪の走行が制御され，1 回の走行中は一定の滑り率に維持される．また，走行性能は定常状態で走行している区間の平均を取り，滑り率に対する推進力，けん引力，走行抵抗，沈下量で評価される [1]．なお，ローバ用の車輪の多くは円筒形状であり 3 次元的なラグパターンを持たず，2 次元解
析でも十分有用な結果が得られると考えている．

4

解析結果と考察4.1

走行後の状態

図 1 に走行後の車輪と土の変形挙動の結果例を示
す．ラグの間に入り込む粒子数が従来より増加し，より詳細な変形状況が表現できている．また，パソコ
ンによる解析では見られなかったが，轍に車輪のラグによる周期的な跡が明確に残されていることが確認できる．そのため，轍の形成については，ラグ形状との関連で要素径を細かくする必要性がうかがえる．

4.2

走行性能の結果

マクロな車輪全体に働く走行性能の実験結果を図2，解析結果を図 3 に示す．B3 ( 弾性 ) 車輪の方がけん引力が高いなど，両者の間で類似した傾向が見受けられ，妥当な解析が実施できたと考えられる．

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

図 1: 走行後の車輪と土の変形挙動 ( B6 車輪，滑り率 60 % )τ xy [kPa]
H, B3
P, B3
R, B3
H, B6
P, B6
R, B6

Force [N]

100

50

40

0

0
-40

-50

図 4: せん断応力 τxy の例 ( B6 車輪，滑り率 60 % )010 20 30 40 50 60 70
Slip [%]

図 2: 実験結果 ( 推進力 H ，けん引力 P ，走行抵抗

4.3

応力状態

図 4 に応力解析の例を示す．ラグを起点に車輪後R ) [2]

方へと土中に局所的な応力鎖が発生していることが明らかである．

H, B3
P, B3
R, B3
H, B6
P, B6
R, B6

Force [N]

100

50

0

10 20 30 40 50 60 70
Slip [%]

謝辞

図 3: 解析結果 ( 推進力 H ，けん引力 P ，走行抵抗

R)[1] K. Nishiyama, H. Nakashima, T. Yoshida, T. Ono,
H. Shimizu, J. Miyasaka, K. Ohdoi. 2D FE-DEM
analysis of tractive performance of an elastic wheel
for planetary rovers．Journal of Terramechanics 64(2016) 23-35．
[2] M. Scharringhausen, D. Beermann, O. Kromer,
L. Richter．Single wheel tests for planetary applications at DLE Bremen．11th European Regional Conference of the International Society for
Terrain-Vehicle Systems – Bremen, October 5-8
2009．

-50
0

参考文献

本研究は，スーパーコンピュータ共同研究制度 (若手・女性研究者奨励枠) によるものである．記して関係各位に厚く御礼申し上げる．

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

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

非線形有限要素法による鉄筋コンクリートの大規模破壊シミュレーション
車谷 麻緒 ∗∗

1

茨城大学 工学部 都市システム工学科

はじめに
土木構造物の設計や維持管理の高度化において，鉄筋コンクリートの力学性状の評価は重要である．鉄筋コンクリートは材料挙動のまったく異なるコンク
リートと鉄筋で構成されており，補強材である鉄筋の配置（配筋）によって異なる破壊モードを示す．鉄筋コンクリートの力学性状を評価するには，材料のモデル化やパラメータを変えずに，配筋の違いのみを考慮して，様々な破壊モードを再現できる数値解
析手法が必要になる．
配筋の違いのみを考慮して，鉄筋コンクリートの様々な破壊モードを再現するには，鉄筋コンクリー
トの平均的な挙動を考えるのではなく，コンクリー

らに，領域分割と MPI 並列化を行うことにより，大規模シミュレーションに応用可能な方法に拡張する．
提案手法の有効性と妥当性を検証するため，せん断補強筋 10 本を有する RC はりの 4 点曲げ試験に本手法を適用し，メッシュサイズによる影響と実験結果の再現性について検討する．

2

数値解析手法
本研究では，鉄筋コンクリートの 3 次元破壊挙動を詳細にシミュレートするために，コンクリートと鉄筋の材料挙動を別々にモデル化する．以下に，コンクリートと鉄筋の材料モデルを簡単に示す．

トのひび割れ進展挙動と鉄筋の塑性変形を別々にモデル化する必要がある．さらに，配筋などの内部構造のみを考慮して，鉄筋コンクリートの破壊挙動を再現するには，メッシュサイズの小さい有限要素モ
デルが必要になるため，メッシュサイズに依存しな

2.1

コンクリートの材料モデル

コンクリートのひび割れ進展挙動のモデル化には，コンクリートの破壊力学に基づく修正 von-Mises 損

いモデル化の適用と解析プログラムの高度化および傷モデル [1] を適用する．このモデルは，コンクリー

数値解析の高速化が必須となる．近年の計算機ハートの破壊進行過程における結合力と開口変位の関係

ドウェアの急速な発達により，大規模な数値シミュを連続体損傷モデルに応用したもので，コンクリー

レーションを行える環境は整いつつある．しかし，そトの破壊力学の特徴である破壊エネルギーの定義ものハードウェア環境を利用して，鉄筋コンクリート

満足するモデルである．これを通常の有限要素法にの破壊シミュレーションを高精度に行える解析手法

導入することにより，コンクリートのひび割れ進展が確立されているとは言い難い．

挙動を再現する．

そこで本研究では，コンクリートのひび割れ進展

多次元での損傷を評価するに際して，ひずみテン

挙動のモデル化にコンクリートの破壊力学に基づくソルの不変量を用いた等価ひずみを用いる．具体的

変形のモデル化に von-Mises 塑性モデルを適用する

く等価ひずみ [3] を採用する．

修正 von-Mises 損傷モデル [1] を適用し，鉄筋の塑性ことにより，鉄筋コンクリートの破壊挙動を 3 次元で詳細に再現できる数値解析手法を構築する [2]．さ36 − Academic Center for Computing and Media Studies, Kyoto Universityには，次式で与えられる修正 von-Mises 基準に基づ

σ2

σ2

Saturation

3D view(Plane strain condition)

Q

σ y0

n

E

σ y0

σ1

Principal strain ε 2

σ y0
σ y0

0

σ y0
ε2

Expanded yield surface
Initial yield surface

k=1

図 2: 非線形等方硬化則に基づく von-Mises 塑性モ
デル
k=2

2.2
0

鉄筋の材料モデル

本研究では，鉄筋の構成モデルに，図 2 に示すよPrincipal strain ε 1

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

図 1: 2 次元および 3 次元での主ひずみ空間におけるデルを適用する．
次に示すように，全ひずみは弾性ひずみと塑性ひ

等価ひずみの等値面ずみに加算分解され，応力はそのうちの弾性成分の

εe =

k−1
1
I′ +
2k(1 − 2ν) 1 2k

√

(

k−1 ′
I
1 − 2ν 1)2

みが寄与するものとする．

+

12k
J′
(1 + ν)2 2(1)

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

はひずみテンソルの第 1 不変量，J2′ は偏差ひずみテ
ンソルの第 2 不変量である．

3 次元および 2 次元の主ひずみ空間に εe の等値面をプロットした例を図 1 に示す．圧縮引張強度比 kε = εe + εp()
σ = c : εe = c : ε − ε pテンソル，εp は塑性ひずみテンソル，σ はコーシー
応力テンソル，c は弾性係数テンソルである．

図 2 に示すような，指数関数で与えられる非線形

等方硬化則を採用することにより，降伏関数 f ≤ 0
は次式で与えられる．

)(
f = σv − σy0 − Q 1 − e−bp

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

(5)

ここで，ε は微小ひずみテンソル，εe は弾性ひずみを用いて引張破壊と圧縮破壊を考慮したことにより，性を表現できることが分かる．

(4)

(6)

ここで，σv は von-Mises の相当応力，σy0 は初期降伏応力， p は相当塑性ひずみ，Q と b は塑性硬化を表すためのパラメータである．塑性流動則は関連流れ則に従うものとし，次式で与えられる．

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

()
Eκ0 heκ0(κ − κ0 )D(κ) = 1 − exp −
κ
Gf

(3)

ここで，κ0 は損傷開始時の等価ひずみである．

dεp = dp

3

∂f
= dpn
∂σ

(7)

数値解析例
本解析手法においては，損傷モデルを用いて，コンクリートのひび割れ進展挙動をモデル化している．
既往の損傷モデルを用いた破壊シミュレーションで京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 37

450
400

Load (kN)

450
100: Experiment: FEM (500,000x4 elements): FEM (1,000,000x4 elements): FEM (2,000,000x4 elements)

φ6

y

200 160

x
z

φ9
0

D16

0

150

Material parameter of steel
E

ν

σ y0

Q

200 GPa

0.3

380 MPa

150 MPaE
29 GPab
15

Material parameter of concrete
κ0
Gf
ν

k

0.2

20

0.0001

0.1 N/mm10Displacement (mm)

図 3: せん断補強筋 10 本を有する RC はりの解析条
件と材料パラメータ

図 4: メッシュサイズの異なる有限要素モデルに対する荷重－変位関係の比較本研究では，異形鉄筋の節が表現できる程度の小さなメッシュサイズを設定し，異形鉄筋の幾何形状を詳細にモデル化している．そのため，コンクリー
トと鉄筋の境界面については，特殊な界面要素や付着モデルを用いずに，境界面に位置する有限要素の応力集中によって，境界面から損傷が発生し，付着は，解析結果がメッシュの寸法や分割に依存することが知られており，これを回避するために特別なモ力が徐々に消失する過程を直接的に再現する．

さらに，図 3 と同条件の RC はり試験体を作成し，デル化を行うことが多い．筆者らの既往の研究 [1, 2]

4 点曲げ試験を行い，解析結果と実験結果とを比較

では，鉄筋の塑性変形が生じない問題を対象に，鉄する．

筋コンクリートのひび割れ進展挙動がメッシュサイ
ズに依存しないことを確認している．以下では，コンクリートのひび割れ進展と鉄筋の塑性変形が生じる RC はりの破壊シミュレーションを対象に，メッ

3.2

解析結果と考察

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

シュサイズの影響と実験結果の再現性について検討を図 4 に示す．メッシュサイズの異なる 3 ケースのする．

結果はほぼ一致しており，メッシュサイズの相違が
荷重－変位関係に与える影響は非常に小さいといえ

3.1

解析対象と条件

解析対象は，図 3 に示すような，せん断補強筋（スターラップ）を 10 本有する曲げ破壊先行型の RC はりの 4 点曲げ試験である．主鉄筋は D16 の異形鉄筋

る．解析結果は実験結果ともほぼ一致しており，本手法による破壊シミュレーションの妥当性が見て取れる．
次に，ひび割れの発生と進展の相違について検討する．モデル表面での等価ひずみの分布を比較したとし，数値解析では異形鉄筋の幾何形状を忠実に有ものを図 5 に示す．提案手法による 3 次元破壊シミュ

限要素メッシュに反映させることとする．問題の対レーションは，ひずみの局所化によりひび割れを再称性を考慮し，1/4 モデルを解析対象とする．1/4 モ

デルの要素数を約 50 万，100 万，200 万要素とし，有限要素には 1 次の四面体要素を用いる．解析条件

は，強制変位 12 mm を 500 ステップで載荷する変位制御試験とする．

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

現する方法であるので，ひずみの分布がひび割れの分布に対応する．結果を見て分かるように，メッシュ
サイズが違っていても，ひび割れの本数や経路はほとんど変化していないことが分かる．図 6 は実験結果におけるひび割れの分布を表している．図 5 と図

500,000x4 elements(a) Crack distribution at load: 50 kN

1,000,000x4 elements(b) Crack distribution at load: 100 kN

図 6: 実験結果におけるひび割れの分布

ションは，RC はりの表面から生じる損傷（ひび割

2,000,000x4 elements(a) Equivalent strain distribution at step: 40/500 (50 kN)れ）に加えて，RC はり内部の鉄筋の表面から発生・進展する損傷（ひび割れ）までを詳細に再現できる方法であることが分かる．

図 8 は，数値シミュレーションに要した計算時間を示している．3 ケースともに，領域分割を行い，128

コアによる MPI 並列を適用している．計算時間はほ500,000x4 elements

ぼ線形的に増加しており，さらに多くのコアを使用することにより，計算時間を短縮することができる．

4

1,000,000x4 elements

おわりに
本研究では，コンクリートのひび割れ進展挙動をモデル化に修正 von-Mises 損傷モデルを適用し，鉄筋の塑性変形に von-Mises 塑性モデルを適用することにより，鉄筋コンクリートの 3 次元破壊シミュレー

ション手法を構築した．さらに，領域分割と MPI 並

2,000,000x4 elements(b) Equivalent strain distribution at step: 180/500 (100 kN)
0.0001

列化を行うことにより，大規模破壊シミュレーショ
ンにも応用可能な方法に拡張した．

0.1

コンクリートひび割れ進展と鉄筋の塑性変形が顕図 5: メッシュサイズの異なる有限要素モデルに対す
るひび割れ分布の比較著に表れる RC はりの破壊シミュレーションに本手法を適用したところ，メッシュサイズに依存することなく，実験結果と同様の結果が得られることを示した．今後は，結果をまとめて，学術誌に投稿し，研6 を比較すると，数値解析から得られたひび割れ分
布は実験とほぼ同様の分布であることが分かる．

図 7 は，モデル内部における損傷（ひび割れ）の分布を 3 次元で可視化した結果である．本研究では，異形鉄筋の幾何形状を忠実に考慮してメッシュ分割を行っているため，提案手法による破壊シミュレー

究成果を発信していく予定である [4]．
謝辞：

本研究を進めるにあたり，京都大学学術情

報メディアセンターならびにクレイ・ジャパン・イ
ンクの方々には，プログラムの高度化・並列化に多大なご協力をいただきました．ここに記して，感謝の意を表します．

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

Load step: 40/500

Computing time (hour)

150

100

50

0

500,000

1,000,000

2,000,000

Number of element

Load step: 180/500(a) 500,000x4 elements

図 8: メッシュサイズの異なる有限要素モデルに対する計算時間の比較

参考文献
Load step: 40/500[1] 車谷麻緒，寺田 賢二郎，加藤準治，京谷孝史，樫山
和男，「コンクリートの破壊力学に基づく等方性損傷モデルの定式化とその性能評価」，日本計算工学会
論文集，pp.20130015，(2013).
[2] 車谷麻緒，根本優輝，相馬悠人，寺田 賢二郎，
「コン
クリートの破壊力学を考慮した鉄筋コンクリートの3 次元破壊シミュレーションとその性能評価」，日本計算工学会論文集，Vol.2016, pp.20160004, (2016).

Load step: 180/500(b) 1,000,000x4 elements[3] J.H.P. de Vree, W.A.M Brekelmans, M.A.J. van Gils,「Comparison of nonlocal approaches in continuum
damage mechanics」, Computers and Structures, 55,
pp.581–588 (1995).
[4] M. Kurumatani, Y. Soma, S. Abe, 「Fracture simulationof reinforced concrete using an isotropic damage model
based on fracture mechanics」, 12th World Congress on
Computational Mechanics, (2016).

Load step: 40/500

Load step: 180/500(c) 2,000,000x4 elements

図 7: 有限要素モデルの内部における 3 次元損傷分布40 − Academic Center for Computing and Media Studies, Kyoto University

津波流体力による構造物の破壊過程を考慮した構造流体連成解析
高瀬
∗

1

慎介 ∗

東北大学大学院

はじめに

工学研究科

し，この数学被覆と物理領域の交わりを物理被覆として定義する．物理被覆には，物理変数を近似する

災害大国である我が国では，数多くの自然災害がための被覆関数と呼ばれる未知パラメータを含む関発生し，大きな社会問題になっている．これらの自数を導入する．被覆関数には，PU 重み関数との積に然災害による被災規模を適切に把握し，被災状況をより構成される基底関数が一次独立性を満たす限り軽減する方策を事前に提案する手法として，数値シ任意の関数形を利用できるが，本研究では，定数項

ミュレーションによる被害予測が行われつつある．

のみを未知パラメータとして用いることにする．

そこで本研究では，これらの自然災害の中でも津
波漂流物および流体力による構造物の破壊に着目する．すなわち，津波流体力および漂流物の衝突による構造物の変形・破壊過程とその過程がさらに２次的な構造物被害を引き起こす災害の解析手法の開発を行う．解析手法として，流体力の計算には安定化

本研究では，FCM の物理要素の基本形状としては

1 次の四面体要素を用いることにする．そして，被覆関数は定数項のみとするので，本研究における FCMの物理要素は，必然的に FEM における 1 次の四面体

有限要素と等価なものとなる．ただし，FCM では流
体の数学領域内に構造物が存在してもよく，流体の有限被覆法，構造物の衝突計算には個別要素法を用

物理要素の境界は数学要素と独立に動くことができ

い，相互の連成解析手法を構築する．現在，安定化

る．この特徴を利用して，構造物の境界位置を正確

有限被覆法による計算は，MPI により並列計算が行に考慮したうえで流れと構造物の剛体運動の相互作われているが，衝突計算部分の個別要素法は，並列用を評価可能な数値解析手法を構築する．なお，流計算がされておらす，その部分が大規模化へのボト体と構造物の境界位置は，構造物表面をゼロとする

ルネックになっている．非並列部分の MPI による並列化が可能になれば，実際の条件に近い津波被害予距離関数である Levelset 関数作成し，流体領域の定義を行っている．

測解析が可能となり，工学的にも有効な被害予測手
法になる．

2

流体の数値解析手法

2.1

有限被覆法

2.2

流体に関する定式化

本研究で用いる支配方程式は，非圧縮粘性流体における Navier-Stokes の運動方程式と連続式で以下の有限被覆法 [1]（FCM) は，近似関数が数学領域 Ω Mと支配方程式が満たされるべき物理領域 ΩP を独立して定義する．そして，FCM では有限要素法（FEM）
における形状関数の代わりに，解析対象となる物理
領域とは独立な PU 条件（Partition of Unity) を満たす近似関数（以下，PU 重み関数）を数学被覆に導入

ように表される．
)( ∂u+ u · ∇u − f − ∇ · σ (u, p) = 0
ρ
∂t

∇·u=0

(1)
(2)

ここで，ρ は密度，u は流速ベクトル， p は圧力，fは物体力ベクトル，σ は応力テンソルである．また，京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 41

Newton 流体を仮定し，構成則には次式を用いる．
σ = −pI + 2µε(u)(3)

ここで，I は 2 階の単位テンソル，µ は粘性係数であり，ε(u) は次式で定義される変形速度テンソルである．

ε(u) =

2.3)
1(∇u + (∇u)T
2

(4)

2.4

本研究では，砕波等の複雑な自由表面形状を表現する必要性があるため，固定メッシュを用いた Euler

的手法の１つである Phase-ﬁeld 法を採用することにする．Phase-ﬁeld 法では，次式に示す保存形式に修正された Allen-Cahn 型移流方程式 [5] を解くことで自由表面位置を決定する．

∂ϕ∂t

安定化有限被覆法

3 次元流れ場の支配方程式 (1) に対して有限被覆法を用いた SUPG/PSPG 法 [2],[3] を適用すると，次式のような安定化有限被覆法による離散化方程式が得
られる．
( h
)
∫
∂u
+ uh · ∇uh − f dΩρ
wh · ρP∂t∫Ω∫+ε(wh ) : σ(uh , ph ) dΩ +qh ∇ · uh dΩ
ΩP

nel ∫
∑

ΩP

{

}

1
τsupg uh · ∇wh + τpspg ∇q
ρ
e=1
{ ( h
)
}
∂u
+ uh · ∇uh − f − ∇ · σ(uh , ph ) dΩ
·
ρ
∂t
nel ∫
∑
+
τcont ∇ · wh ρ∇ · uh dΩ
+

e=1

+

∫

ΩeP

ΩeP

h

Γg

h

p̄(u − ū ) · w = 0hPhase-ﬁeld 法を用いた界面捕捉(5)

ここで，ΩP は Navier-Stokes 方程式の物理領域，uh ，

ph は，それぞれ速度と圧力の有限要素近似式，wh ，qhは，それぞれ運動方程式と連続式に対する重み関数の近似式，p̄ はペナルティー定数，ūh はペナルティーϵ+ u · ∇ϕ = ∇ · (δ(∇ϕ) − Fa ) ,δ∇ϕ
Fa = ϕ(1 − ϕ)
|∇ϕ|

(6)

ここで，ϵ は移動度，δ は界面領域の代表幅である．
また，ϕ は Phase-ﬁeld 変数を表し，気体であれば 0.0，

液体であれば 1.0，自由表面上であれば 0.5 をとるものとする．そして，各要素における流体の密度 ρ と粘性係数 µ は，液体（水）と気体（空気）の密度 ρl , ρgと粘性係数 µl , µg ，Phase-ﬁeld 関数 ϕ を用いて次式のように求められる．

ρ = ρl ϕ + ρg (1 − ϕ)µ = µl ϕ + µg (1 − ϕ)

(7)
(8)

式 (6) に対して，SUPG 法に基づく安定化有限要素法を適用し，解析を行っている．詳しい定式化については，文献 [6] を参照されたい．

3

構造の数値解析手法

3.1

個別要素法を用いた剛体要素

球要素を用いる個別要素法 [7] では，球要素を剛境界面での速度である．安定化項である第 4,5 項は結することで複雑な剛体モデルを表現することが可

要素ごとに定義される不連続量であるため，要素ご能である．このとき，剛体を構成する球要素がそれとの領域 ΩeP での積分の総和となる．ここで，e はぞれ個別に接触判定を行っており，球要素ごとに計流の卓越に対して安定化を施す SUPG 項，および圧する力に換算することで剛体の運動を表現している．

力振動を回避するための PSPG 項であり，第５項はすなわち，剛体要素 g の重心に関する次式の運動方

自由表面の数値不安定性を回避するための衝撃捕捉

程式を解く．

要素番号，nel は要素数である．また，第 4 項は移（Shock-Capturing) 項 [4]，ペナルティ項は構造壁面
上での Dirichet 型境界条件を考慮するために付加したものである．また，τsupg ，τpspg ，τcont は，すべて安定化パラメータである．

算された接触点と接触力の情報を剛体の重心に作用d(mg vg )
= Fg
(9)
dt
d(Ig ωg )
= Mg
(10)
dt
ここで，mg は剛体の質量，vg は剛体の速度ベクトル，Ig は慣性モーメントテンソル，ωg は剛体の角速度ベクトルである．重心に作用する力 Fg 及びモーメ

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

図 1: Cohesive model による DEM 剛体モデルの連結

図 2: Cohesive model のイメージ

ント Mg は剛体を構成する球要素 i に対して接触しを用いて連結することで破壊挙動を表現する．本研ているすべての要素に関する接触力の総和として以下のように書き下せる．

Fg =であり，次式で表現される．

i
∑

Mg =

究で用いた Cohesive model は，多項式を用いたもの(Fn + F s )

(11)

i
∑

(N × F s )(12)

このとき，xi を球要素 i の重心座標，xg を剛体の重心座標として N = (xi + ri n) − xg である．なお，ri は球要素 i の半径，n は球要素 i の重心から接触してい( )[( ) ( )2 ]


27
δ
δ
δ


1
−
2
,
σ
 4 c δf
δf + δfσcoh = 



 0,

δ ≤ δf ,(16)

δ > δf ,ここで，σcoh は結合力，σc は結合力の上限値，δ は開口量，δ f は完全な破壊に至る開口量を意味する．

る球要素の重心へと向かう法線方向の単位ベクトル

このモデルのイメージ図を 図-2 に示す．なお，本研である．また，本研究では, 球要素は剛体モデルの表

究では，負荷時による破壊に着目しているおり，簡面にのみ配置し，計算負荷を軽減している．剛体の便のため，除荷時も式 (16) を用いるモデル化を行っ

姿勢管理には，クォータニオン（四元数）を導入し

ている．

ている．
式 (9) から，時間ステップ n と時間ステップ n + 1

における剛体の速度ベクトルを vng ，vn+1
とすると，
g

vn+1
= vng + (Fg /mg )∆tg(13)となる．ここで，∆t は微小時間増分量である．また，
式 (10) から，時間ステップ n と時間ステップ n + 1 にとすると，おける剛体の角速度ベクトルを ωng ，ωn+1
g

4

数値解析例
数値解析例として，建物に津波が襲来した場合の構造流体連成解析結果を示す．DEM の並列計算には，DEM の解析で用いている探索用のバックセルを

図-3 に示すように領域分割をし，MPI による並列計算を行っている．解析結果として，建物が津波によ(14)

り崩壊する様子を 図-4 に示す．この図に示すように，となる．また，剛体重心の位置座標 xg の更新は，剛津波波力により，壁面が壊されていく様子が見てと体の速度ベクトルを用いて次式により更新する．

トルネックになっていた DEM の計算部分が高速化

−1 n
n
ωn+1
= ωng + (I−1gg Mg − Ig ωg × Ig ωg )∆t

xn+1= xng + vng ∆tg(15)

れる．次に，並列化効率を 図-5 に示す．いままでボ

され，変更前に比べて，約 3.2 倍の高速になってい

ることがわかる．

3.2 Cohesive model による破壊の表現
本研究では，図-1 に示すように，剛体モデルの表面を潜在的な不連続面として，この面を Cohesive model5

おわりに
本研究では，構造物の破壊過程も考慮可能な構造

流体連成解析手法を構築し，その並列化性能向上の京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 43

図 3: バックセルを用いた領域分割法

図 4: 解析結果

ため，DEM を用いた剛体要素の並列化について検討を行った．DEM の並列化により高速計算が可能となり，より実現象に近い解析が可能となった．しかし，今回の解析では，バックセル領域を分割後，領域分割の変更は行わないため，剛体要素が動くことにより，各ノードでの計算負荷のばらつきが多くなり，計算効率の低下が見られているので，今後は，動的負荷分散を考慮した DEM の並列計算が課題となる．

参考文献[1] Terada, K., Asai, M. and Yamagishi, M.: Finite Covermethod for linear and nonlinear analyses of heterogeneous solids, International Journal for Numerical Methods in Engineering, Vol.58, Issue 9, pp.1321-1346,
2003.
[2] Brooks, A.N., Hughes, T.J.R.:
streamlineupwind/Petrov-Galerkin formulations for convection
dominated ﬂows with particular emphasis on theincompressible Navier-Stokes equations, Computer
Methods in Applied Mechanics and Engineering, 32,
pp.199-259, 1982.
[3] Tezduyar, T.E.: Stabilized ﬁnite element formulations
for incompressible ﬂow computations : Advanced in
Applied Mechanics, 28, pp.1-44, 1991.

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

図 5: 並列化効率[4] Aliabadi,S., and Tezduyar, T.E.: Stabilized-ﬁniteelement/interface-caputuring technique for parallel
computation of unsteady ﬂows with interfaces, Compute Methods in Applied Mechanics and Engineering,
190, pp.243-261, 2000.
[5] N. Tanaka, J. Matsumoto and S. Matsumoto: Phaseﬁeld model-based simulation of motions of a two-phaseﬂuid on solid surface, Journal of Computational Scienceand Technology, Vol.7 No.2, pp.322-337, 2013.
[6] 高瀬慎介，森口周二，寺田賢二郎，小山直輝，金子賢治，車谷麻緒，加藤準治，京谷孝史：構造物の破壊を考慮した有限被覆法に基づく構造流体連成解析，
土木学会論文集 A2(応用力学），Vol.71,No.2（応用力学論文集 Vol.18), I 203-I 212, 2015.
[7] Cundall, P. A. : A Computer model for simulating progressive, Large Scale movement in blocky rock system,
Proceedings of ISRM Symposium, pp.11-18, 1971.

産業集積の空間的同期と都市規模・空間パターンに関する研究
森知也
京都大学経済研究所

1 緒言
log (人口規模)｢集積｣は今日の経済立地の代表的形態であり、先進国
及び途上国における先進地域では、人口の 80％以上が都市に集中する。この事実は、国内に限らず、国際地域における経済地理を理解する上で、経済集積原理の解明が鍵であることを意味している。経済集積の特徴的な性質は、その規模の多様性である。例えば、2010年の日本では、人口の 90% 以上が住む 200 余りの都市の人口は、平均が 55 万人であるのに対し、3500 万人の東京から 2 万人の倶知安まで幅広い。
都市規模の多様性は、都市の空間パターンとともに、明確な秩序によって特徴づけられることが知られている。具体的には、都市群は、大都市の周囲に小都市群が集積する｢空間結束性｣を示し、それが再帰的に生起することによる空間的フラクタル構造を持ち、更に、フラクタルにおける階層間で、都市規模分布が共通の冪乗則を示すことが、日米独のデータを用いて示れている [1]。
図 1,2は、アメリカとドイツで同定された地域経済圏における都市規模分布を、それぞれ示している。図2には、ドイツの都市規模分布に加えてシカゴ経済圏の都市規模分布も示してあり、それは、ドイツ全域を含むケルン経済圏の都市規模分布とほぼ一致していることが分かる。この事実は、地域経済圏における都
市規模分布の共通冪乗則が、国を超えて成立している可能性を示唆している。
本研究では、規模の経済の程度が異なる多数の産業を含む仮想経済において、多様な規模の経済集積が形成され、上述の、都市規模･空間分布における、フラク
タル構造と共通冪乗則が発現することが、安定均衡の一般的な性質であることを、モンテカルロ・シミュレーションによって示すことを目的としている。現実経
済と同定度の都市数、及び、都市規模の多様性を実現するためには、産業数 × 都市数に相当する数十万次元の高次非線形最適化問題を解く必要があり、本共同研
究では、そのための解析プログラムの高度化･並列化を行った。既に、昨年度の高度化共同研究において、1次元円環立地空間の下で産業毎の計算を並列化し、シ
ステム D を用いて 7-8 倍程度の高速化を実現していNew York
Los Angeles
ChicagoDallas110log (人口規模順位)

Philadelphia
BostonDetroit
San Francisco
Phoenix
Seattle100log (人口規模)

図1. アメリカの地域経済圏における都市規模分布(2007 年)

110

log (人口規模順位)

100

図2. ドイツの地域経済圏における都市規模分布 (2008年)

る。本年度は、これに加えて、都市毎の計算について
並列化を行い、更なる高速化を試みるとともに、立地
空間を 2 次元周期境界空間に拡張した。

2 モデルと数値解析における課題
分析には、多産業・多地域に拡張した Pﬂüger [2] モデルを用いた。異なる産業は異なる種類の消費財を生産し、各種消費財はさらに多数の差別化された財バラ
エティで構成される。各産業の集積は、財の多様性を

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

嗜好する消費者、生産における規模の経済、及び、輸送費用の存在から生ずるフィードバック効果によって形成され、集積の規模と空間的周期は、フィードバ
ック効果の大きさにより異なる。産業の集積はその生産財の消費者/労働者の集積を伴って都市を形成する。
大都市は、共通の消費者を介して、複数の産業が共集積することにより形成される。
仮想経済では、立地空間上に一定間隔で配置した地点群からなる立地空間を仮定し、所与人口の労働者 (= 消費者) を各産業・地点にランダムに配置した
初期状態から、労働者の産業間・地点間分布を所与とした “短期均衡” の達成と、労働者の効用水準改善を目的とした移住を繰り返し、最終的に各労働者が移住の動機を持たない “長期均衡” に到達する。各短期均衡の導出は、産業数 × 地点数の未知数に対して、産業毎に (地点数 × 地点数) の行列演算を繰り返す
do-loop を複数含む。今回のプログラム高度化は、これら do-loop の OpenMP と MPI を用いたハイブリッド並列化・周期境界型立地空間から生ずる対称性を活かした演算の効率化を中心に行われた。

計算時間間の場合と同様に、並列化による高速化は、基本的には産業単位の並列化のみが有効であった。

3.2

を含むプログラム高度化による計算効率の改善を示している。OpenMP を用いたごく単純な並列化を含むオリジナル･プログラムは、システム D を用いて、各シミュレーションを 1 プロセス･28 スレッドの下で実行した。図は、ノード/コア数と実行時間の関係を示している。いずれも、(地点数, 産業数)= (512, 256)の下で、5 万回の短期均衡計算の所要時間で評価している。オリジナル･プログラムに対して、高度化プロ
グラムでは、7 ～ 20 倍の速度改善があった。一方で、並列化による高速化は概ね産業毎の並列化によって実現されており、都市毎での並列化による計算速度の改善は殆ど得られなかった。
5000

計算時間

オリジナル
産業単位の並列化
産業及び都市単位の並列化

2000
1000
0

1×28 2×14 2×28 4×14 4×28 8×14 8×28 16×14 16×28 32×14
ノード数 × コア数

図3. 1 次元立地空間の下での高度化･並列化による計算速度の改善 (短期均衡 5 万回)

2 次元立地空間の下では、オリジナル･プログラムに対して 3 ～ 6 倍の速度改善があった。1 次元立地空246 − Academic Center for Computing and Media Studies, Kyoto University

1×28 2×14 2×28 4×14 4×28 8×14 8×28 16×14 16×28 32×14
ノード数 × コア数

シミュレーション結果
実際の経済分析への応用では、アメリカの貿易デー
タから推計された約 14000 品目に関する価格マーク
アップ率を用いて、産業構造を特定している。具体的には、この分布から、所与数の品目をランダムに抽出し、抽出した個々の品目の価格マークアップ率を個々の産業の生産品目の価格マークアップ率とみなす。価格マークアップ率は規模の経済の指標と解釈でき、この分布の形状が、均衡における都市規模分布の形状を規定する。
現状、2 次元立地空間の下での均衡計算は十分に高速化されていため、本年度は 1 次元立地空間の下でのブートストラップ標本の蓄積に特化した。具体的には、仮想経済においてアメリカと同程度の都市の数･規模の多様性を実現できる (産業数, 地点数) = (256,1024) の場合を基本ケースとした。
図 5は、安定均衡における典型的な人口の空間分布を示している。輸送費水準によって、人口集積の程度は異なるが、人口の空間分布がこのようなフラクタル構造を持つ性質は共通している。この例は、財の輸送費用が比較的低い場合で、図が示すように、単一の都市に圧倒的な人口集中が起こる一極型構造が実現することが知られている。なお、輸送費が高い場合には、ほぼ同様規模の「最大都市」が複数が存在する「多極型となる。
現時点で、一極型の場合について安定均衡のブート
ストラップ標本を 823 生成し終えている (各パラメタ
ー設定下で 1000 標本の生成を行う)。計算にはシステ
ム D を用いており、4 プロセス･28 スレッドで計算した。図 6は、均衡導出に要した 1 標本当りの計算時間の分布の分布を示しており、平均は 9.36 時間であった。なお、新システム A を用いた場合、2 プロセス×68 スレッドの下で平均 28.97 時間と、システム Dに比べて平均で 3 倍の計算時間を要している。
図 7は、上述の米独の実証分析と同様の方法で同定した、仮想経済における各地域経済圏の都市規模分布

3.1 計算性能
図 4は、1 次元立地空間の下で、産業単位での並列化

3000

オリジナル
産業単位の並列化
産業及び都市単位の並列化

図4. 2 次元立地空間の下での高度化･並列化による計算速度の改善 (短期均衡 5 万回)

3 結果

4000

700
600
500
400
300
200
100
0

人口規模

6e+10

4 結語
本共同研究では、現実の空間経済において観察される、都市規模分布の空間的フラクタル構造と共通べき乗則の創発を再現する数値解析プログラムの高度化･並列化を行った。特に、計算の並列化と効率化によって、1 次元立地空間下では大規模なモンテカルロシミュレーションが可能になった。今後は、2 次元立地空間における計算性能の改善を模索するとともに、1 次元立地空間の下で異なる輸送費水準下でシミュレー
ションを進め、現実に観察される都市規模･空間分布の秩序形成を理論的に再現し、そのメカニズムの解明
を試みる。

4e+10

2e+10

0
0

200

400
600
立地点

800

1000

図5. 仮想経済における都市の位置と人口規模.1

平均9.36時間

References

シェア

.08

[1] Tomoya Mori, Tony E. Smith, and Wen-Tai
Hsu. Cities and space: Common power laws
and fractal structures. Unpublished manuscript, Institute of Economic Research, Kyoto
University, 2016.
[2] Michael Pﬂüger. A simple, analytically solvable, chamberlinian agglomeration model. Regional Science and Urban Economics, 34(5):565–
573, September 2004.

.06
.04
.02
0

0

20

計算時間 (時間)4060

図6. 計算時間を示している。図から明らかなように、都市規模分布
は地域経済圏間で共通の冪乗則に従っている。但し、地域経済圏間の共通冪乗則の成立については、今後、各輸送費水準の下で生成した個々のブートストラッ
プ標本について、[1] による実証研究と同様の方法で統計的に検証し、安定均衡の一般的性質として、厳密に結論付ける予定である。

1e+11

610
985
221
449
771
831
331

log(人口規模)

1e+10
1e+9
1e+8
1e+7
1e+6
1e+5
1e+4
1

10
100
log(人口規模順位)

1000

図7. 仮想経済の地域経済圏における都市規模分布3
京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 47

「メッシュ状透過性平行平板間乱流の
直接数値シミュレーション」に関するプログラム高度化
横嶋 哲 ∗∗

静岡大学 学術院工学領域 数理システム工学系列

はじめに

1

∇ · u =0.

「メッシュ状透過性平行平板間乱流の直接数値シ
ミュレーション」は，京都大学 学術情報メディアセ
ンターのプログラム高度化支援事業（2015 年度）に採択され，続いてスーパーコンピュータ共同研究制(2)

ここで u は流速ベクトルを，τ は流体の応力テンソル
(
)
T
(3)
τ = −p′ I + µf ∇u + (∇u)
を，p′ は流体圧力を示す．また，fIB は固体境界表面度（大規模計算支援枠（2015 年度・後期））として上で粘着・不透過条件を課すために IB 法で導入されも採択された．後者（大規模計算）は前者（プログる単位質量当たりの体積力（仮想外力）を表す．fIB

ラム高度化）の作業環境としても利用された．本報の決定方法は 2.2 節で示す．なお，式 (3) を式 (4) にでは前者の概要を報告する．後者の概要は文献 [1] を代入して整理すれば，次式を得る：

参照されたい．

∂u+ (u · ∇)u = − ∇p + νf ∇2 u + fIB .
∂t

ここで p ≡ p′ /ρf ，νf = µf /ρf は流体の動粘性係数を

高度化前のプログラムの概要

2

(4)

表す．

プログラム高度化対象は，非圧縮なニュートン流

IB 法では式 (4) および (2) を流体領域だけでなく体の 3 次元非定常解析（直接数値シミュレーション）

固体を含む計算領域全体で解く．固体境界が流れに及ぼす影響を IB 法で扱う利点は，流体解析用の計コードであり，特徴としては

算格子点を固体境界形状とは独立に設定できる点に• 有限差分法による離散化ある．

• 3 方向周期境界条件の適用

2.2• 埋め込み境界法（IB 法）[2, 3] の導入• 高速フーリエ変換を利用した，線形連立方程式の効率解法の導入

ように表す：(5)

ここで rhsn+1/2 は式 (4) の移流項および応力項を時支配方程式

刻 tn と tn+1 = tn + ∆t の間で評価したものを表す．

熱的に一様な密度 ρf ，粘性係数 µf の流体の支配
方程式は以下のように表現できる：∂u1+ (u · ∇)u = ∇ · τ + fIB ,∂tρf

流体の運動量保存則 (4) を時間離散化して以下のun+1 − unn+1/2= rhsn+1/2 + fIB.
∆tが挙げられる．

2.1

埋め込み境界法による境界条件の適用(1)48 − Academic Center for Computing and Media Studies, Kyoto University

u の定義位置のどこかで u = u(d) の拘束条件を課す
場合には，外力 fIB を
n+1/2

fIB

=

u(d) − un− rhsn+1/2∆t(6)と与えればよい．ただし拘束条件を課す位置は必ずしも式 (4) の評価位置とは一致しない．

のように ri = (xi − Xi )/∆E に関する 1 次元関数の掛け合わせで表される．関数 ψ にはいくつかの選択

一般に IB 法では固体境界上に点（以下，IB 点と肢があるが，ここでは幅 3∆E の次式 [4]呼称）を配置し，IB 点の集合で固体表面形状を表現


1
3|r|


6 (5 −√−3(1 − |r|)2 )+ 1),(− √ψ(|r|) =1

−3|r|2 + 1 ,

3 1+



0.

する．Uhlmann[2] の IB 法では，速度境界条件を課すための外力をこの IB 点上で評価する：
n+1/2

FIB

=

U(d) − Un− RHSn+1/2 .
∆t(7)

以下では IB 点上で評価される物理量を大文字で示
すこととする．U(d) は境界条件で，境界が静止していれば U

RHS

n+1/2

(d)

n

= 0 となる．式 (7) は Ũ = U +∆t の表現を用いれば
n+1/2

FIB

=

1
3
2 ≤ |r| ≤ 2|r| ≤ 12otherwise

(13)

(i)

を採用した．すなわち IB 点（Xj ）と流速定義点（xl ）で情報のやり取りが生じるのは，互いの距離が全ての方向に 1.5∆E 未満のときに限られる．式 (10)および (11) は，これを用いてオイラー系とラグラン

U(d) − Ũ∆t(8)と書ける．この中間流速 Ũ は流速 u の定義位置ではũ − un= rhsn+1/2∆t(9)のように求められる．

ジュ系で情報転移を行っても流体に働く力やトルクの総量が変わらないという特長を有する [5]．

2.3支配方程式の離散化

支配方程式の離散化の概略は以下のようである：

以上より，IB 点上での中間流速 Ũ を ũ から，お
n+1/2

n+1/2

を FIB

式 (4), (2) の時間積分については，粘性項を 2 次精度

からそクランク・ニコルソン法で陰的に，移流項を 3 次精

れぞれ評価できれば，式 (9) および (8) を用いるこ

度ルンゲ・クッタ法（RK3）[6] で陽的に扱った．未とで u = u(d) の拘束条件を任意の位置で課すことが知量である流速ベクトルと圧力はスタガードに空間できる．Uhlmannn の IB 法ではラグランジュな IB配置し，完全保存型の 2 次精度差分法 [7] によって空よび流速定義位置での外力 fIB(i)

点（Xj ）とオイラーな流速定義点（xl ）での情報間離散化を行った．そのようにして得られた離散方

交換を離散デルタ関数を用いて以下のように行う：
) ∑
(
)
(
(i)
(i)
=
U Xj
u(xm )δh xm − Xj ∆3E ,(10)

程式に部分段階法を適用して解いた．以下では上述

mfIB (xl ) =

∑∑
p

q

(

) (
)
δh xl − X(p)∆Vq(p) .
FIB X(p)qq(11)

ここで ∆E は流体解析用のオイラー系のデカルト格(p)子間隔を，∆Vq は固体 p の j 番目の IB 点が受け持つ流体領域の体積を表す．Kempe and Fröhlich[3]は，Uhlmann が提案した上述の fIB の決定法（フォーシング・スキーム）がフォーシング誤差を伴うことを指摘し，その誤差の低減法を提案した．本研究での方法で支配方程式を時間方向にのみ離散化し，プ
ログラム内で実際に扱われている離散式を導出する．
式 (4), (2) を時間方向に離散化すると次式を得る：()
uk − uk−1
2βk=−∇pk + βk νf ∇2 uk + uk−1∆tρfk
− γk Nk−1 − ηk Nk−2 ,+ fIB∇ · uk =0.

(14)
(15)

ここで N は移流項 (u · ∇)u を表す．上付きの添え字

は RK3 のサブステップ（k = 1, 2, 3）を表し，u0 は

古い時刻の，u3 は新しい時刻の流速をそれぞれ表す．

も Kempe and Fröhlich の修正フォーシング・スキー

他の物理量も同様である．β, γ, η は RK3 の係数で，ムを採用するが，その詳細については支配方程式のβ1 = 4/15, γ1 = 8/15, η1 = 0, β2 = 1/15, γ2 =時間離散化と併せて 2.3 節で述べる．

5/12, η2 = −17/60, β3 = 1/6, γ3 = 3/4, η3 =−5/12 のように与えられる．

離散デルタ関数 δh (x − X) はδh (x − X) = ψ(r1 )ψ(r2 )ψ(r3 )/∆3E

式 (14) および (15) に部分段階法を適用し，式 (10),(12)(8), (11) を組み込めば，流体挙動に関する支配方程京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 49

periodicity imposed in x2 as well as in x1 & x3
x2

x1

L3

x3

L2
P

P+∆P

Flow
D
W

S

(b) another view of the system.

L1(a) ﬂow to be studied.
図 1: 検討対象流れの概要．
式は以下のように扱われる：

ソン方程式については，それぞれ上述の差分法に基ũ − uk−1
2βk=−∇pk−1 + 2βk νf ∇2 uk−1∆tρf

づいて得られた離散方程式を波数空間上で解き，得− γk Nk−1 − ηk Nk−2 ,
(
) ∑
(
)
(i)
(i)
Ũ Xj
=
ũ(xm )δh xm − Xj ∆3E ,
m

FkIB

(

(i)
Xj

)

(

(i)

U(d) Xj)((i)− Ũ Xj

られた数値解を物理空間に戻した．上記の操作に必

(16)
(17)

)

,
(18)
∆t
(
)
(
)
∑∑
k
fIB
δh xl − X(p)∆Vq(p) ,(xl ) =
FkIB X(p)
q
q
=

p

q

∗

(
)
u − ũ
k=βk νf ∇2 u∗ − uk−1 + fIB
,
∆t
∗
∇·u
,
∇2 ϕ k =2βk ∆tu k − u∗
2βk
=−
∇ϕk ,
∆t
ρf
pk =pk−1 + ϕk − βk ∆tνf ∇2 ϕk .

要な 3 次元正/逆フーリエ変換には FFTW を利用した．このため，空間離散化は全ての方向に幅 ∆E の
等間隔デカルト格子に基づく．

3

検討対象流れ：メッシュ状透過性
平行平板間乱流(19)

2 節でその概要が示されたコードがターゲットとす(20)

る流れの概略を図 1 に示す．図 1(a) が本研究で実際に対象とした系を表す．境界接線方向に加えて，境

(21)
(22)
(23)

k
以上は強制外力 fIB
を Uhlmann[2] に倣って決定した場合であるが，Kempe and Fröhlich[3] は Uhlmannの方法では IB 点上で U(d) と U∗ にズレが生じることを指摘し，修正フォーシング・スキームを提案した．

Kempe and Fröhlich の方法では，式 (20) で求まった u∗ を式 (17) を用いて IB 点に移し，次にそこで得()
られる F′IB = U(d) − U∗ /∆t を式 (19) によって u′の定義位置に戻し，最後に u∗ ← u∗ + ∆t · fIB
によっ

k
k
′
← fIB
+ fIB
とする．
て u∗ の値を修正し，併せて fIB

彼らはこの操作を数回程度反復することでフォーシング誤差が低減されることを示しており，本研究では 3 回反復後に式 (21) 以後のプロセスに進んだ．

式 (20) のヘルムホルツ方程式および式 (21) のポア

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

界法線（x2 ）方向にも距離 L2 の周期性が課されており，2 枚のメッシュ状平板は同一境界を表すことに留意されたい．すなわち，図 1(a) は，図 1(b) のように無数のメッシュ状境界が等間隔に並んだ系の幾何学的最小ユニットにあたる．
次節で概要を示すプログラム高度化作業のための
具体的な検討対象として，平板間距離 L2 の半分 δを基準として，主流方向に L1 = 3.84δ ，横断方向に

L3 = 1.92δ の領域を確保し，メッシュ厚 D を δ/150，
メッシュ間空隙 S をゼロとした．このときメッシュ状透過性平行平板間乱流は通常の不透過な平行平板間乱流に帰着する．主流 x1 方向に摩擦レイノルズ数Reτ ≡ Uτ δ/νf = 150 となる一定の圧力差を課すことで流れを駆動した．この直方体領域を全ての方向に ∆ = 1.0(νf /Uτ )，すなわち壁座標で 1 の等間隔格子で離散化した．すなわち，x1 , x2 , x3 の各方向の計算格子点数はそれぞれ 576, 300, 288 となる．

表 1: 高度化による計算速度向上．例えば “p32t2” はプロセス数が 32，プロセス当たりのスレッド数が 2 を示
す．P3DFFT での 2 次元領域分割は，p64t1 の場合に主流 x1 方向と横断 x3 方向に 8 × 8 分割，p32t2 の場合に 8 × 4 分割とした．FFTW では横断方向の 1 次元分割である．
オリジナルコード
並列度
経過時間
速度向上比

4

911870 sec
1.0

FFTW（1 次元領域分割）

P3DFFT（2 次元領域分割）

p64t1 4node

p32t2 4node

p64t1 4node

p32t2 4node

28999.55 sec
31.4

22116.95 sec
41.2

21685.39 sec
42.0

18683.11 sec
48.8

プログラム高度化の概要と性能数（本検証では 576）を超えるような大規模並列環

評価も明瞭に認められる．この結果，オリジナルコード

プログラム高度化として，2 節でその概要が示さ境下だけではなく，高々4 ノードの本検証においてと比較して約 49 倍の高速化が確認できた．

れたコードに対して，MPI によるプロセス並列とOpenMP によるスレッド並列を組み合わせたハイブ

謝辞：プログラム高度化の実施にあたって，クレリッド並列化を施した．高度化前のコードでは，式イ・ジャパン・インクの鈴木 幸朗 様，および高橋(20) のヘルムホルツ方程式および式 (21) のポアソン方程式の数値解法に伴う 3 次元高速フーリエ変換に由樹 様から多大な協力を得ました．記して謝意を表します．

60% 近い処理時間を要することが確認されており，
この点への対処が最も重要である．他方で固体境界(i)形状を表現する点集合（IB 点）の座標 Xj および，
そこでの流速 U や外力 FIB については，全てのランクで同じ値を保持することとした．
上述の高速フーリエ変換の高速化として，MPI 並列に領域分割を導入する場合，現在利用している FFT
ライブラリである FFTW では，1 次元分割のみ可能
である．これでは，当該 1 次元方向の計算格子数よりも MPI で使用するプロセス数が大きいような大規模並列化が行えないことになる．よって本高度化では，2 次元分割が可能な FFT ライブラリを調査した．

Fortran 対応，r2c/c2r 対応，MPI と OpenMP のハ
イブリッド並列対応といった基準を全て満たすものの中から，開発歴，使用実績，ライセンス形態やメン
テナンス状況を考慮して，今回は P3DFFT を新たに導入することとした．得られたハイブリッドコードの性能評価を京都大学のスーパーコンピュータシス
テム B の 4 ノード（64 コア）上で行った．主要な結果を表 1 に示す．FFTW では最も計算格子点数が多い主流 x1 方向について 1 次元領域分割を，P3DFFT
では x1 と横断 x3 方向についての 2 次元領域分割が施された．

P3DFFT の導入の効果は，当初に問題が懸念されたような，分割数が 1 次元領域分割方向のメッシュ

参考文献[1] 横嶋 哲: 「メッシュ状透過性平行平板間乱流の直接
数値シミュレーション」に関する大規模計算， 京都大学 学術情報メディアセンター 全国共同利用版 広報, Vol.14 No.2, 2015 （印刷中）.
[2] Uhlmann, M.: An immersed boundary methodwith direct forcing for the simulation of particulate flows, J. Comput. Phys., Vol.209, pp.448–476,
2005.
[3] Kempe, T. and Fröhlich, J.: An improved immersed boundary method with direct forcing forthe simulation of particle laden flows, J. Comput.
Phys., Vol.231, pp.3663–3684, 2012.
[4] Roma, A.M., Peskin, C.S. and Berger, M.J.:
An adaptive version of the immersed boundary
method, J. Comput. Phys., Vol.153, pp.509–534,
1999.
[5] Peskin, C.: The immersed boundary method, Acta
Numerica, Vol.11, pp.1–39, 2002.
[6] Spalart, P.R., Moser, R.M. and Rogers, M.M.:
Spectral methods for the Navier-Stokes equations
with one infinite and two periodic directions, J.
Comput. Phys., Vol.96, pp.297–324, 1991.
[7] Morinishi, Y., Lund, T.S., Vasilyev, O.V. andMoin, P.: Fully conservative higher order finite difference schemes for incompressible flow, J. Comput. Phys., Vol.143, pp.90–124, 1998.

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

乱流輸送コード GKNET への運動論的電子の実装
今寺 賢志
京都大学大学院 エネルギー科学研究科 エネルギー基礎科学専攻

磁場閉じ込め核融合プラズマでは、電子ダイナミックスが深く関与した捕捉電子モード(TEM)が乱流輸送過程に重要な役割を果たす一方で、電子を運動論的に解く必要があることから、そのシ
ミュレーションを実施するためには、大規模計算に向けて高並列化されたコードの開発は不可欠である。そこで本研究では、我々が開発した大域的乱流輸送コード GKNET に対してハイブリッ
ド並列の実装等を行った。本稿では、それらの実装方法と得られた高速化率、およびそのコードを用いて行った TEM シミュレーションの結果について報告する。

1

研究背景/目的的な内部輸送障壁の形成に成功した[3, 4, 5]。この結果を踏まえ我々は、「乱流構造のポロイダル対称

乱流を抑制することで輸送係数を低下させ、断性を破ることが乱流の安定化に寄与し、結果とし

熱層/粒子遮蔽層として作用する「内部輸送障壁」て輸送障壁の形成をもたらす」との着想に至った。

の形成は、磁場閉じ込め核融合プラズマプラズマを高性能化する上で重要な鍵を握っているが、第一原理に基づいた乱流輸送シミュレーションでは再現されていなかった。これは、大規模な計算機資源を必要とする大域的乱流輸送コードの開発・実行が容易でないこと、さらには径電場に代表される背景分布と乱流の相互作用が物理的に解明されていないことに起因している。
そこで我々は、新たに大域的乱流輸送コード「GKNET (GyroKinetic Numerical Experimen- tal

Tokamak)」を開発した[1, 2]。本コードでは、外部

図 1: 運動量源を導入した場合のイオン温度分布からの熱ソース/シンク、ならびに衝突効果が適切に導入されており、熱ソースによる温度分布の急

そのために重要となってくる物理的要素が、電峻化によって駆動される乱流(熱源駆動型乱流)に子ダイナミックスを第一原理に沿って解いた場合

加え、新古典輸送過程を経て形成される径電場もに得られる運動論的電子効果である。図 2 は、断適切に評価することが可能である。

熱電子を仮定した従来のシミュレーションで明ら

このコードを用いて、熱源駆動型乱流輸送シミかとなった ITG 乱流の安定化/不安定化メカニズュレーションを行った結果、力学的平衡によってムを模式的に示したものである。ITG 乱流では、もたらされる径電場は、イオン温度勾配(ITG)モー

径電場が不安定化に寄与していたため、外部からドと呼ばれるイオン系乱流のポロイダル対称性を十分に大きな運動量ソースの導入を必要としてい回復させる方向に作用し、不安定化に寄与するこた。一方、捕捉電子モード(TEM)呼ばれる電子ダとが明らかとなった。さらにその径電場を変調さ

イナミックスが深く関与したモードが励起されたせるために運動量ソースを導入した結果、図 1 に場合、乱流構造が反転するために、径電場はむし

示すように、順方向のトロイダル回転による強制ろ安定化に作用し、さらに、ITG 乱流と比較して52 − Academic Center for Computing and Media Studies, Kyoto University

より強く励起されると考えられている自発回転[6]を最適化することでその安定化効果が更に増幅されることで、自発的な内部輸送障壁が形成されることが期待できる(図 3)。

2

ハイブリッド並列の導入
本研究課題では、申請時のオリジナルコードに対して主に以下の 3 点に関する作業を行った。
(1) ハイブリッド並列の実装

全ての DO ループにおいて、ハイブリッド並列を実装した。主な利点としては裾値データの通信量の削減と、プロセス・スレッドの性能の良い分割数の組合せを選択できる自由度の増加が挙げられ、特に 105コア程度の大規模計算において有効であると考えられる。
図 2: 先行研究で明らかとなった乱流構造の対称性による ITG 乱流の安定化/不安定化メカニズム(2) スカラーチューニング
ループ順序や配列順序の入れ替えによってキャッシュチューニングを行った。
(3) プロセス並列に伴う領域分割法の変更
オリジナルコードでは5 次元(x, y, z, v, µ)領域分割を用いており、特にµ方向については優先的に分割を行っていたが、その制約を解除して自由に領域分割できるようにコードを変更した。
図 4 は、(A) オリジナルコード、(B) (1)+(2)を行

ったコード、(C) (1)+(2)+(3)を行ったコードの各パートの計算時間を示している。”Vlasov”は 5 次元の双曲線型偏微分方程式である Vlasov 方程式の図 3: 本研究で解明を目指す乱流構造の対称性による

ソルバー部を指しており、”Field”は楕円型偏微分

TEM 乱流の安定化/不安定化メカニズム

方程式である準中性条件のソルバー部 (速度空間の積分に伴う MPI_ALLREDUCE を含む) を指し
しかしながら一方で、電子を運動論的に解いたている。まず(A)と(B)を比較すると、Field 部で計場合、質量比(イオンの質量/電子の質量)の平方根

算時間が大幅に削減されていることがわかる。こに比例して CFL 条件が厳しくなるため、計算量は

れは、準中性条件を解く際の 3 重ループの順序と莫大なものとなり、TEM 不安定性を扱える大域的配列順序が対応しておらず、その結果ストライド

ジャイロ運動論コードは世界的にもまだ極めて希アクセスが発生していたことが原因と考えられる。

少である。

また Vlasov 部でも若干の計算時間の削減が得ら

5

そこで本研究では、10 コアを超える大規模計算れており、特に 128 コアの場合に顕著であったこ

において良並列スケーリングを確保するために、とから、5 次元裾値データの通信時間が短くなっ

従来の MPI 並列をハイブリッド並列(MPI 並列たことに起因しているものと思われる。一方で、

時間で電子を運動論に沿って時間積分できる大域µ方向の積分のために使用している+OpenMP 並列)に拡張することで、現実的な計算的ジャイロ運動論コードを開発し、ITG/TEM 不安

MPI_ALLRUDUCE にかかる通信時間は(B)のケー
スの方が増加しており、128 コアではその比率が定性と温度分布形成の因果関係を明らかにするこ

26%に達していることから、大規模計算においてとを目的とした。

非常に強いボトルネックになることが推測される。
そこで、µ方向の領域分割を避けるために(3)の改善を更に行った。具体的にはμ方向の領域分割

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

図 4: (A)オリジナルコード、(B)Hybrid 並列の実装とスカラーチューニングを行ったコード、(C)Hybrid 並列の実装とスカラーチューニング、領域分割法の変更を行ったコードの Vlasov 部、
Field 部にかかる計算時間とその和の割り当てを最も効率改善に優れていたz方向に回した。その結果、MPI_ALLREDUCE にかかる
時間は大幅に削減され、さらに Vlasov 部についてもキャッシュミスが削減されたことから計算時間が短くなっていることがわかる。

に対応しており、(B)は電子温度勾配が急峻な場合(R 0 /LTi = 2.22, R 0 /LTe = 6.66)で TEM モードに対応している。実周波数が正(ITG)から負(TEM)に変化することで、モードの回転角(ブロッホ角)が反転することが確認され、当初想定していた現象

実際の計算時間を(A)と(C)を比較した場合、16と合致していることがわかった。また、その回転コアのケースで 2.18 倍、128 コアのケースで 1.58角の反転が自発回転と呼ばれる乱流に起因したプ倍の高速化に成功し、当初予定していた通りの成

ラズマ回転も反転させることが確認され、その制果が得ることができたと言える。

御によって内部輸送障壁の自発形成が促進される

3

ITG/TEM シミュレーション
2 節で説明したハイブリッド並列されたコードを用いて、線形 ITG/TEM シミュレーションを行ことが期待できる。

4

総括と今後の課題
本稿では、大域的乱流輸送コードに対してハイっ た。図 5 (A) はイ オン温 度が急 峻な場合

ブリッド並列とスカラーチューニング、領域分割

LTi :イオン温度勾配の特性長、LTe :電子温度勾配の倍、128 コアのケースで 1.58 倍の高速化に成功し(𝑅𝑅𝑅𝑅0 /𝐿𝐿𝐿𝐿 𝑇𝑇𝑇𝑇𝑖𝑖𝑖𝑖 = 6.66, 𝑅𝑅𝑅𝑅0 /𝐿𝐿𝐿𝐿 𝑇𝑇𝑇𝑇𝑒𝑒𝑒𝑒 = 2.22, R 0 :装置の大半径、特性長)に得られたモード構造であり ITG モード

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

法の変更を行うことで、16 コアのケースで 2.18
た。また、そのコードを用いて、線形 ITG/TEM

scale transport with internal transport barrier
formation in magnetic confinement fusion

plasmas”, 18th International Congress on Plasma[5]Physics, Taiwan (2016).

Kenji Imadera, Jiquan Li and Yasuaki Kishimoto,“ITB formation in flux-driven gyrokinetic ITG
turbulence”, Proceedings of 26th Fusion Energy

図 5: (A)イオン温度勾配が急峻な場合と(B)電子温
度勾配が急峻な場合におけるバルーニング構造

シミュレーションを行った結果、ITG/TEM モード[6]

Conference (2016).

Yasuhiro Idomura, “A new hybrid kinetic electron

model for full-f gyrokinetic simulations”, Journalof Computational Physics, 313, 511 (2016).

の再現に成功し、モードの回転角が反転することが確認された。この結果から、TEM 乱流に対しては径電場の安定化メカニズムが増幅される可能性が示唆され、内部輸送障壁の自発的な形成の再現が期待される。
今後の課題としては、熱源駆動型 ITG/TEM シ
ミュレーションを行い、上述の点について検証することが挙げられ、そのために(1)半陰解法の実装
による CFL 条件の緩和、(2)ハイブリッド電子モデ

ル[6]による CFL 条件の緩和も併せて行う予定である。
Reference[1]

Kenji Imadera, Yasuaki Kishimoto, Kevin

Obrejan, Takuya Kobiki and Jiquan Li, “Global
profile relaxation and coupled with E × B

staircase in toroidal flux-driven ITG turbulence”,

Proceedings of 25th Fusion Energy Conference[2](2014).

Kevin Obrejan, Kenji Imadera, Jiquan Li andYasuaki Kishimoto, “Development of a global
toroidal gyrokinetic Vlasov code with new real

space field solver”, Plasma and Fusion Research,[3]

10, 3403042 (2015).

Kenji Imadera, Yasuaki Kishimoto, Kevin

Obrejan, Takuya Kobiki and Jiquan Li,“Non-local

profilerelaxation

andbarrierformation in toroidal flux-driven turbulence”, The

7th IAEA Technical Meeting on “Theory of[4]Plasma Instabilities”, Italy (2015).

Kenji Imadera, Jiquan Li and Yasuaki Kishimoto,“Flux-driven gyrokinetic simulation on multiple

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

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

地震前兆における山岳地表面を伝搬する
電磁プラズマ波の特性解明
藤井
∗

1

富山大学

雅文 ∗

大学院理工学研究部

はじめに
近年、地震の前兆現象として電磁波の異常な長距離伝搬現象が観測されている。この現象は地殻活動に伴って地表に出現する電荷が原因で誘起されることが推察され、さらに電荷の発生メカニズムが半導体における電子や正孔の生成と同様な機構であることが解明されつつある。このように、従来の地震学における測地的な地震予測とは異なり、地震の前兆として生じる電磁気的な現象を捕えようとする研究が進展している。本研究ではこのような電磁気現象のうち、地表面電荷と電磁波、特にラジオ放送波との相互作用を大規模な数値解析により解明することを試みた。
同様な物理現象に金属表面の自由電子による光散乱がある。金属に光が照射されると表面の電子は光

工学領域

た電磁波はそれほど大きな作用をおよぼさないで、一部は地中に浸透し、一部は反射する比較的穏やかなふるまいをする。
このような電磁波と地表面電荷の相互作用を大規
模数値解析によって、特に山岳地帯における実際の複雑な地形を考慮することにより精密に解析する。
このような解析によって得られる知見に基づき電磁波の観測を実施すれば、観測結果の変化によって応力変化などの地殻活動を捕えることが可能になると考えられ、このような電磁波の異常現象を明らかにすることが本研究の目的である。これに関連する論文はすでに米国地球物理学連合の学術誌に掲載され[1]、同学会ホームページにおいて注目の論文として取り上げられた [2]。本稿はこれらの一部を抜粋し加筆したものである。

の電場から受けるクーロン力によって振動しようとする。このとき、自由電子密度が十分大きいと高い振動数の外部電場の動きに追随することができ、このような状況では光は金属内部深くまでは浸透できず反射される。しかし、自由電子密度が小さいと電子は外部電場の動きに追随できず、その物質はもはや金属ではなく誘電体と同様になり光は物質の中に侵入する。すなわち、表面電荷密度の大小により物
質は外部からの光あるいは電磁波に対して金属のよ
うにふるまうか誘電体のようにふるまうかが決まる
こととなる。
同じことが地表面に出現する電荷に対して示すことができる。すなわち、地震の前兆として地表面に電荷が多く存在していれば近傍の電磁波と結合して強く散乱したり回折させる。その反対に地表面に特に多くの電荷が存在していなければ地表面に到達し56 − Academic Center for Computing and Media Studies, Kyoto University2

理論背景
地球上の岩石はプレートと呼ばれる数１０ｋｍの厚さの岩盤を構成しており、これがマントルの流動によってゆっくりと移動する際に特定の箇所に大きな応力が作用する。地中深くの火成岩は高温で溶融した状態で水分子を取り込み、上方に移動し冷えて固まる。この際に水素や酸素原子が岩石中の特定の部位に入り込み、ペルオキシ結合を形成することが示されている。このような水素や酸素を取り込んだ火成岩に応力が作用すると特定のペルオキシ結合が破断され、それに伴ってちょうど半導体の電子と正孔のような電荷が発生する。特に火成岩に応力が作用するとＰ型半導体に類似した状態となり正孔が生成され岩石中を移動して地表に出現する。また、深部の高温の岩石は熱によって励起された伝導電子を

生じ、これらが地中の主に活断層等に沿って伝搬しと電波の散乱を解析するためには１ｋｍ四方程度の地表面に到達する。地表面の異なる場所に到達した地形をモデル化する必要がある。このため数値解析

正孔と電子は再結合するまでは地表を移動しその部は大規模問題となり大型計算機を使用することが必

分が静電気を帯びたような状態になる。このような要となる。

地中の岩石における電子と正孔の生成が理論と実験によって示されている [4]。
地表に出現した電荷が自然界において興味深い現象を引き起こすことについて様々な推察と実際に観察された現象が報告されている。それらの詳細は他

3

数値解析

3.1

手法

の文献にゆずり、ここではその１つであるラジオ放

今回の数値解析には有限差分時間領域法（ＦＤＴ

送波との相互作用を検証する。ラジオ放送波はＡＭ

Ｄ法）を用いた。本手法は電磁波のマクスウェル方

ラジオの長波帯からＦＭラジオの超短波（ＶＨＦ）

程式を差分近似して逐次時間を追って計算するもの

帯が良く知られている。長波帯の電波は周波数が約であり、解析領域を分割することにより効果的な並３０ｋＨｚから３００ｋＨｚであり、上空に存在す

列計算が可能であり、したがって大規模数値解析にる電離層によって反射されやすく地表面との間で反適している。また、プラズマ振動は周波数に依存す射を繰り返して遠方まで伝搬する性質がある。一方、る性質を持ち、この周波数特性を微分方程式の形で超短波帯は周波数が約３０ＭＨｚから３００ＭＨｚ

正確に有限差分法に取り込むことが可能となっていであり比較的周波数が高く、上空の電離層を透過しる。筆者は本手法をナノメートルサイズの金属粒子てしまい通常１００ｋｍ程度より遠方へは届きにくに生じる局在表面プラズモンの解析に適用しその妥い性質がある。

当性は確認済みである [3]。

地震前兆時に地表に電荷が出現するとそれらの電解析対象には図１ (a)-(c) に示すように富山県と長荷は付近を伝搬する電波の電場から力を受け振動す

野県にまたがる北アルプス山岳地帯のうち立山連峰る。これは、例えて言えば湖面の上に風が吹くときの剱岳を選んだ。同峰は標高２９９９ｍであり約２水面が波立つのと同じような現象として理解できる。

５００ｍ以上は森林限界を超え樹木はなく、花崗岩

地表面の電荷が波立つ現象は、表面電荷のプラズマが切り立ち露出している構造である。この外観を図

振動である。光学領域では金属表面に光を照射した

１ (d) に示す。このため山体のモデル化には典型的

ときに自由電子が振動する表面プラズモンとしてよな岩石の電気定数として一様に比誘電率６、導電率

く知られている。電波の領域でも物質表面に電荷が10−3 [S/m] を適用した。解析は図１ (c) に示す直線存在するところに電波が伝搬すると同様の表面プラ

ＡＢ（５ｋｍ）に沿った２次元解析と正方形Ｃ（６０

ズマ波が励振される。

０ｍ×６００ｍ）の範囲における３次元解析を行っ実際、地震や地殻活動による表面電荷がこのようた。このために必要な計算機資源は超並列計算機Ｇな地表面のプラズマ波、あるいは地表面プラズモン

Ｂ８０００による３次元解析について最大６４ノーを生じさせるような表面電荷密度に達することは可ド、１０２４プロセス、記憶容量３．８ＴＢ、計算能である。このことは上述の理論および実験から推時間は１解析当たり約６時間であった。

察されており、超短波帯の電波との相互作用を観測

デジタル標高モデルは国土地理院による５ｍ水平すればこれによって地殻に作用する応力の変化を検分解能標高地図をもとに双３次スプライン補間法を出することが可能となる。すなわち、遠方のＦＭラ用いＦＤＴＤ解析に必要な空間分解能０．２ｍの解ジオ放送波を常時観測しておき、地殻への応力変動

析モデル（図２）を生成した。このモデルの山頂にがあれば受信されるラジオ放送波の強度が変化する

向けて周波数７０ＭＨｚの電波を西から東（ｚ軸のことにより検知できる。このような現象を実際の地正方向）へ伝搬させている。

形を考慮して解析し、正確に電波の伝搬特性を理解することが重要となっている。超短波の波長は数ｍであるが、地形の変化や電波の伝搬距離は数１００
ｋｍにおよぶ。したがって、少なくとも地形の起伏

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

3.2

結果および考察

今回の解析の主な結果である山頂付近の電波と地
表面プラズマ波との相互作用を、図３には西から観察した様子、図４には上方から見下ろした様子を示している。図３ (a) および図４ (a) には山体の誘電率
分布を黒／黄色のコントラストによって示す。地表面に電荷が存在しない場合、すなわち図３ (b)、図４(b) では電波は山頂で弱く回折していることが分かる。電波は山頂の岩石中に侵入してから外部へ出ていく様子が分かる。これに対して図３ (c)、図４ (c)に示す地表に電荷が存在する場合、すなわち地殻の応力が変動していると考えられる状況では、電波は山頂の表面プラズマ波と強く相互作用し、山頂に引

(a)

(b)

き付けられるように反対側まで地表面に沿って伝搬し、さらに地表の凹凸によって散乱される様子が判明した。このような状況では通常電波が到達しない地域に向けて電波が再放射され、高感度の受信機で受信できることが推測できる。
以上のように、電波と地表面電荷の相互作用を大規模な数値解析によって解明することが可能となった。このような結果に基づき様々な地点から出射されるラジオ放送波の受信強度を常時観測することに(c)より、これまで不明であった地殻活動と電波伝搬の関係を推察することが可能になり、ひいては地震につながる地殻の応力変動を検出できる可能性が高まることを期待している。

(d)

参考文献

図 1: 解析領域の位置と地形。(a) 日本地図中、赤の[1] M. Fujii, ”A new mode of radio wave diﬀraction via the terrestrial surface plasmon on mountain range” Radio Science, v.51, pp.1396-1412,
doi:10.1002/2016RS006068, Sep. 2016.

の赤い領域の拡大図。本解析で考慮した北アルプス[2] L. Crane, ”Earthquakes could funnel radio waves
to dark zones in mountains”, American Geophysical Union Eos, 97, doi:10.1029/2016EO059987,
Sep. 2016.
[3] M. Fujii, ”Fundamental correction of Mie’s scattering theory for the analysis of plasmonic resonance of a metal nanosphere” Phys. Rev. A, vol.89,
no.3, 033805, Mar. 2014.
[4] F. Freund, ”Charge generation and propagationin igneous rocks”, J. of Geodynamics, v.33, p.543570, 2002.

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

長方形は日本の北アルプス山岳地帯。(b) 図 (a) 中の領域。赤い四角形は解析領域の位置。(c) 図 (b) の赤い領域（東西に６ｋｍ、南北に４ｋｍ）の劔岳のディジタル標高モデル。頂上をＰで示す。(d) 西側からの剣岳頂上の眺望。図の領域は水平方向幅約３ｋｍ、垂直方向は標高約２２００ｍから３０００ｍの高度差約８００ｍ。視点は剣岳の真西よりやや北方向に寄る。図 (a)-(c) は国土地理院地図、図 (d) は
著者による撮影。

3000

3000

2950

Altitude [m]h [m]; Altitude

2950

2900

2900
2850

2850

300
200
100

x [m]; South to North

0-100-200 -300

-200

-100

0

100

200

z [m]; West to East(a)

図 2: 剣岳頂上の３次元解析モデル。水平方向ｘ軸（南から北）ｚ軸（西から東）、山頂位置がｘｚの原点。標高ｈはｙ軸（下から上）方向。

(a)

(b)

(b)

(c)
図 4: (a) 劔岳３次元モデルにおける頂上付近の比誘(c)
図 3: (a) 劔岳３次元モデルにおける頂上付近の比誘電率分布。黒の領域は空気 ϵr1 = 1、黄色の領域は山体 ϵ∞ = 6 であり、標高 2940m における水平面上の電率分布。黒の領域は空気 ϵr1 = 1、黄色の領域は山

電界分布. (b) および (c) ３次元解析結果。西から電体 ϵ∞ = 6 であり、x = 0 における垂直面上の電界波が入射。(b) 電界分布 |E|、表面電荷が存在しない分布。 (b) および (c) ３次元解析結果。西から電波場合。 (c) 同、表面電荷が存在する場合、ただし実が入射。(b) 電界分布 |E|、表面電荷が存在しない場

効プラズマ周波数 fp′ = 408 MHz。

合。 (c) 同、表面電荷が存在する場合、ただし実効プラズマ周波数 fp′ = 408 MHz。

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

「メッシュ状透過性平行平板間乱流の
直接数値シミュレーション」に関する大規模計算
横嶋 哲 ∗∗

1

静岡大学 学術院工学領域 数理システム工学系列

はじめに「メッシュ状透過性平行平板間乱流の直接数値シ

ミュレーション」は，京都大学 学術情報メディアセ

x2
x1
x3
P

P+∆P

ンターのプログラム高度化支援事業（2015 年度）にD採択され，続いてスーパーコンピュータ共同研究制

W

度（大規模計算支援枠（2015 年度・後期））として∆(a) flow to be studied.

も採択された．後者（大規模計算）は前者（プログ
ラム高度化）の作業環境としても利用された．本報では後者の概要を報告する．前者の概要は文献 [1] を参照されたい．

2

P

P+∆P

検討対象流れ：メッシュ状透過性
平行平板間乱流検討対象であるメッシュ状透過性平行平板間乱流（メッシュ・チャネル乱流）の概略を図 1 に示す．図(b) another view of the system.
図 1: 検討対象流れの概要．

1(a) が本研究で実際に対象とした系を表す．全ての乱流（図 2 左）に，厚み D をゼロ，幅 W と間隔 S方向に周期境界条件が課されていることに留意されを同時にゼロに近づけることで，理想的な透過性 [3]たい．すなわち，図 1(a) は，図 1(b) のように無数を備えた完全透過チャネル乱流 [4] に，それぞれ帰着のメッシュ状境界が等間隔に並んだ系の幾何学的最する．すなわちメッシュ・チャネル乱流は，その幾小ユニットにあたる．

何学形状を系統的に変化させることによって，これこの流れの流体力学的な位置付けの詳細について

ら 2 つの理想的な系の結びつきを明らかにし，境界は文献 [2] を参照されたい．以下ではその概要を述べの透過性が壁乱流に及ぼす影響の理解に大きく貢献

る．図 2 の中央の絵は図 1(a) と同じ系を，x2 軸方向することが期待できる．

にメッシュ平板間距離の半分だけずらして描画した

ここではメッシュ平板間距離を x2 軸方向の周期長もので，2 枚のメッシュ状平板は同一境界を表すこさ L2 と一致させ，δ = (1/2)L2 とする（すなわち δとに留意されたい．メッシュ・チャネル乱流ではメッ

は平板間距離の半分）．x1 および x3 方向にそれぞ

シュの幾何学形状をメッシュ厚 D ，メッシュ幅 W ，れ L1 = 3.84δ ，L3 = 1.92δ の領域を確保し，x1 方

メッシュ間隔 S の 3 つのパラメータでコントロール

向に摩擦レイノルズ数 Reτ ≡ Uτ δ/νf = 150 となる

する．この場合，間隔 S をゼロとすることで，メッ
シュ・チャネル乱流は完全不透過な古典的チャネル

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

一定の圧力差を課すことで流れを駆動した．

メッシュ形状について，ここでは境界粗度の影響

Classical channel flow

Mesh-screen type channel flow

Perfectly-permeable channel flow

periodicity

periodicity on u2 & px2
x1

u1=u2=u3= 0

u1=u3= 0

x3

Flow

D
W

S

S→0 & W→0

S=0D

図 2: 古典的チャネル乱流および完全透過チャネル乱流 [4] に対する，メッシュ・チャネル乱流 [2] の位置付け．

S = 0 で古典的チャネル乱流，W → 0 かつ S → 0 の極限で完全透過チャネル流にそれぞれ帰着．
2500

表 1: 検討対象としたメッシュ状境界の幾何学形状．

ϕ

S+

W+

S/δ

W/δ

M25-12

0.25

12

12

0.08

0.08

M25-24
M25-36

0.25
0.25

24
36

24
36

0.16
0.24

0.16
0.24

M44-16
M44-32
M44-48

0.44
0.44
0.44

16
32
48

8
16
24

0.11
0.21
0.32

0.05
0.11
0.16

M56-18
M56-36

0.56
0.56

18
36

6
12

0.12
0.24

0.04
0.08

M56-54

0.56

54

18

0.36

0.12

2000
1500

Reb

Case

1000
500
0

⋄ φ = 0.25
◦ CC
0

20

⊳ φ = 0.44
• PPC

S+

40

⊲ φ = 0.5660

図 3: バルク・レイノルズ数の時間平均値（プロット）と標準偏差（エラーバー）の境界幾何学形状に対すを議論から排除するために D = 0 とし，空隙サイズ

る依存性．ただし標準偏差は 10 倍に拡大して表示．

S と空隙率（平板境界上を空隙が占める割合）ϕ ≡S 2 /(S + W )2 を表 1 のように変化させた．以下では壁変数（Uτ および動粘性係数 ν ）で規格化された物

3

理量を添え字 + を用いて f + のように表す．
メッシュ状平板が流れに及ぼす影響は埋め込み境界

計算結果と考察
図 3 は各ケースから得られた，断面平均流速 Ub に基づくバルク・レイノルズ数 Reb ≡ Ub δ/ν の時間平

法 [5, 6] で表現した．計算格子は全ての方向に ∆ =均値と標準偏差を表す．ただし標準偏差は 10 倍に拡

2 で等間隔に配置した．比較のため，古典的チャネル
乱流および完全透過チャネル乱流もメッシュ・チャネ

大して表示した．これより，平板境界の空隙率 ϕ がル流と同じ計算領域でシミュレートした．ただしこ

すなわち抵抗が増すことがわかる．また，透過性がの 2 ケースのみ，境界法線方向の計算格子サイズを増して流量が大きく低下した場合にも，その揺らぎ

1 粘性長で等間隔とした．用いた支配方程式やその離散化手法などについては文献 [1] を参照されたい．

の規模はほとんど弱まらない．

+

増す，またはメッシュ空隙 S が狭いほど流量が減り，次に瞬時流れ特性のメッシュ形状に対する依存性を図 4 に示す．図 4 には，(i) 速度勾配テンソルの第2 不変量の等値面で抽出された瞬時渦構造，(ii) いく
つかの 2 次元断面（x1 − x2 および x2 − x3 ）上に主

流方向変動流速 u1 − ū1 のカラーコンター（赤/青は京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.15 No.2（ Mar. 2017）− 61

図 4: 瞬時流れ特性のメッシュ形状に対する依存性．個々の瞬時流れスナップショットは (i) 速度勾配テンソルの第 2 不変量の等値面による渦構造の可視化，(ii) 主流速変動成分 u1 − ū1 のカラーコンターを x1 − x2 および

x2 − x3 断面上で表示（赤:瞬時局所的に平均流 ū1 より高速; 青:平均流より低速），(iii) 瞬時流速 u2 のカラーコンターを x2 = 0 平面上で表示（赤:吹き出し (u2 > 0); 青:吸い込み (u2 < 0)），を重ねて示す．
その領域が平均流 ū1 より速い/遅い流体で占められれ構造は，完全透過チャネル乱流で観察されたものていることを意味），および (iii) メッシュ状境界面と多くの共通点を有する．

上（x2 = 0）ではメッシュ空隙での瞬時流速 u2 のカの境界上での
乱れ強度の境界法線方向成分 urms
2

ラーコンター（赤は吹き出し（u2 > 0），青は吸い

値はメッシュ空隙を通じた流体の往来の活性の度合込み（u2 < 0）を意味），を併せて表示した．空隙いを表し，透過性境界の重要な特徴量のひとつと言率 ϕ = 0.25 の場合にはメッシュ空隙サイズ S の値

える．ここではメッシュ状境界から計算格子 1 つ分に依らず，古典的チャネル乱流と同様に主流方向に（∆+ = 2）離れた位置における urms
2 （ただしケース

伸びた管状の渦構造が遍在し，その周囲でスウィー
プやイジェクションを誘起する様子が見られる．空CC および PPC では境界上の値を使用）の境界幾何学形状に対する依存性を図 5 に示す．空隙率 ϕ が高隙率 ϕ の値が 0.44, 0.56 と増すにつれて管状渦は減

いほど，そして空隙サイズ S が小さいほど吹き出し/

少し，メッシュ境界の縁に張り付いたひだのような吸い込み強度が強まる傾向が認められる．図 3 と併構造が多く出現する．この傾向は空隙サイズ S が小せると，境界での吹き出し/吸い込み強度が強いほどさくなるにつれて強められる．管状渦の減少に伴っ抵抗が大きい．これは境界空隙を通して流体の往来て，スウィープやイジェクションが流路横断方向にが可能となり，より広い断面内で流体混合とそれに一様に生じる傾向も認められる．すなわち，空隙率

伴う運動量交換が活性化されたことで，より大きなを高めて空隙サイズを狭めたケースで生じた瞬時流抵抗が生じたと理解できる．

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

1.5

参考文献
[1] 横嶋 哲: 「メッシュ状透過性平行平板間乱流の直接
数値シミュレーション」に関するプログラム高度化，
京都大学 学術情報メディアセンター 全国共同利用版
広報, Vol.14 No.2, 2015 （印刷中）.

urms,+
2,x2 ≃0

1.0[2] 横嶋 哲: 乱流中の透過性境界の抵抗特性， ながれ,第 35 巻 第 4 号, pp.307–312, 2016.

0.5

0⋄ φ = 0.25

⊳ φ = 0.44

◦ CC

• PPC

0

20

S+

40

⊲ φ = 0.5660

図 5: メッシュ境界空隙を過ぎる吹き出し/吸い込み
速度 u2 の乱れ強度の境界幾何学形状に対する依存性．

4

おわりに
本報ではこれまで最も研究が積み重ねられてきた

壁乱流と言える平行平板間乱流（古典的チャネル乱流）の平板に理想的な透過性を付加した系とみなせ[3] Perot, B. & Moin, P.: Shear-free turbulent boundary layers. Part 1.: Physical insights into nearwall turbulence, J. Fluid Mech., Vol.295, 199–227,
1995.
[4] Yokojima, S.: Eﬀect of wall permeability on
wall-bounded turbulent ﬂows, J. Phys. Soc. Jpn.,
Vol.80, 033401, 2011 .
[5] Uhlmann, M.: An immersed boundary methodwith direct forcing for the simulation of particulate ﬂows, J. Comput. Phys., Vol.209, pp.448–476,
2005.
[6] Kempe, T. and Fröhlich, J.: An improved immersed boundary method with direct forcing forthe simulation of particle laden ﬂows, J. Comput.
Phys., Vol.231, pp.3663–3684, 2012.

る完全透過チャネル乱流（図 2 右）[4]，およびその一般化であるメッシュ・チャネル乱流（図 2 中央）[2]の検討例を紹介し，境界の透過性が流れや抵抗特性に及ぼす影響を検討した．
メッシュ・チャネル乱流においてメッシュ空隙サイズ S をゼロに近づけると古典的チャネル乱流に，Sとメッシュ素材幅 W を同時にゼロに近づける（非常に細い素材を用いて非常に目の細かなメッシュを組む）と完全透過チャネル乱流に，それぞれ漸近する
ことが数値実験によって示された．これは，完全透過チャネル乱流で観察された，古典的チャネル乱流とは質的にも大きく異なる力学機構 [4] が極限理想条件下でのみ生じる特異な存在ではなく，実在し得ることを意味する．
本報で示したメッシュ・チャネル乱流の直接数値
計算例は，計算領域および境界極近傍での空間解像度の影響をいくらか受けている可能性がある．これ
らの問題の克服，メッシュ厚み D の影響の把握，レイノルズ数依存性の理解など課題は多く，本研究はまだ緒に就いたばかりと言える．
謝辞：大規模計算の実施にあたって，クレイ・ジャ
パン・インクの鈴木 幸朗 様，および高橋 由樹 様から多大な協力を得ました．記して謝意を表します．

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

会議参加報告

SC16 参加報告
平石 拓 ∗∗

学術情報メディアセンター

ハイパフォーマンスコンピューティング（HPC）
分野の国際会議 SC161 が 2016 年 11 月 13 日から 18日にかけて米国ユタ州ソルトレイクシティで開催されました．

SC は，最近では登録参加者数が 10,000 人を超える超大規模な会議です．学術会議として HPC 分野のトップカンファレンスに位置付けられるだけでな
く，企業や大学・研究所等が出展する exhibition（見本市）も催される「HPC の祭典」的な性格もありま

図 1: SC16 の京都大学ブース

す．また，世界のスーパーコンピュータの演算性能を競う Top5002 の結果が発表される 1 年に 2 度の機会のひとつとしても知られています．
筆者を含む京都大学の教員も数名が参加し，今後の研究や次期スーパーコンピュータの調達などを見据えた情報収集および exhibition 展示などを行いました．京都大学の展示（図 1）では，学術情報メディアセンターが 2016 年 10 月にサービスを開始した新スーパーコンピュータ Camphor 2（Cray XC40），

12 月にサービスを開始した Laurel 2 (Cray CS400
2820XT)，Cinnamon 2 (Cray CS400 4840X) の紹介ならびに XC40 のコンピューティングボードの展示（図 2）を行いました．また，メディアセンターのスーパーコンピューティング研究分野，ビジュア
ライゼーション研究分野の研究成果である並列言語・
並列実装技術や可視化技術などをパネルや大型ディスプレイを用いて行いました．
京都大学以外の国内研究機関でも，理化学研究所の「京」や「ポスト京」，筑波大学・東京大学の Oakforest-

PACS，産業技術総合研究所の ABCI（AI Bridging
Cloud Infrastructure），東京工業大学の TSUBAME3.0
といった稼働中のシステムや開発プロジェクトの紹介が行われていました．
1 http://sc16.supercomputing.org
2 https://www.top500.org/

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

図 2: Cray XC40 のコンピューティングボード

SC16 での Top500 の 1 位は中国の Sunway Taihulight（神威太湖之光）が 2 回連続で獲得し，4 位までのランキングにも変動がありませんでした．6 位には筑波大学・東京大学の新システム Oakforest-PACS
がランクインしています．京都大学の Camphor 2 は33 位（国内では 4 位）にランクインしました．また，性能あたりの消費電力を競う Green500 では，Camphor 2 は 9 位にランクインしています．
次回の SC17 は，11 月 12 日～17 日に米国コドラド州デンバーで開催されます．SC17 にも引き続き参加し，ブース出展を行う予定です．

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

サービス再開日時

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

0:00

2016/04/01

13:00

13.00

2016/06/07

9:00

2016/06/08

9:30

24.50

2016/08/15

17:00

2016/10/01

0:00

1111.00

2016/05/19

サービス再開日時
4:21

2016/05/19

8:35

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

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

707.0017,16841,672

10,668,200

9,510,460

801.5

62

%

5月

739.77

17,434

68,466

12,322,300

10,913,500

871.5

59

%

6月

695.5020,25567,29315,000,500

12,459,900

940.0

72

%

7月

744.0021,757

79,478

16,640,600

14,630,600

940.0

72

%

8月

353.005,172

26,9527,840,680

7,027,830

940.0

70

%

9月

0.00

0

0

0

0

0.0

0

%

計

3239.27

81,786

283,861

62,472,280

54,542,290

748.8

56

%

ノード
利用率

※システム更新に伴い、8 月後半および 9 月はサービス休止

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

4月

5月

6月

7月

8月

9月



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

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



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

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

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

サービス再開日時

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

0:00

2016/04/01

13:00

13.00

2016/06/07

9:00

2016/06/08

9:30

24.50

2016/08/15

17:00

2016/10/01

0:00

1111.00

障害発生日時

サービス再開日時

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

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

707.0043,914

250,318

3,341,310

253,938

539.1

62

%

5月

744.00111,675

267,246

4,000,350

307,622

561.8

60

%

6月

695.50156,066

230,240

3,388,670

264,601

587.0

53

%

7月

744.00134,761

330,487

4,024,320

319,988

588.9

63

%

8月

353.0040,454

314,590

2,323,610

182,678

587.8

59

%

9月

0.00

0

0

0

0

0.0

0

%

計

3243.50

486,870

1,392,881

17,078,260

1,328,827

426.1

49

%

ノード
利用率

※システム更新に伴い、8 月後半および 9 月はサービス休止

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

4月

5月

6月

7月

8月

9月



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

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



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

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

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

サービス再開日時

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

0:00

2016/04/01

13:00

13.00

2016/06/07

9:00

2016/06/08

9:30

24.50

2016/08/15

17:00

2016/10/01

0:00

1111.00

障害発生日時

サービス再開日時‐なし

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

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

707.00808

893

23,024

1,019

10.1

12

%

5月

744.00474

3,592

30,962

8,641

11.0

25

%

6月

695.50

1,789

8,330

86,21250,227

12.0

35

%

7月

744.00

1,697

9,378

96,48256,275

12.0

43

%

8月

353.00861

5,573

53,897

38,761

12.0

34

%

9月

0.00

0

0

0

0

0.0

0

%

計

3243.50

5,629

27,766

290,577

154,923

9.5

25

%

ノード
利用率

※システム更新に伴い、8 月後半および 9 月はサービス休止

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

4月

5月

6月

7月

8月

9月



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

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



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

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

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

サービス再開日時

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

0:00

2016/04/01

13:00

13.00

2016/06/07

9:00

2016/06/08

9:30

24.50

2016/08/15

17:00

2016/10/01

0:00

1111.00

障害発生日時

サービス再開日時‐なし

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

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

707.0014,07020,2752,453,220

117,803

416.0

33

%

5月

744.0017,105130,961

5,794,490

284,382

416.0

67

%

6月

695.5016,087139,291

5,784,630

144,771

416.0

69

%

7月

744.0023,843140,199

6,102,690

215,745

415.9

74

%

8月

353.0013,74365,1733,350,530

143,932

415.8

72

%

9月

0.00

0

0

0

0

0.0

0

%

計

3243.50

84,848495,899

23,485,560

906,633

346.6

52

%

ノード
利用率

※システム更新に伴い、8 月後半および 9 月はサービス休止

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

4月

5月

6月

7月

8月

9月



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

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



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

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

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

サービス再開日時

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

0:00

2016/04/01

13:00

13.00

2016/06/07

9:00

2016/06/08

9:30

24.50

2016/08/15

17:00

2016/10/01

0:00

1111.00

障害発生日時2016/08/09

サービス再開日時

16:52

2016/08/10

04:57

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

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

4月

707.00

3,513

2,653

1,159,560

26,778

233.5

70

%

5月

744.00

1,451

4,764

1,622,990

17,877

258.9

84

%

6月

695.50

3,004

6,556

1,397,190

22,762

288.0

79

%

7月

744.004,470

28,511

1,977,640

90,546

287.4

83

%

8月

340.92

939

7,295

896,538

62,367

287.8

84

%

9月

0.00

0

0

0

0

0.0

0

%

計

3231.42

13,377

49,7797,053,918

220,330

225.9

67

%

ノード
利用率

※システム更新に伴い、8 月後半および 9 月はサービス休止

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

4月

5月

6月

7月

8月

9月



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

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



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

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

汎用コンピュータシステムのサービス状況
1.

ホスティング・ホームページサービス利用状況

VMホスティングサービス
600

595

591

WEBホスティングサービス

588

585

582

579

ストリーミングサービス

550
500
450
400
350
300

290

290

277

275

272

271

250
200
150
100
50
0

2

2

2

3

2

2

平成28年4⽉

平成28年5⽉

平成28年6⽉

平成28年7⽉

平成28年8⽉

平成28年9⽉(平成28年4⽉から平成28年9⽉)

大型計算機システム利用承認件数について

平成28年9月末現在、大型計算機システムの利用件数は、2,894件となっています。

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

センター利用による研究成果―平成 27 年度計算機利用結果報告書から―
6.

【1301 医用生体工学・生体材料学】
1.

Kanehiro Matsuyama, Yasuhiro Ishidou,

homogeneous shear turbulence laden with

Yong-Ming Guo, Hironori Kakoi, Naohiro

finite-size particles： Journal of Turbulence,

Shinohara, Yuhei Yahiro, Ichiro Kawanura,Vol.16, No.10, pp.979 - 1010, 2015

Takao Setoguchi, Shingo Maeda, Setsuro【4402 気象・海洋物理・陸水学】

Komiya： Prediction of the Long-term Course

7.

Shige： Development of a non-hydrostatic

Finite element Analysis and Mid-term

atmospheric model using the Chimera gridRadiological Evaluation： The Orthopaedic

method for a steep terrain： AtmosphericResearch

Science Letters, Vol.114, No.November 2015,

Society

(ORS)

2015

Annualpp.109 - 114, 2015

石堂康弘, 松山金寛, 瀬戸口啓夫, 栫博則, 廣津

武村一史, 石岡圭一, 重尚一： 重合格子法により地形表現を行う気象モデルへの雲微物理過く有限要素解析によるセメントレスステムの長程の実装： 第 17 回非静力学モデルに関するワ

期経過予測： 第８８回日本整形外科学会学術総ークショップ予稿集, 20159.

武村一史, 石岡圭一, 重尚一： 重合格子法を

松山金寛, 石堂康弘, 郭永明, 他： 中期的 X 線

用いた非静力学モデルへの雲微物理過程の導学的評価に基づいたセメントレスステムの有限入： 日本気象学会 2015 年度秋季大会予稿集,

要素解析と長期的 X 線学的評価の検証： 第４2015
10. 武村一史, 石岡圭一, 重尚一： 重合格子法を６回日本人工関節学会, 2016
4.

8.

匡隆, 郭永明, 小宮節郎： X 線学的評価に基づ会, 20153.

Kazushi Takemura, Ishioka Keiichi, Shochiof the Cementless Femoral Stem Based on

Meeting, 20152.

M. Tanaka, D. Teramoto ： Modulation of

Kanehiro Matsuyama, Yasuhiro Ishidou,

用いた非静力学モデルによる山岳波の数値実

Yong-Ming Guo, Hironori Kakoi, Naohiro

験： 日本気象学会 2015 年度春季大会予稿集,

Shinohara, Yuhei Yahiro, Ichiro Kawamura,

2015

Takao Setoguchi, Shingo Maeda, Setsuro【4801 機能材料・デバイス】

Komiya ：Finite Element Analysis of

11. Genta

Cementless

Femoral

Stem

Based

on

Kawaguchi,

Tokutaro

Komatsu,

Mitsuhiko

Maesato,

Hiroshi

Kitagawa,

Radiological

Tatsuro Imakubo, Andhika Kiswandhi, David

Evaluations ： The Orthopaedic ResearchGraf, James S Brooks ： Unconventional

Society (ORS) 2016 Annual Meeting, 2016

Magnetic and Resistive Hysteresis in an

Mid-termandLong-term【2202 自然災害科学】

Iodine ‐ Bonded Molecular Conductor ：5.

Asano, K., H. Sekiguchi, T. Iwata, M. Yoshimi,

Angewandte Chemie International Edition,

T. Hayashida, H. Saomoto, and H. Horikawa：

Vol.54, No.35, pp.10169 - 10172, 2015

and12. Jared M Taylor, Tokutaro Komatsu, Shunattenuation in the Osaka sedimentary basin,

Dekura, Kazuya Otsubo, Masaki Takata,western Japan, during the 2013 Awaji Island

Hiroshi Kitagawa： The Role of a Threeearthquake

JournalDimensionally Ordered Defect Sub-lattice on

International, Vol.204, No.3, pp.1678 - 1694,

the Acidity of a Sulfonated Metal-Organic

2016

Framework ：Modelling

of

：

wave

propagation

Geophysical【4304 数理物理・物性基礎】

JOURNALOF

THE

AMERICAN CHEMICAL SOCIETY, Vol.137,

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

を用いた CFRP 薄膜の変形挙動の解析： 第 40

No.35, pp.11498 - 11506, 2015
13. Tokutaro Komatsu, Jared M. Taylor, Hiroshi

回複合材料シンポジウム 講演論文集, 201522. 福造晃，西川雅章，北條正樹： 不連続 CFRTP

Kitagawa ：Design of a Conducting

Metal-Organic

Framework:

Orbital-Levelの塑性力学特性の微視的解析： 第 1 回材料

Matching

MIL-140A

Derivatives ：WEEK 若手学生研究発表会 講演論文集, 2015

INORGANIC CHEMISTRY, Vol.55, No.2,

23. 古賀貢史，橋爪良輔，新玉重貴，西川雅章，北in

條正樹： CFRP 内部のボイドの観察と幾何形

pp.546 - 548, 2015
14. Genta

Kawaguchi,

Tokutaro

Komatsu,

Hiroshi Kitagawa ：Mitsuhiko

Maesato,状を考慮した力学特性の有限要素解析： 日本機

TatsuroImakubo,械学会 M&M2015 材料力学カンファレンス

Interplay between

講演論文集, 2015

spin-density wave and 3 d local moments

24. Masayuki Uranagase, Ryosuke Matsumoto：with random exchange in a molecularTension-compression asymmetry in uniaxial

conductor： Physical Review B, Vol.93, No.7,

deformation of a magnesium bicrystal with

pp.075124 - , 2016[-1100] symmetric tilt grain boundary ：Computational Materials Science, Vol.118,【4905 工学基礎】
15. Shugo

Yasuda,Synchronized

Ryoichi

Yamamoto ：Molecular-Dynamics

pp.124 - 132, 2016
【5004 流体工学】

simulation for the thermal lubrication of a

25. Pillai Abhishek, Tomoaki Kitano, Ryoichi

polymeric liquid between parallel plates ：Kurose, Satoru Komori： Direct Numerical

Computers and Fluids, Vol.124, pp.185 - 189,

Simulation of Combustion Noise in Open2016

Hydrogen Diffusion Jet Flames： In Proc. of【5001 機械材料・材料力学】
16. 永田恵輔，西川雅章，北條正樹： シェル要素にInternational Gas Turbine Congress 2015
(IGTC2015), 2015

よる CFRP 膜面の展開挙動の解析： 第 6 回日

【5401 金属物性】

本複合材料会議（JCCM-6) 講演論文集, 2015

26. T. Ichitsubo, S. Okamoto, T. Kawaguchi, Y.

17. 橋爪良輔，新玉重貴，北條正樹，西川雅章： ボ

Kumagai, F. Oba, S. Yagi, N. Goto, T. Doi, E.

イドを含む CFRP 構造要素の強度解析手法のMatsubara ： Toward "rocking-chair type"

検討： 第 6 回日本複合材料会議（JCCM-6) 講

Mg-Li dual-salt batteries： J. Mater. Chem. A,

演論文集, 2015

Vol.3, pp.10188 - 10194, 201518. 福造晃，西川雅章，北條正樹： 不連続 CFRTP

27. S. Okamoto, T. Ichitsubo, T. Kawaguchi, Y.

の塑性特性の実験と解析： 第 6 回日本複合材

Kumagai, F. Oba, S. Yagi, K. Shimokawa, N.

料会議（JCCM-6) 講演論文集, 2015

Goto, T. Doi, E. Matsubara： Intercalation

19. 古賀貢史，橋爪良輔，新玉重貴，西川雅章，北and push-out process with spinel-to-rocksalt

條正樹： X 線 CT によるボイドの定量化と複合

transition on Mg insertion into spinel oxides

材剛性・強度への影響の解析： 日本機械学会関

in magnesium batteries： Advanced Science,

西学生会平成 26 年度学生員卒業研究発表会

Vol.2, pp.1500072-1 - 1500072-9, 201528. Y. Hinuma, Y. Kumagai, F. Oba, I. Tanaka：

講演論文集, 2015
20. Masaaki Nishikawa, Ichiro Taketa, Kei Ueda,

Categorization of surface polarity from a

Micromechanical Finite

crystallographic approach： Computational

Element Analysis on Compressive Strength

Materials Science, Vol.113, pp.221 - 230, 2016of Carbon Fiber-Reinforced Plastics： Proc. of

29. K. T. Butler, Y. Kumagai, F. Oba, A. Walsh：

14th Japanese-European Symposium on

Screening procedure for structurally andComposite Materials, 2015

electronically matched contact layers for

Masaki Hojo ：21. 永田恵輔，西川雅章，北條正樹： モード合成法

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

high-performancesolar

cells:

hybridperovskites ： J. Mater. Chem. C, Vol.4,【5803 生物物理学】

pp.1149 - 1158, 201637. Hiraku

30. A. A. Belik, W. Yi, Y. Kumagai, Y. Katsuya, M.
Tanaka, F. Oba ：(Tl1-xScx)ScO3:

LiNbO3-type oxide

High-pressure

synthesis,

Oshima,Tomohiko

Hayashi,Masahiro Kinoshita ： Essential Roles ofWater in Actin-Myosin Binding： BiophysicalJournal, Vol.110, No.Issue 3, Supplement 1,

crystal structure, and electronic properties：

pp.300a - , 2016

Inorganic Chemistry, Vol.55, pp.1940 - 1945,【6801 化学系薬学】

2016

38. Masato

Saito,

NobuyaTsuji,

Yusuke【5501 化工物性・移動操作・単位操作】

Kobayashi, Yoshiji Takemoto ：Direct31. 岡本幸也, 高木洋平, 岡野泰則： 高シュミット

Dehydroxylative

Coupling

Reactionof

混合溶液における相分離ドメインと乱流渦構造

Alcohols with Organosilanes through Si-Xの相関： 日本流体力学会年会 2015 講演論文集,

Bond Activation by Halogen Bonding： Org.

2015

Lett., Vol.17, No.12, pp.3000 - 3003, 201532. N. Adkar, T. Yamamoto, Y. Takagi, Y. Okano,

39. Yusuke Kobayashi, Ryuta Kuramoto, Yoshiji

S. Dost： Proc. of the joint conference ofTakamoto ： Catalytic asymmetric formal

ISPS-6synthesis of beraprost： Beilstein J. Org.

and

ITTW2015 ：Numerical

Simulation of the Effect of the Free Surface

Chem, Vol.11, pp.2654 - 2660, 2015

Shape on Marangoni Convection in a
Ring-Shaped Liquid Film with Low Prandtl
Number, Vol.Fl-ST(A)-P07, 2015
33. N. Adkar, T. Yamamoto, Y. Takagi, Y. Okano：Numerical investigation of the effect of theliquid film volume on thermocapillary flow ina thin circular liquid film ~ various Prandtl
number fluids ~： 化学工学会第４７回秋季大会, Vol.P122, 2015
34. T. Yamamoto, Y. Takagi, Y. Okano, S. Dost：
Numerical

investigationofoscillatory

thermocapillary flows under zero gravity in acircularliquidsurfaces ：

film

with

concavefreePhys. Fluids, Vol.28, No.3,

pp.032106 - , 2016
【5603 地球・資源システム工学】
35. Yumi Katasho, Yunfeng Liang, Sumihiko
Murata,Yasuhiro

Fukunaka,

Toshifumi

Matsuoka, Satoru Takahashi： Mechanisms
forEnhanced

Hydrophobicity

byAtomic-Scale Roughness： Scientific Reports,
Vol.5, No.13790, pp.1 - 12, 2015
36. Asuka Fukutome, Haruo Kawamoto, Shiro
Saka：

Gas-Phase

ReactionsofGlyceraldehyde and 1,3-Dihydroxyacetone as
Models for Levoglucosan Conversion during
Biomass Gasification： ChemSusChem, 2016

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

資料：大型計算機システム利用負担金 別表

別表１ スーパーコンピュータシステム
提供サービス
利用負担額
コース

タイプ

セット

エントリ

タイプA
タイプB
タイプC
タイプE

基本

パーソナル

タイプA1

グループ

基本
基本
基本
基本
最小
追加単位

タイプA2

最小
追加単位

タイプA3

最小
追加単位

タイプB1

最小
追加単位

タイプB2

最小
追加単位

タイプB3
タイプC1
タイプC2
タイプE1
タイプE2
タイプE3
タイプA
タイプB

大規模ジョブタイプCタイプE
専用クラスタ

-

最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位

ストレージ容量追加
ライセンスサービス

システム バッチ

12,600 円/年
100,000 円/年
100,000 円/年
100,000 円/年
100,000 円/年
230,000 円/年
115,000 円/年
276,000 円/年
69,000 円/年
690,000 円/年
345,000 円/年
240,000 円/年
120,000 円/年
288,000 円/年
72,000 円/年
720,000 円/年
360,000 円/年
150,000 円/年
150,000 円/年
180,000 円/年
90,000 円/年
280,000 円/年
140,000 円/年
336,000 円/年
84,000 円/年
840,000 円/年
420,000 円/年
23,000 円/週(7日）
11,500 円/週(7日）
24,000 円/週(7日）
12,000 円/週(7日）
15,000 円/週(7日）
7,500 円/週(7日）
28,000 円/週(7日）
14,000 円/週(7日）
720,000 円/年
360,000 円/年
10,000 円/年
20,000 円/年

システム資源

経過時間 ストレージ
無料
（時間）
（TB)
利用者数

最大1ノード相当((36コア、128GBメモリ)×1)
1
最大4ノード相当((68コア、16+96GBメモリ)×4)
最大4ノード相当((36コア、128GBメモリ)×4)
B
168
最大1ノード相当((72コア、3072GBメモリ)×1)
C
最大2ノード相当((10コア、32GBメモリ+ 1MIC)×2)
E
4ノード((68コア、16+96GBメモリ)×4)
優先
2ノード((68コア、16+96GBメモリ)×2)
8ノード((68コア、16+96GBメモリ)×8)A
準優先
2ノード((68コア、16+96GBメモリ)×2)
8ノード((68コア、16+96GBメモリ)×8)
占有
4ノード((68コア、16+96GBメモリ)×4)
4ノード((36コア、128GBメモリ)×4)
優先
2ノード((36コア、128GBメモリ)×2)
8ノード((36コア、128GBメモリ)×8)
B
準優先
2ノード((36コア、128GBメモリ)×2)
8ノード((36コア、128GBメモリ)×8)
336
占有
4ノード((36コア、128GBメモリ)×4)
1ノード((72コア、3072GBメモリ)×1)
優先
1ノード((72コア、3072GBメモリ)×1)
C
2ノード((72コア、3072GBメモリ)×2)
準優先
1ノード((72コア、3072GBメモリ)×1)
4ノード((10コア、32GBメモリ+ 1MIC)×4)
優先
2ノード((10コア、32GBメモリ+ 1MIC)×2)
8ノード((10コア、32GBメモリ+ 1MIC)×8)E
準優先
2ノード((10コア、32GBメモリ+ 1MIC)×2)
8ノード((10コア、32GBメモリ+ 1MIC)×8)
占有
4ノード((10コア、32GBメモリ+ 1MIC)×4)
8ノード((68コア、16+96GBメモリ)×8)A
占有
4ノード((68コア、16+96GBメモリ)×4)
8ノード((36コア、128GBメモリ)×8)
B
占有
4ノード((36コア、128GBメモリ)×4)
168
2ノード((72コア、3072GBメモリ)×2)
C
占有
1ノード((72コア、3072GBメモリ)×1)
8ノード((10コア、32GBメモリ+ 1MIC)×8)
E
占有
4ノード((10コア、32GBメモリ+ 1MIC)×4)
8ノード((36コア、128GBメモリ)×8)
B
4ノード((36コア、128GBメモリ)×4)
ストレージ容量10TBの追加につき可視化ソフト（AVS,ENVI/IDL)およびプリポストウェアの1ライセンスにつき
B

A

共有
共有
共有
共有
共有

0.2
2.0
2.0
2.0
2.0
16.0
8.0
19.2
4.8
32.0
16.0
16.0
8.0
19.2
4.8
32.0
16.0
16.0
16.0
19.2
9.6
16.0
8.0
19.2
4.8
32.0
16.0

-

32.0
16.0

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコースを、
年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除した後、利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． 大型計算機システムの全ての利用者は、上記表のサービスの他、次のサービスを受けることができる。
1） 大判プリンタサービス
2） その他、大型計算機システムが提供するサービス、機器の利用
3． 上記表の大規模ジョブコース、ストレージ容量追加、ライセンスサービスの申請には、スーパーコンピュータシステムの利用者であることが必要である。
4． 「共有」： 当該カテゴリのユーザ間で一定の計算資源を共有するベストエフォートのスケジューリングを行う。
「準優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/4の計算資源が確保されることを保証する。
「優先」： 定常稼働状況において記載値（以上）の計算資源が確保されるように優先スケジューリングを行う。
また、稼働状況によらず記載値の1/2の計算資源が確保されることを保証する。
「占有」： 稼働状況によらず記載値の計算資源が確保されることを保証する。
5． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。
6． グループコース及び専用クラスタコースの利用者番号は利用者あたり年額5,000円を負担することで追加できる。
7． 機関・部局定額制度
他機関又は学内における部局（『国立大学法人京都大学の組織に関する規程』第3章第2節から第11節で定める組織をいう。）の組織が、その組織単位でグループコースサービスを利用申請する場合の利用負担額は、別表１に規定するの1.5倍の額とする。なお、利用負担額が年額150万円未満の場合は100人、年額150万円を超える場合は、150万円毎に100人までの利用者を認める。ストレージは、1.5倍の容量とする。
8． スパコン連携サービス
学術情報メディアセンターのスーパーコンピュータシステムと密な連携により、学内における部局の組織が計算サーバ等を設置する
場合、下記の負担額を支払うものとする。
冷却方式
水冷
空冷

利用負担額
11,300 円/月
13,200 円/月

利用負担額算定単位
水冷冷却方式の計算サーバ等の定格電力 1kWにつき
空冷冷却方式の計算サーバ等の定格電力 1kWにつき

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

-

8
4
16
4
16
8
8
4
16
4
16
8
8
8
16
8
8
4
16
4
16
8

-

16
8

別表２（汎用コンピュータシステム）
区分

利用負担額

単位

VM ホスティングサービス

72,000 円／年 1 仮想マシンにつき

ホームページサービス

6,000 円／年 1 ドメイン名につきストリーミングサービス

6,000 円／年 1 申請につき備考
1. 利用負担額は、総額表示である。
2. 上記表の汎用コンピュータシステムのサービスを利用するためには、大型計算機シ
ステムの利用者であることが必要である。
3. VM ホスティングサービスにおいて、下記の負担額を支払うことにより CPU、メモリ、ディスクを増量することができる。
区分

利用負担額

単位

CPU 増量18,000 円／年 2 コアにつき(最大 8 コアまで)

メモリ増量

18,000 円／年 8GB につき(最大 64GB まで)

ディスク増量18,000 円／年 200GB につき(最大 1,000GB まで)

4. VM ホスティングサービスにおいて VMware を用いる場合は、下記の負担額を支払うことにより VMware の利用及び CPU、メモリ、ディスクを増量することができる。ただし、システム資源が限られているためサービスの提供を限定することがある。
区分

利用負担額

単位

VMware 利用

72,000 円／年 1 仮想マシンにつき

CPU 増量36,000 円／年 2 コアにつき(最大 8 コアまで)

メモリ増量

36,000 円／年 8GB につき(最大 64GB まで)

ディスク増量18,000 円／年 200GB につき(最大 1,000GB まで)

5. ホームページサービス及びストリーミングサービスにおいて、下記の負担額を支払うことにより公開スペースの上限を拡大することができる。
区分

利用負担額

公開スペース 20GB プラン

3,000 円／年

公開スペース 50GB プラン

9,000 円／年

6. 利用負担額は、当該年度(4 月から翌年 3 月まで)の利用に対して年額として算定す
るが、年度途中から利用を開始する場合には月数に応じて減額する。

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

※民間機関利用

別表３ スーパーコンピュータシステム
システム

A

B

C

E

経過時間
（時間）

ストレージ
（TB)

無料
利用者数

8ノード(68コア、16+96GBメモリ)×8)

336

19.2

16

1,104,000 円/年

12ノード(68コア、16+96GBメモリ)×12)

336

28.8

24

1,656,000 円/年

16ノード(68コア、16+96GBメモリ)×16)

336

38.4

32

2,208,000 円/年

8ノード(36コア、128GBメモリ)×8)

336

19.2

16

1,152,000 円/年

12ノード(36コア、128GBメモリ)×12)

336

28.8

24

1,728,000 円/年

16ノード(36コア、128GBメモリ)×16)

336

38.4

32

2,304,000 円/年

2ノード(72コア、3072GBメモリ)×2)

336

19.2

16

720,000 円/年

3ノード(72コア、3072GBメモリ)×3)

336

28.8

24

1,080,000 円/年

4ノード(72コア、3072GBメモリ)×4)

336

38.4

32

1,440,000 円/年

8ノード(10コア、32GBメモリ+ 1MIC)×8)

336

19.2

16

1,344,000 円/年

12ノード(10コア、32GBメモリ+ 1MIC)×12)

336

28.8

24

2,016,000 円/年

16ノード(10コア、32GBメモリ+ 1MIC)×16)

336

38.4

32

2,688,000 円/年

システム資源利用負担額

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコースを、
年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除した後、利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。

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

資料：サービス利用のための資料一覧―
1.

サービス利用のための資料一覧―

スーパーコンピュータシステム・ホスト一覧
システム A：camphor.kudpc.kyoto-u.ac.jp



システム B・C：laurel. kudpc.kyoto-u.ac.jpシステム B（SAS 利用時）：sas.kudpc.kyoto-u.ac.jp

システム E：camellia.kudpc.kyoto-u.ac.jp※ ホストへの接続は SSH(Secure SHell) 鍵認証のみ、パスワード認証は不可

2.

問い合わせ先 & リンク集


情報環境機構のホームページ
http://www.iimc.kyoto-u.ac.jp/



学術情報メディアセンターのホームページ
http://www.media.kyoto-u.ac.jp/



スーパーコンピュータシステムに関する問い合わせ先

利用申請などに関する問い合わせ先【情報環境支援センター】
E-mail : zenkoku-kyo@media.kyoto-u.ac.jp / Tel : 075-753-7424
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/

システムの利用など技術的な問い合わせ先【スーパーコンピューティング掛】
E-mail : consult@kudpc.kyoto-u.ac.jp / Tel : 075-753-7426
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/contact.html



ホームページ・ホスティングサービス関する問い合わせ先【クラウドコンピューティング掛】
E-mail : whs-qa@media.kyoto-u.ac.jp / Tel : 075-753-7494
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/whs/

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

編集後記

編 集 後 記今号の SC16 参加報告記事を担当しました．この会議には，日本の大学・企業・報道関係者が数多く参加しており，日本語の参加報告記事も数多く公開されています．「SC16 ニュース」などでウェブ検索
すると，写真つきの読み易い記事もいくつか見つかります．スパコン業界のトレンドを掴むのにもよ
いかと思いますのでご興味のある方は読んでみてください．
(副部会長)

京都大学学術情報メディアセンター全国共同利用版広報 Vol．15，No．2
2017
年 ３月
3 月３０日
xx発行
日発行
２
０１７年
広報編集部会
編集者 京都大学学術情報メディアセンター

深沢 圭一郎（部会長）

全国共同利用版広報編集部会

平石 拓（副部会長）

発行者 〒606-8501 京都市左京区吉田本町

檀原 正憲

南雲 円

京都大学学術情報メディアセンター

尾形 幸亮

高見 好男

Academic Center for Computing and Media Studies

元木 環

Kyoto University
Tel. 075-753-7400
http://www.media.kyoto-u.ac.jp/
印刷所 〒616-8102 京都市右京区太秦森ヶ東町 21-10
株式会社エヌジーピー

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

表紙デザイン：谷 卓司（ティアンドティ･デザインラボ）