巻頭言

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

牛島 省

本号では，京都大学学術情報メディアセンターにおいて実施された平成 29 年度の共同研究報告について特集いたします．平成 29 年度は，
「若手・女性研究者奨励枠」では 21 件，「プログラム高度化共同研究」では 4 件の共同研究が実施されました．
学術情報メディアセンターでは，スーパーコンピュータ共同研究制度として，「若手・女性研究者奨励枠」，「大規模計算支援枠」，「プログラム高度化共同研究」という３種類の研究支援事業を実施しています．課題募集は年度開始前後に行われます．また，年度途中で追加募集が実施される場合があります．応募された課題の採否は，スーパーコンピュータ共同研究企画委員会にて審査されます．本共同研究制度への応募方法や申込み時期等の詳細につきましては，http://www.iimc.kyoto-u.ac.jp/ja/services/comp/support/kyodo/
をご参照下さい．
上記の共同研究制度のうち，「若手・女性研究者奨励枠」の応募資格は，(1) 40 歳未満の若手研究者 (学生を含む, 性別は問わない)，あるいは，(2) 女性研究者 (年齢は問わない)，と平成 27 年度より拡大されており，スーパーコンピュータを利用することで学術的にインパクトがある成果を創出できる課題に対して，計算機利用負担金の全額または一部を本センターが負担しています．さらに，「若手・女性研究者奨励枠」は，平成 29 年度より JHPCN (学際大規模情報基盤共同利用・共同研究拠点) の活動の一環となり，優れた
課題は，JHPCN の萌芽型共同研究へ推薦されることになりました．推薦された課題は，
JHPCN シンポジウムにてポスター発表を行うことが可能となるなどの特典がありますので，是非この共同研究制度をご活用ください．本号に掲載された 21 件の「若手・女性研
究者奨励枠」の研究成果報告では，各種の流体および流体構造連成問題，破壊力学，量子化学，分子動力学計算，MRI データ解析，などに関する多様な研究が行われ，全国 9 大
学等の若手・女性研究者の研究をサポートいたしました．なお，JHPCN の萌芽型共同研
究には 6 件の課題が採択されました．
また，「プログラム高度化共同研究」では，京都大学のスーパーコンピュータをグループ
コースまたは専用クラスタコースでご利用の研究グループ (JHPCN および HPCI にて京都大学のスーパーコンピュータを利用する研究者も含む) を対象に，プログラムの制御構
造・データ構造の改良による性能チューニングや並列化手法の改良，問題分割・負荷分散方式などの改良による大規模計算プログラムの高度化・高性能化などを支援しています．
平成 29 年度は，(1) 地盤の変形・破壊を扱う有限要素解析，(2) 高分子分離膜の大規模解
析，(3) 非均質異方性材料に対するボクセル有限要素解析，(4) マイクロ波放電式中和器の 3 次元プラズマ粒子計算，に関する 4 件のプログラム高度化支援を行いました．これらの共同研究では，並列 FFT ライブラリの最適化，並列 IO の導入や MPI 通信の改善によるスケーラビリティの改善，行列解法の高速化などが実施され，約 36 倍の高速化が実現
した課題もあります．このように，多分野にわたる解法の高度化・高性能化が行われており，計算効率を向上させる際の参考になると思います．
今後も皆様の研究，教育にご活用いただけるようにセンター教職員も尽力していきます
ので，ご利用・ご支援のほど、よろしくお願い申し上げます．

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

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

高精度視覚質感記憶の心理学的基盤と神経機構の解明

藤道 宗人
京都大学大学院人間・環境学研究科共生人間学専攻

1 問題・目的
質感とは物体そのものがもつ表面の性質をさす。
ヒトは物体がもつ様々な質感を正確に把握することで，その物体に正しく働きかけることができる。
それではヒトの質感認知はどのような脳内メカニ
ズムによって成立するのだろうか。
これまでの研究によりヒトの質感知覚に腹側視覚経路が関与することが明らかとなった。ところ
が，腹側視覚経路の中でどの領域が関連するのか

2.1 参加者
24名の京都大学生および大学院生 (男性15名，女性 9 名)が実験に参加した。平均年齢は 22.9 歳(SD=2.8)であった。全ての参加者は正常な視力を有していた。実験は京都大学こころの未来研究センターの MRI スキャナー及び関連実験設備を用いて行った。

2.2 ローカライザー課題については見解が分かれている。例えば腹側視覚領域の細分化のために顔，シーン，物体，そし経路の中でシーンの認知に重要な役割を担う海馬て色条件を設定した。さらにそれぞれの条件に属傍回が関与するという知見[1]や，顔の認知に関わする画像をスクランブルにした 4 条件の計 8 条件

る紡錘状回が質感の認知にも関連することを示すを設定した。課題はブロックデザインで行われ，知見[2]がある。このように，質感知覚に腹側視覚各条件に属する刺激が集中的に連続呈示 (1 ブロ経路が関与することは明らかとなっているが，具ック)された。参加者はブロックの最初に画像が呈体的な領域の特定には至っていない。さらに，こ示されたらボタン押しによって反応した。

れらの研究は認知の一側面である知覚処理に着目しており，質感の視覚記憶メカニズムは明らかになっていない。例えばヒトの知覚と視覚記憶に関わる脳領域には共通性が見られるが[3]，質感認知においてはどうだろうか。本研究は質感を知覚して記憶する脳内メカニズムを検討することを目的
とした。

2 方法
実験は頭頂間溝の同定および腹側視覚経路内の領域を細分化するための課題 (ローカライザー課
題)と質感の視覚記憶メカニズムを検討するための課題 (視覚記憶課題)で構成された。

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

2.3 視覚記憶課題
実験は粗さ記憶課題と光沢感記憶課題で構成された。参加者は継時呈示される 2 枚のサンプル球
体画像を見て，どちらの課題を行っているかを識別する必要があった（弁別フェイズ）。例えば 2枚のサンプルが粗さの次元で変化していれば，粗さ記憶課題であった。2 枚のサンプル球体画像の呈示後に保持しておくべきサンプルを指示した（保持フェイズ）。
参加者は試行の最後に呈示されるテスト球体画像と保持した画像を比較して，粗さの程度が小さい，または光沢感の強い方の画像を選択することが求められた。照明環境に依存しない質感の視覚記憶を検討するために，2 枚のサ

頭頂間溝はどちらのフェイズでもチャンスレベル

ンプルとテスト球体画像の照明環境は常に異なっ以上だった (t(19) = 2.26, p < .05; t(19) = 4.23, pていた。

< .01)。さらにこれらの領域とフェイズに交互作用があることも明らかになった (F(1, 19) = 4.76,

3 解析

p < .05, ηp2 = .2)。これらの結果は，保持期間中に果たす役割が領域によって異なっていることを示

ローカライザー課題によって同定した頭頂間溝している。すなわち，ヒトが質感を記憶するとき

および腹側視覚経路内の各領域において，課題情には，記憶すべき情報を腹側高次視覚野内で弁別

報 (粗さ・光沢感記憶課題)が保持されているかをして，その情報を頭頂間溝で保持している可能性

マルチボクセルパターン解析 (MVPA)によってが示唆された。

検討した。具体的には，まず保持期間中の脳活動に対して MVPA を適用し，質感の視覚記憶に関連する領域を同定した。そして保持期間中のタイムコースにそって MVPA 成績を算出して，それらの領域が弁別フェイズと保持フェイズでどのように変化するのかを検討した。

4 結果・考察保持期間中の平均 MVPA 成績を図 1 に示す。
腹側視覚経路の中で顔の認知に関わる領域（以下，
腹側高次視覚野）でのみチャンスレベル以上の成図 1 保持期間中の平均 MVPA 成績。黒点線はチャンスレベルを示す。

績だった(t(19) = 2.53, p < .05)。この結果は紡錘状回が質感認知に関連するという結果[2]を支持す
るものとなった。これは着目した質感の特徴レベルの違いに起因すると考えられる。具体的には，
本研究では粗さと光沢感，[2]では革や鉄など質感
カテゴリ内におけるそれぞれの性質を題材としていた。一方[1]では物体特徴の 1 つとして質感（テクスチャ）を取り上げて，テクスチャと色あるいは形を比較した。ただし，紡錘状回と海馬傍回は相補的関係にあることも知られているため，紡錘状回は低次の質感情報に，海馬傍回は高次の質感

図 2 フェイズごとの MVPA 成績。黒点線はチャンスレベルを示す。

情報に関連する可能性も考えられる。加えて頭頂間溝もチャンスレベル以上で課題識別が可能であった (t(19) = 4.30, p < .01)。したがって，質感知

5 参考文献覚と視覚記憶に共通する腹側高次視覚野をベース[1] Cant, J. S., et al. (2009). Exp Brain Res, 192,

として，視覚記憶を担う頭頂間溝[4]との関連のも

391–405.

とで質感の視覚記憶が成立すると考えられる。

[2] Hiramatsu, C., et al. (2011). NeuroImage, 57,

腹側高次視覚野と頭頂間溝のフェイズごとの482–494.

MVPA 成績を図 2 に示す。MVPA 成績の推移は[3] Harrison, S. A., & Tong, F. (2009). Nature, 458,

領域間で異なる様相を示していた。腹側高次視覚632–658.

野では弁別フェイズでのみチャンスレベル以上の[4] Katherine, C. B., & Yaoda, X. (2016). Nature

MVPA 成績だった (t(19) = 2.97, p < .01)。一方のNeurosci, 19, 150–157.

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

阿蘇山の複雑地形における高精度メッシュを用いた領域大気モデルの感度解析及び性能評価
石田 桂
熊本大学大学院先端科学研究部

1 はじめに2 領域大気モデル

平成 24 年 7 月九州北部豪雨において AMeDAS

本 研 究 で は 領 域 大 気 モ デ ル と し て The

阿蘇乙姫観測所で 816.5mm という記録的な降水

Weather Research and Forecasting Model

量が観測された．阿蘇山は中央火口丘を中心に外(WRF)を用いた．計算領域にはネストした複数の輪山に囲まれたカルデラ地形を形成している（図ものを使用した．最内の計算領域の解像度を南北

1）．外輪山とカルデラ内部との標高差は 200～

25km，東西 18km のカルデラを覆うように設置

600m 程で，外輪山の斜面は切り立っており阿蘇する．最内側の計算領域の解像度に合わせ，他のカルデラは非常に複雑な地形となっている．平成

計算領域の解像度及び大きさの調節を行った．鉛24 年 7 月九州北部豪雨における阿蘇カルデラ内

直方向の分割数は 58 とした．ただし，阿蘇山のに位置する阿蘇乙姫観測所での記録的な降水量は，

標高が 1,592 m であり，外輪山の標高はそれ以下

高湿度な大気が強い西南西の風に乗り西側の外輪であるため，
高度2,000 m 以下を細かく分割した．

山開口部からカルデラ内に侵入し，外輪山に衝突

領域大気モデルでの再現計算では，大気再解析デしたときの地形効果により生じたと考えられていータから初期条件及び境界条件を得る必要がある．

る．そこで，本研究では水平・鉛直方向に高解像

大気再解析データとして ERA Interim 及び度メッシュを用いることにより，領域大気モデル

NCEP FNL を用いた．それぞれ，水平解像度はによる平成 24 年 7 月九州北部豪雨の再現計算を

おおよそ 75 km，100 km である．WRF はデフ行い豪雨時の大気の流れを解析すること最終目標

ォルトの標高データは最も詳細なものでも 1

として，領域大気モデルの感度解析及び性能評価

km 程度の解像度である．そこで，阿蘇カルデラを行う．

の複雑な地形を反映させるために，国土地理院発行の数値地図 50m メッシュ標高データにより
WRF にデフォルト標高データを置き換えた．
WRF にはパラメタリゼーションとして雲微物
理過程，積雲過程，地表面過程，放射過程，惑星

外輪山

阿蘇山

境界層過程などがあり，それぞれに複数のスキー
ムが用意されている．本研究では領域内の降水量に特に影響が大きいと思われる雲微物理過程及び惑星境界層過程に注目する．
上記の設定条件の元，
初期・境界条件として用いる大気再解析データ，

図 1．阿蘇カルデラ（国土地理院アナグリフを地位院地図を用いて描
画．鉛直方向を 5 倍に強調．
）

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

最内側の計算領域の解像度，及びパラメタリゼー
ションオプションに関し感度解析を行う．
領域大気モデルによる再現計算を行うために，

まず，計算領域を準備した．最内の計算領域の解の再現計算の結果を図 2 に示す．解像度 100 m で像度を 100 m と 200 m の 2 セットを準備し，国の再現計算でも平成 24 年 7 月九州北部豪雨時に

土地理院発行の数値地図 50m メッシュ標高デー

阿蘇カルデラ内で生じた強い降水帯の位置を正確

タにより，モデル内の標高データを置き換えた．

に再現することはできなかった．一方で，実際のしかしながら，200 m の解像度では複雑な阿蘇カ降雨帯とは位置がずれているものの，解像度 100

ルデラ地形の再現性が低かったため，100 m の計

m の計算領域を用いた再現計算では解像度 400

算領域を採用した．結果，解像度を 25.6 ,m，6.4

m の計算領域では再現できない強い降水帯の再km，1.6 km，400 m，100 m とするネストした 5

現ができている．このことから，本豪雨時に阿蘇

つの計算領域を用いることとした．

カルデラの複雑地形が豪雨の降水量に影響を与えた可能性，また，高解像度の計算領域を用いるこ

3 感度解析及び性能評価とにより降水帯の再現性を向上できる可能性が示された．

 2 種類の大気解析データ（ERA Interim 及び次に平成 24 年 7 月九州北部豪雨の風の流れをNCEP FNL）を用い，領域大気モデルのパラメタ

図 3 に示す．豪雨時の風の流れは，当初考えられ

リゼーションオプションに関し感度解析を行った．

ていた外輪山開口部に流入したものだけでなく，

ここで，5 つの計算領域を用いた計算は非常に多カルデラの複雑な地形が風の流れを変化させ，大くの時間がかかるため，25.6 km，6.4 km，1.6 km

気の収束が発生し，本豪雨が生じさせた可能性がの 3 つの計算領域を用い，広域の降水量の再現性示された．

から各パラメタリゼーションオプションのスキー
ムの絞り込みを行った．ERA Interim を初期・境

4 おわりに

界条件に用いた計算では降水量の再現性が非常に低かったため，NCEP FNL に注目することとした．

 降水帯の再現性は不十分であるものの，スーパーコンピュータの使用により可能となる詳細な

 3 つの計算領域を用いた感度解析の結果から，再現計算により，複雑な地形の影響を受ける豪雨

3 組のパラメタリゼーションオプションの組み合

事象のより現実に近い再現が可能となることが示わせを選定し，5 つ全ての計算領域を用いた再現された．また，そのときの複雑な大気の流れを解計算を行った．ここで，全ての計算領域に関し初析することが可能となった．今後，領域大気モデ期条件を大気再解析データから得る通常の方法でルの設定条件などを改良していくことにより，豪は，解像度のギャップから計算が不安定になり，雨時の降水帯の再現性を高めることができると考再現計算を行うことはできなかった．そのため，えられる．

400 m までの再現計算を行い，その計算結果から
100 m の計算領域の初期・境界条件を取得した．
 結果の一例として，同じパラメタリゼーショ
ンオプションを用いた解像度 400 及び 100 m で

図 2． 400 m（上），100 m（下）の計算領域での
平成 24 年 7 月九州北部豪雨の降水量の再現計算結果

図 3．豪雨時の風の流れ（鉛直方向を 10 倍に強調）

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

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

1

はじめに
土木構造物の防災・減災の観点から，建設材料として使用される鉄筋コンクリートのひび割れ進
展挙動を定量的に予測することが重要な課題となっている．鉄筋コンクリートでは，鉄筋とコンクリートの相互作用やひび割れによる応力の再分配などにより，応力状態が複雑に変化するため，ひび割れが開口するだけでなく，形成されたひび割
れ面が閉口し，摩擦接触が生じる．そのため，鉄筋コンクリートのひび割れ進展挙動を精度よく再現するためには，ひび割れの開口挙動に加えて摩図1 損傷と摩擦を考慮したせん断応力

擦接触まで考慮した数値解析手法を開発する必要がある．

を表現する．接触方向の応力は次式で表される．

そこで本研究では，準脆性材料の破壊力学に基 'x > 0

'x = (1−D)C11'x + (1−D)C12'y(2)づく等方性損傷モデル [1] に摩擦接触のモデル化

 'x ≤ 0

'x = C11'x + (1−D)C12'y(3)を導入し，摩擦接触を伴うひび割れ進展挙動を再ここで，'x > 0 はひび割れ面が開口している状態，

現できる損傷モデルを提案する．そして，既往の'x ≤ 0 はひび割れ面が接触している状態である．

一面せん断試験の実験結果と同条件の解析結果を接触している際に生じる摩擦力は，次式で表さ

比較し，提案モデルの妥当性を検証する．

2

れるクーロンの摩擦法則により算出する．

f = n =  |'x|(4)

ここで，n は垂直応力， は摩擦係数である．

数値解析手法

摩擦応力と局所座標系のせん断応力を比較し，ひび割れの開口挙動を表現するため，準脆性材

すべりを判定する．すべらない場合，ひび割れ面

料の破壊力学に基づく等方性損傷モデルの定式化が固着した状態となるため，せん断応力は弾性応を使用する．構成則は次式で表される．

力とする．すべる場合，図1 に示すような損傷と = (1−D)C

(1)

摩擦を考慮したせん断挙動を定式化する．局所座

ここで， は応力ベクトル， はひずみベクトル，

標系のせん断応力は次式で表される．

C は弾性係数マトリックス，D は 0 ～ 1 の値を取

f ≥ |'xy|

'xy = C33'xy

る損傷変数である．

f ≥ |'xy|

'xy = (1−D)C33'xy + sgn('xy)Df (6)(5)

接触方向に対応した局所座標系の構成則を使用

ここで，f ≥ |'xy| はひび割れ面が固着している状することで，ひび割れ面における摩擦接触を定式

態，f ≥ |'xy| はひび割れ面に沿ってすべっている化する．接触方向の剛性を保持させることで接触状態である．

6 − Academic Center for Computing and Media Studies, Kyoto University
図3 実験と解析のせん断応力－せん断変位関係
図2 有限要素モデル

3

妥当性の検証3.1 解析条件
提案モデルの妥当性を検証するため，van derPluijm の一面せん断試験 [2] の実験結果と解析結果を比較する．図2 に試験体寸法を基に作成した有限要素モデルを示す．要素には四角形要素を用いた．節点数は 2790，要素数は 2720 である．材
料パラメータは，レンガが，E = 17 GPa， = 0.15，
モルタルが，E = 5 GPa， = 0.2，界面が，E = 1 GPa，

 = 0.2，k = 10，Gf = 0.01 N/mm，0 = 0.0006， = 0.8

図4 変形および損傷分布 (垂直荷重：1.0 MPa)

4

おわりに

である．実験ではレンガとモルタルの界面が破壊

本研究では，ひび割れの開口挙動に加えて摩擦したことから，解析では界面のみ損傷が生じる条接触を考慮した損傷モデルを提案し，その妥当性件を与えた．水平方向に強制変位 0.8 mm を 400を検証した．今後は，鉄筋コンクリートの破壊挙

step で与え，垂直方向に 0.1, 0.5, 1.0 MPa の圧縮荷動に提案モデルを適用し，大規模シミュレーショ

重を一定に与えた 3 ケースの解析を実施した．

ンを実施することで，摩擦接触を伴うひび割れ進

3.2 解析結果

ある．

図3 に実験と解析におけるせん断応力－せん

展挙動を詳細に再現し，有効性を検討する予定で

参考文献

断変位関係を示す．提案モデルは，実験結果にお[1] 車谷麻緒，寺田賢二郎，加藤準治，京谷孝史，樫

ける破壊から摩擦接触に至る一連の力学挙動を精

山和男：コンクリートの破壊力学に基づく等方性損

度よく再現できることが分かる．

傷モデルの定式化とその性能評価，日本計算工学会

図4 に垂直荷重を 1.0 MPa 与えた際の変形およ論文集，Vol.2013, pp.20130015, 2013.

び損傷分布を示す．界面に沿ってすべるように変[2] R. van der Pluijm: Out-of-plane bending of masonry

形しており，損傷した要素において接触を再現で

behaviour and strength, PhD thesis, Eindhoven

きることが分かる．

University of Technology, The Netherlands, 1999.

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 7高効率・高耐久性色素増感太陽電池の実現に向けた新規ポルフィリン
色素の開発と電子構造の解明
東野	
  智洋
京都大学大学院工学研究科分子工学専攻

1

緒言

は酸化チタン表面に対して非常に高い吸着能を示すことから[1]、本研究ではトリエトキシシリル基

色素増感太陽電池において、色素を酸化チタンを導入したポルフィリン色素 ZnPSi1、ZnPSi2、上に固定する吸着基はエネルギー変換効率だけでおよび ZnPSi3 を設計、合成した（Figure 1）。[2]

なく、セルの長期耐久性においても重要な役割を果たす。増感色素の多くはカルボン酸部位で酸化

2

結果と考察

チタン表面に吸着しているが、その吸着力が十分でなく、酸化チタンから脱着しやすいことが素子の耐久性を下げる一因となっているため、より強く吸着できる吸着基を用いることで長期耐久性の向上が期待できる。特にトリアルコキシシリル基2.1 新規ポルフィリン色素の色素太陽電池性能と長期耐久性評価
参照色素 YD2-o-C8 を用いてセル作製条件の最適化を行ったところ、色素増感太陽電池の変換効率は ZnPSi1、ZnPSi2、ZnPSi3 を用いたセルでそれぞれ 2.2 %、4.7%、2.3%となった。カルボン酸
吸着基をもつ YD2-o-C8 を用いたセルでの変換効率は 7.7 %であり、シリル吸着基を導入した色素はいずれも変換効率が低下することがわかった。
しかし ZnPSi2 を用いたセルでは、これまでに報告されているシリル吸着基を有するポルフィリン
色素のなかでも最も高い変換効率を達成した。まFigure 2. Plots of the normalized η value (average of
two independent experiments) as a function of
Figure 1. Molecular structures of porphyrin dyes.

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

illumination time for the DSSCs based on ZnPSi2(solid line) and YD2-o-C8 (dashed line).

た、ZnPSi2 と YD2-o-C8 を用いて光照射下における色素増感太陽電池の長期耐久性試験を行った結果、250 時間の光照射後の変換効率が、YD2-o-C8を用いたセルではセル作製直後の65%程度にまで低下したのに対し、ZnPSi2 を用いたセルでは作成
直後の 85%程度の値を保持することがわかった（Figure 2）。このことから、シリル吸着基はカルボン酸吸着基よりも酸化チタン膜への強い吸着能を有し、色素増感太陽電池の長期耐久性向上に有効な吸着基であるといえる。

2.2 理論計算
ポルフィリン色素の最安定構造およびそのフロ
ンティア軌道の電子構造について知見を得るために、Gaussian09 プログラムを用いて密度汎関数法（ DFT ） に よ る 理 論 計 算 を 行 っ た（B3LYP/6-31G(d)）
。Figure 3 にポルフィリン色素
ZnPSi1、ZnPSi2、および ZnPSi3 の最高被占軌道（HOMO）と最低空軌道（LUMO）における軌道分布を示す。一般に、アンカー部位の LUMO の軌道分布が大きいほど、励起状態の色素と酸化チタンの 3d 軌道との間の電子カップリングが大きくなることが知られている。しかしシリル基上にほとんど軌道が見られないことから、色素が酸化チタン上で傾き、空間を介した相互作用によって電子注入が起こると考えられる。そのため、セル性能の違いは分子の形状によって傾き角が変化し、電子移動の起こりやすさに差が生じたためと推察される。このように、本計算結果は実験結果の理論的解釈の一助となった点で意義がある。

3

Figure 3. Selected Kohn−Sham orbitals for (a) ZnPSi1,(b) ZnPSi2, and (c) ZnPSi3 obtained by DFT
calculations with the B3LYP/6-31G(d) level. To
simplify the calculations, hexyl groups on thediarylamino groups, octyloxy groups on the meso-aryl

参考文献

groups, and triethoxysilyl groups were replaced with

[1] (a) K. Kakiage, Y. Aoyama, T. Yano, T. Otsuka, T.
Kyomen, M. Unno, M. Hanaya, Chem. Commun. 2014,

methyl,methoxy,

andtrimethoxysilylgroups,respectively.

50, 6379–6381; (b) B. J. Brennan, M. J. Llansola
Portolés, P. a Liddell, T. A. Moore, A. L. Moore, D.
Gust,

Phys.

Chem.

Chem.

Phys.

2013,

15,

16605–16614.
[2] T. Higashino, S. Nimura, K. Sugiura, Y.
Kurumisawa, Y. Tsuji, H. Imahori, ACS Omega 2017,
2, 6958–6967.

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

気液混相撹拌操作に対する数値解析山本 卓也
東北大学大学院環境科学研究科

1 緒言
材料を混ぜ合わせるという操作は各種工業プロrpm, 400 rpm で回転させることで撹拌した。撹
拌翼の形状は 4 枚パドル翼であるが、溶融金属を最終的に対象とするため翼を分厚くした。

セスで必要な操作であり、化学プロセス、金属生

数値解析に利用した支配方程式は産プロセス、生物プロセス等の幅広い分野で研究

Navier-Stokes 式、連続式であり、気液界面の追されてきた。この「混ぜる」プロセスは材料によ

跡には代数的 Volume of Fluid(VOF)法[2]を用いって混ぜ方が異なり、対象に応じた混合方法が開た。また、表面張力モデルには CSF モデル[3]を発されてきた。

利用した。本解析条件では翼端直径基準の撹拌

溶融金属に対しても「混ぜる」操作は頻繁に行

Reynolds 数が 18000, 24000 であった。
このため、われている。
金属分野では撹拌翼を用いた「機械」

乱流モデルを利用し、Large Eddy Simulation を

撹拌と電磁効果による「電磁」撹拌の２種類が用行った。気液界面変動、そして乱流挙動を追跡す

途とスケールによって使い分けられているが、よるため Taylor スケールのおおよそ半分程度とな

り大規模な装置に対しては機械撹拌が用いられる。

るように計算格子を分割し、総格子セル数は約

ここで、溶融アルミニウムの機械撹拌操作に注目

400 万であった。また、これらすべての計算はオ

すると、溶融アルミニウム中から不純物を取り除ープンソースであるOpenFOAM を用いて計算をくために、塩素ガスを吹き込む、もしくはフッ素行った。

系粒子を混ぜることで不純物を取り除く[1]。しかしながら、アルミニウムは反応性が良いために、3 結果と考察撹拌操作を行うと気液界面上で空気中の酸素や水蒸気と反応し、酸化皮膜を形成する。この酸化皮

まず、計算結果を検証するために、気泡巻き込

膜は不純物となり、界面変動に伴って溶融アルミみに対する臨界回転速度から比較した。上記条件

ニウム中へと巻き込まれるが、その巻き込まれ挙において実験では 390 rpm において気液界面か動については明らかになっていない。

ら気泡が撹拌翼に巻き込まれる現象が観測された。

本研究では、上記の溶融アルミニウム中へと酸

一方で、数値シミュレーション結果からは 300

化物の巻き込みの基礎現象解明を目的として、機rpm では気泡が巻き込まれなかったが、400 rpm

械撹拌操作における気液界面変動を解明した。溶では気泡が巻き込まれた。これは実験による観測

融アルミニウムでは実験との比較検証が全くできと一致した。さらに、時間平均化した気液界面位ないため、今回の研究では対象は水であるとした。

置を実験と数値シミュレーションで比較したところ界面形状が良好に一致した。
これらの結果より、

2 数値解析手法

数値シミュレーション結果は妥当であると判断した。

対象とする撹拌槽は内径 192 mm、高さ 400

300rpm, 400 rpm における時間平均化した界mm の円筒容器であり、水を 330 mm まで満たし

面形状を Fig. 1 に示す。ここで、気液界面位置はた。撹拌翼を深さ 40 mm の位置に設定し、300

VOF 法に利用した液体率aが 0.5 の位置で可視化

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

を行った。
時間平均化すると軸対称の形状になり、ても、翼前面部分の上に位置する気液界面は上側

回転数の増加に伴って気液界面変形が大きくなる。

へと変形され、気液界面が撹拌翼後方へ向かって伸びている。この変形挙動は、流体中の圧力分布から説明できる。翼前面部分では流体を押し出す
ために高圧となっており、後面部分では減圧状態になっている。この高圧部分では気液界面が入れないため、翼前面部分では気液界面が上昇した。

Fig. 1 Time-averaged free surface shape at (a)300 rpm and (b) 400 rpm.
続いて、回転速度 300 rpm における時間変化する気液界面変動を Fig. 2 に示す。それぞれの画像

一方で、翼後面では圧力が下がっているため、気液界面を引き込むことができたと考えられる。最終的には翼後面で気泡として分断されることで、気泡として巻き込まれるということがわかった。

間の時間刻みは 0.03 秒である。撹拌槽の縁に沿って大きく盛り上がる部分が三箇所存在し、それらがゆっくりと回転している。この盛り上がり部分の回転速度は撹拌翼の回転速度よりも非常に遅いことが Fig. 2 より読み取れる。このような表面の盛り上がり部分が回転する現象は 400 rpm においても確認された。撹拌槽における研究では、撹拌翼から吐出される流れによって渦がゆっくりと振動するマクロ不安定性と呼ばれる現象が報告されている。今回、マクロ不安定性によって渦が大域的にゆっくり振動したことに起因して気液界面の盛り上がり部分がゆっくり回転したと想定される。

Fig. 3 Snapshots of free surface shape at 400
rpm when the aeration of the free surface
occurs: (a) Experimental observation and (b)numerical simulation.

4 結言
本研究では、撹拌操作中の気液界面変動に着目して数値シミュレーションを行った。気液界面は
撹拌槽内マクロ不安定性によって盛り上がり部が回転する。さらに、気泡を巻き込む場合は撹拌翼
後面へ向かって気液界面が伸び、そこで気泡として分断されることがわかった。

引用文献[1] M. E. Schlesinger, Aluminum recycling, Second
edition, CRC press (2014).
[2] H. G. Weller, Technical report, OpenCFD Ltd.,
TR/HGW/04 (2008).
Fig. 2 Time variation of free surface shape at[3] H. Rusche, Ph. D. Thesis of Imperial Collage of

300 rpm.

Science, Technology and Medicine, London (2002).

400 rpm での気泡が巻き込まれる瞬間の気液界
面挙動を Fig. 3 に示す。瞬時場の気液界面形状も実験と数値シミュレーションで良く一致している。

謝辞
本研究制度（若手・女性研究者奨励枠）を活

このことからも今回のシミュレーションの健全性用させて頂いたことを、この場を借りて厚く御礼が示されている。実験、数値解析結果の何れにし申し上げる。

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

iPS ĕĝºíŕ$ÒĻx4üû ÈÂÝlĐÅwħØÿĿŊįŃņ ă§āĉą ðır»£Ë x§ő

iPS ĕĝ$éĹŕ$÷x4üÄÈÂÝlĐÅw$ÍjħØ4ġĕĝŕ#¸°Ő4{&5Î¶s ĢċŒ«$[$ĥē4ã,ÈÂÏç #âĴÒĻ"ÈÂÏç$Łú$î$ĀĦ4³ŕê ŎIOUHSèd4÷2 Ě^ĵ$5Î¶s$¦¼û"4¾ 2ĕĝ#o×2 p $ /1iPS ĕĝŕ#¶÷2

2ENU@ ¥""1/1Ŗ ı" iPS ĕĝ4ŀ#Ĥķ2~Ğ· Ą

32

1 緒言
 iPS ĕĝ$÷x#ôò$ïĘĜ#/2ōęŕ 0KPSHġ ğwx3ºíŕ#A=UQ6FK2ī+ "32VĠ#Ğ·©ĕĝ%5Î¶s e÷2òæXpx1)* ªŏ#ĕĝ åþ2 #/1YV"#"2ñ#Ě4÷2ÈÂ
%Ě^ĵ¦¼û#"5Î¶s úõ2œ "2Õāĉ%Ü"đ«$ŕê# 2ÈÂÝlĐÅw4ħØ2½ç4Łúĕĝ#
25Î¶s.ªŏ#Ç
2ĕĝÍ4ďoÒĻ"ÈÂÖa#$Ùĩ4ġ


2 解析手法 èd$ÊĽÏĆ¬%ĸė¬ Navier-Stokes ÏĆ¬1ĐÅwħØ#%IOUHS$ĺwÏĆ¬4÷


dus
(1)
m
=
(FD + ( rf - rs )g)dt
rf (uf - us ) 2(2)
FD = CD Ap
2ë  f, s %33èd Đ4ģm%Đ$ıŀu %Ķ«FD %¬(2)ģ32èdÁÀρ %¡«C D%ÀshÍ%¿°ŏĈAp2ÈÂĚ Đ %k®·Ģċ2 `ĐÍ%208iĐ±%1 mmĐ èd$âľ%1.005qÔÖa Đ4čłň#ĽęħØő%ý±63 mmŖ36 mm $mĎÈ

å

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

ÂÝ#ý±17 mmVŖ$ÈÂĳ4|1^¯ò2Ě¯ò% Fig. 1(a)#ĄVáĚĲ
Ķ«%60 rpmĐ èd$âľ%1.04 WćŔ$ø"2wđ«4ÃIOUHSèd ¬(3)ģ
32É·èd4÷
æ 1.5 ´ 10 -5γ! £ 20çη = ç 3 ´ 10 -5 ´ γ! -1 20 < γ! < 60ç 0.5 ´ 10 -5
60 £ γ!
è(3) η  γ! %wđ«{'5ÎĶ«4ģŕê$5ÎĶ«#¢2wđ«x% Fig. 1(b)#Ą ÊĽÏĆ¬$ŋÌx#%ÓŅdĈç4÷Ñł#¢% Euler Ňħç4Ċł#¢%[ßĒ«Z´¨p4÷sTĶ«ĸ»% PIMPLE ç#/ġ :UKSCUA$èdħØCJH8962OpenFOAM 4;ADL7B÷


3 結果と考察 Figure2(a) #AHU<Aĵb4÷ŉ$Đ#e÷25Î¶s$Òj$Ñłx4Ą
Fig. 2(b) #ªŏ#Ç2ĐÍ$Ñłx4ĄIOUHSèd%đ« )!
5Î¶s$Òj% ªŏ#Ç2ĐÍ%¤"2p *É·èd$#%5Î¶s$ÒjªŏÇĐÍ -#Zđ«IOUHSèd/1-c¾032
p Đ$ÊĽÏĆ¬ /'èdÁÀ5Î¶s$¬4_X#Ą¶vs= FD + Fg
¶tµ A
FD = f s (v f - vs )
2 Ds
2τ = FD
3
ms

(4)

(a)

(5)
(6)

 m  v  t  µ %ıŀĶ«Ñłđ«2* FDFg A  D %èdÁÀľsâģŏĈ /'Đ±4ģ 1ë  fs %èd Đ433¹2¬(4)/1Đ
4éĹ2s%èdÁÀ$+1¬(5)/1è
dÁÀ%đ«#âf2,đ« )!ªŏÇĐÍ ì¥2*¬(5)(6)/15Î¶s%đ«#âf2,đ« )!
5Î¶s "2 řIOUHSèd%5Î¶s4¾ Đ4é  %Ō2VÏřÉ·èd4÷2 
5ÎĶ«$"ÈÂĚ1$+đ«4X2¦¼û"5Î¶s$t4ń2 ě 032*kdû#%đ« ,Đ$åþ-¾ 2
 "0řªŏÇĐÍ%w 1řĖĨŀ $>SKQ%Yzp2Ś]²%ĨďÑ
ł4./1àĂ"ĖĨŀ4³ř0#Ĭĕ"Įĭ4ġ µĥ 2

(b)



 


Fig. 1(a) Numerical configuration of one-level model.
(b) The kinematic viscosity of culture fluid.
(a)

4 結言
 iPS ĕĝ$ºíŕ$÷x(ºíŕÑ$ĕĝÅw$ÍjħØ4ġIOUHSèd{'É·èd4÷éĹŕ# 5Î¶sªŏ#Ç2ĕĝÍ4ĪgÈÂĚ^ĵ­5Î¶s2  œ%É·èd4÷2 ħä32 2Ą3]²%řĚ¯ò$ÒĻx-ġ/1vó
û"ÈÂÖa4ÆĔ2
謝辞 Õāĉ%VļČāĉŁúç\ÐÕyùāĉŁ
úÞÛŗAMEDŘ$nõyù$öÚx#ĕĝĤķTt§?AGM$ŁúKR@9<Hř /' MEXT ąāİ 15H04173$u»4}

(b)
Fig. 2. Time variations of (a) maximum shear stresses;(b) the number of particles on the bottom.

参考文献[1] Adamo, L. et al., Nature, 459, 1131-1135(2009).
[2] Yamamoto, T. et al., J. Chem. Eng. Japan., 51,
423-430 (2018).

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

パワーデバイス用半導体製造装置設計、最適化に関する数値解析
堀内 鷹之大阪大学大学院基礎工学研究科物質創成専攻化学工学領域

1 緒言

算を 3 段階に分けて行う。これらの計算についてはオープンソースである OpenFOAM[7]に導入し解シリコンカーバイド(SiC)は高い熱伝導率や広い析した[3,8]。結晶成長は種結晶面内への拡散が律速

バンドギャップを持つことから、シリコンに取って段階とし、成長速度 vg は界面の法線ベクトル n と溶代わるパワーデバイス用半導体として期待されて液内の炭素濃度 C を用いた n  C の値で評価した。

いる。シリコンカーバイドは昇華法による製造が主

結晶成長分布の均一度は以下の式(1)で表される標流であるが、過飽和を駆動力とした溶液成長法によ準偏差 γ を用いた。

る製造は転移欠陥の少ないより高品質な結晶をよ 4 v r   v 2 dr g, rep

 g
γ 0

4


0 dr



り低コストで製造することが出来るため注目されている[1]。ここで、Si 溶液への炭素溶解度が小さい
ために成長速度が極めて遅く、高品質結晶成長を維

0.5v g, rep(1)

ここで vg,rep は 0 ≦ r ≦4 mm における平均の結晶

持しつつ高速成長化させることが要求されている。

成長速度であり、変数の上付きバーは 100 秒間の時

三谷ら[2]は Si 溶液に Cr や Cr/Al を添加剤として間平均値を意味する。定義から γ の値が小さいほど用い、炭素溶解度を向上させることで高速化を図っ

均一であることが示される。

ているが、多量の添加剤が必要となり、結晶に取りまた、ParaView[9]のトレーサー粒子機能を用い込まれない制御が不可欠である。また、山本ら[3]はて 2 秒毎に仮想粒子を液中に流し、液中の全粒子数

るつぼ内に温度差のある系を数値解析した結果、種

ntot と種結晶付近に到達した粒子数 nr を用いた式(2)結晶下部流れは自由表面上の温度差に起因するで表される雑晶付着確率 ψd を求めた。

Marangoni 対流による下降流が支配的となり、結晶付与しても Marangoni 対流支配は変化しないこと

 n dtψ 
 n dtが明らかとなった。
そこで本研究では Marangoni 対

るつぼ径として溶液領域の半径 R を 8 mm から

成長速度分布が不均一化した。更に、種結晶回転を流を抑制するためにるつぼ径変化やるつぼ回転を検討した。

2 数値解析手法
解析領域は既報[3]と同一であり、概略図を Fig. 1に示す。炉内の計算には Derby らのグループによって開発された Integrate Process Model(IPM)を用いた[4-6]。この手法は、磁場は温度、Si 溶液流れに依

100

d

0
100

0

r

(2)

tot

22.5 mm まで変化させた。るつぼ回転は R = 22.5mm の系で考慮し、回転速度 ΩC = 2 s-1 まで計算した。ΩC = 2 s-1 の系については種結晶回転 ΩS = 6 s-1 を加えた系も検討した。

3 結果と考察
Figure 2 に γ、ψd に対する(a) R、(b) ΩC の影響を示す。横軸には種結晶半径 RS (= 5 mm)との比 R/RS

存しない仮定をおくことで、(i) 軸対称を仮定した電及び回転 Reynolds 数を用いて無次元化した。R/RS

磁場解析、(ii) 輻射を考慮した三次元炉内熱解析、の減少に伴って、成長速度の均一度は向上する一方(iii) 軸対称を仮定した Si 溶液内移動現象解析と計で、雑晶付着確率は 3 倍程度に増大することがわか

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

った。これより R/RS が 3 より小さい系では高品質な結晶を製造できないと考えられる。一方で、高回転
Reynolds 数では成長速度の均一度向上及び雑晶付
着確率低下となることがわかった。しかし、Fig. 3 に示す溶液内速度ベクトル分布と過飽和度分布から、
種結晶回転を加えても、結晶中心部で淀み領域が生じてしまうことが明らかとなった。従って、より液深の浅い系で結晶回転の pump-up 効果を生じやすくすることで改善がみられると考えられる。

4 結言
溶液成長法によるSiC 結晶成長時のるつぼ径及び
るつぼ回転の影響に関する数値解析を行った。結果

Figure 1 Numerical domainより小さいるつぼを用いると雑晶付着確率が高くなってしまうため高品質結晶成長に適さないことがわかった。また、るつぼ回転は全体的な高品質結晶に繋がるが、種結晶中心部に局所的に淀み領域を生じさせるため更なる改善が求められる。

引用文献[1] K. Murayama et al, J. Cryst. Growth, 468, 874878 (2017).
[2] T. Mitani et al, J. Cryst. Growth, 468, 883-888(2017)[3] T. Yamamoto et al, J. Cryst, Growth, 470, 75-88(2017)[4] P. M. Gresho and J. J. Derby, J. Cryst. Growth,85, 40-48 (1987).
[5] J. J. Derby et al, J. Cryst. Growth, 97, 792-826(1989)[6] T. Munakata and I. Tanasawa, J. Cryst.
Growth, 206, 27-36 (1999).

Figure 2 Effect of R, ΩC on γ and ψd[-][7] OpenFOAM, (see https://www.openfoam.com/)[8] T. Yamamoto et al, J. Cryst. Growth, 474, 50-54(2017)[9] ParaView, (see https://www.paraview.org/)
0.1 m/s

Figure 3 Time-averaged velocity vectors

and supersaturation degree in the case of
R = 22.5 mm, ΩC = 2 s-1 and ΩS = 6 s-1

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 15カスケード選択型分子動力学シミュレーションの開発と	 	
  タンパク質構造変化予測への適用	 原田 隆平
筑波大学計算科学研究センター

1 はじめに	  タンパク質機能発現に関係する構造遷移
は，通常の分子動力学 (molecular dynamics:
MD) シミュレーションが到達可能な時間ス
ケールと比較して遥かに長時間で観測される「レアイベント」である. 故に, MD シミュ
レーションの時間スケールの制約から, 従来
法によりレアイベントを再現・予測することが困難な場合が多い. 時間スケールの制約を打開するため, 様々計算手法が開発されている. また, レアイベント自体が確率過程であるため, 長時間 MD シミュレーションを実行
してもレアイベントを再現できる保証はな
く, 不確定性が伴う. これまで我々は, 新しい分子シミュレーション技術として, タンパ
ク質機能発現に重要な構造遷移であるレアイベントを再現する計算手法を開発した. 本
稿では, 構造探索手法の代表例である「Parallel cascade Selection MD (PaCS-MD)」1-2とその拡張型 3 について概説する.

2 PaCS-MD	  PaCS-MD は, レアイベントを再現•予測する構造探索法であり, タンパク質に限らず様々なシステムへの適用研究が展開されている. PaCS-MD は, 長時間 MD シミュレーシ
ョンに代えて初期構造の異なる短時間 MD
シミュレーションを独立•並列に実行する「構造リサンプリング」のサイクルを繰り返すことで, 効率的な構造探索を実現する. 図
1 に PaCS-MD のフローチャートを示す.
	  PaCS-MD は, 構造遷移の始構造と終構造が既知である場合に, 構造間の遷移経路を探

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

図 1: PaCS-MD のフローチャート
索することが出来る. 図 1 に示す様に, 始構
造から終構造へ遷移する可能性が高い候補
構造を反応座標に基づき選択し, 短時間 MD
シミュレーションをリスタートする構造リ
サンプリングのサイクルを繰り返すことで,構造分布が終構造の近傍に近づいていく.
様々な反応座標が考えられるが, 最も単純な例として, 終構造が分かっている場合には,終構造から測定した平均自乗距離(root-meansquare deviation: RMSD)が考えられる. 各サ
イクルで得られるトラジェクトリを RMSDの値でランキングした後, 値の小さい順に複数個選択し, 短時間 MD シミュレーションの初期構造として構造探索を繰り返していくことで, 徐々に終構造へ近づいていく.

3 拡張型 PaCS-MD	  PaCS-MD を効率的に実行するには, タン
パク質の構造遷移を記述する適切な反応座標を指定しなければならない. 従来のPaCS-MD では, 遷移経路を記述する反応座標として RMSD を選択してきた. RMSD を反応座標に選択した場合, RMSD の値が常に小さくなる方向へ構造リサンプリングが進む.
しかし, 一時的に RMSD が大きな(終構造から遠い)準安定状態を経由して終構造へ遷移する場合は, RMSD による構造変化の記述が適さない. 具体例として, タンパク質構造が局所的に一度変性し(RMSD が増大し), 再度
折り畳んで終構造へ至る構造遷移などが該当する. 故に, 適用範囲を拡大するためには,
PaCS-MD の反応座標を RMSD に決め打ちすることなく, 反応座標の候補を複数用意し,遷移経路の探索過程でより最適な反応座標を逐次的に選択していく戦略を採用する. 図
2 に, 拡張型 PaCS-MD のフローチャートを示す.具体的には, 経路探索の途中でそれぞれの反応座標方向の勾配を計算し, 一番勾配が小さい反応座標方向を最適な反応座標軸と見なし, 構造リサンプリングを継続する.
このスキームにより, PaCS-MD の構造リサンプリングの過程において, 適切な反応座標が動的に選択されていく. 勾配がフラットになってきたら, その反応座標軸方向の探索は終了したと判定され, 残りの反応座標の中で一番勾配が小さい反応座標軸に関して探索が進む. これを繰り返すことにより, より複
雑な遷移過程を効率的に抽出できる.
	  方法論の検証として, 加水分解酵素 T4 リ
ゾチームの Open/Closed 構造遷移に適用した.
T4 リゾチームの構造遷移は, 構造間のRMSD を反応座標に指定することにより容易に抽出可能であるが, 慣性半径や末端間距
離, 溶媒露出面積などを反応座標に指定した
場合, 効率的に構造遷移を抽出することが難しい. そこで, 予め反応座標を RMSD に指定することなく, 慣性半径, 末端間距離, 溶媒

図 2: 拡張型 PaCS-MD
露出面積の3つを複合的に反応座標に指定し,構造サンプリングの過程において反応座標
方向の勾配を参照にすることで動的に指定する拡張型 PaCS-MD を実行したところ, 慣
性半径, 末端間距離, 溶媒露出面積を単体で反応座標に指定した従来の PaCS-MD と比較して極めて効率的に Open/Closed 構造遷移を抽出できることが分かった.

4 まとめ
	
  本研究では, PaCS-MD1,2 における反応座標の指定方法を一般化し, より適切な反応座標
を動的に指定する拡張型PaCS-MD3 を開発し,実用性を検証した. 今後は, 複雑な構造変化を経て生体機能を発揮するタンパク質に適用し, 生体機能の解析を進める.

謝辞	 	
  本研究では, 京都大学学術情報メディアセンターの共同研究制度(若手奨励枠)を活用させて頂きました.

[参考文献]	
 [1] R. Harada, A. Kitao, J. Chem. Phys., 139,
035103 (10 pages) (2013)[2] R. Harada, A. Kitao, J. Chem. Theory
Comput., 11, 5493-5502 (2015)[3] R. Harada, Y. Shigeta, J. Chem. Theory
Comput., 14, 3332-3341 (2018)

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

粒子との接触を伴う液体挙動の直接数値解析
自由界面―固体物体間相互作用を考慮した多相直接数値解析井上 幹允東京理科大学大学院理工学研究科機械工学専攻

1

緒言

OpenFOAM を用いて水平基板上での液滴の濡れ拡がり挙動を再現し，先行研究と比較することで近年，濡れ現象を理解し応用・制御する動きが活

OpenFOAM の妥当性を検証することを目的とした．

発化している．その上で重要なものが，固体との接触を伴った液体の挙動である．この液体と固体の接

2

計算方法

触は非常に単純な現象と思えるが，液体表面の変形，固液気接触界線の移動，接触で生じる液体‐固体間の今回，解析対象として，水平基板に対して液滴を相互作用など様々な物理現象が複合したものであ落下させ，液滴を基板上に濡れ拡がらせる系を導入

るため，極めて複雑なものとなる．界面の挙動に着した．図 1 に解析系の概略図を示す．解析において目した研究は従来広く行われており，特に，水平基y 軸に関して軸対称な系を扱うことによって三次元板上での液滴挙動に関する研究は今でも活発に行を仮定した．解析系の寸法，物性に関しては表 1 にわれている．その中で，この研究の先駆けとなった

示す．支配方程式は以下の三式を適用した．

のが，Tanner(1)である．この研究は時間と液滴の濡れ拡がり半径に着目したものある．この結果，濡れ拡
がり半径を R，時間を t とした時，液滴が不揮発性で，基板に対して完全な濡れかつ表面張力が支配的な状態では，R ~ t1/10 となることを示した．これに対して，Lopez(2)は，重力が支配的な領域においては時間と濡れ拡がり半径の間に，R ~ t1/8 の関係が成り立つことを明らかにした．そして，近年では，基板に

𝝏𝝏𝝏𝝏
𝝏𝝏𝝏𝝏

= 𝟎𝟎

𝝏𝝏𝝏𝝏𝝏𝝏
𝝏𝝏𝝏𝝏

𝝏𝝏𝝏𝝏
𝝏𝝏𝝏𝝏

(1)

+ 𝛁𝛁 ∙ (𝝆𝝆𝝆𝝆𝝆𝝆) = −𝛁𝛁𝒑𝒑 + 𝛁𝛁 ∙ 𝝉𝝉 + 𝝆𝝆𝝆𝝆 + 𝑭𝑭

(2)

+ 𝛁𝛁 ∙ (𝜶𝜶𝜶𝜶) + 𝛁𝛁 ∙ (𝒖𝒖𝒓𝒓 𝜶𝜶(𝟏𝟏 − 𝜶𝜶)) = 𝟎𝟎

(3)

対して液滴を滴下した直後の挙動に注目した研究

支配方程式としては以上の三式を用いた．式(1)は連が行われている．Mitra and Mitra(3)は液滴滴下直後か続の式であり，u は各方向の速度を表している．式

ら長時間に渡って，液滴の濡れ拡がり半径と時間の(2)は運動方程式であり，は密度，p は圧力， は粘関係を観察したところ，両者の関係が R ~ t から R ~性項，g は重力，F は界面張力を表している．式(3)t1/2 なった後に R ~ t1/10 になることを明らかにした．

は VOF 法に基づいた移流方程式であり，界面の移本研究では，この液滴の濡れ拡がりに関して，水流をこの式によって考慮している．

平基板を用いない場合における液滴の濡れ拡がり挙動に着目した．実験では，濡れ性などの実験条件の厳密な設定が困難であるため，数値解析によって液滴の濡れ拡がり挙動を観ることとした．本研究では，汎用ソフトである OpenFOAM を用いて解析を行った．まず，OpenFOAM を本研究に導入することが適切かどうかを判断する必要がある．そこで，

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

Fig. 1 Schematic of the computational domain.

Table 1 Computational condition and physical property.

フ中に二種類の近似曲線を示す．近似曲線はそれぞれ R ~ t1/8, R ~ t1/10 を表している．この曲線と濡れ拡がり半径の時間変化を比較すると，概ね傾向が一致することがわかった．この傾向を先行研究である
Mitra and Mitra(3)の傾向と比較した．この先行研究では，針から液体を供給し続け，下壁面に液滴を濡れ拡がらせるという実験を行い，時間と濡れ拡がり半径の関係について議論した先行研究である．先行研式中の は格子中の液体含有率，ur は格子内の流体

究では液体を供給し続けているのに対し，本研究での相対速度を示している．この式によって，界面のは決まった量の液滴を滴下させるとい系であるた

再構築を可能にしている．

め，多少の違いはあるが，比較に際し，ここに大き3な違いはないと考えられる．比較の結果，本研究で結果と考察得られた結果はこの先行研究と良好な一致を示す
ことがわかった．

 Fig.2 は液滴の濡れ拡がりの様子を示した図であ この結果，本研究で焦点を当てている微小空間のる．図に示すように，今回は液滴が基板に接触したスケールにおいて OpenFOAM を用いた解析を行う

瞬間の時刻を t = 0 と定義した．この図より，液滴がことが適切であることを，先行研究と比較すること

基板を濡れ拡がる様子が確認できる．この結果をもによって定量的に示された．

とに，液滴の濡れ拡がり半径を追跡した．

Fig. 3 Temporal variation of the radius of spreading droplet.
Red and blue curves indicate approximate curve and green lineis the transition point.

4

結言

Fig. 2 These pictures indicate the volume ratio of the liquid at

本研究では汎用ソフトである OpenFOAM を用い

each time. Red and blue regions indicate liquid and gasて，水平基板上の液滴の濡れ拡がりの再現を試みた．

respectively.

その際，特に，時間と濡れ拡がり半径に着目した．
その結果，時間と濡れ拡がり半径の挙動に関して，Fig. 3 は時間と濡れ拡がり半径の関係を示した結果

先行研究と良好な一致を得ることができた．

であり，横軸に時間，縦軸に濡れ拡がり半径を示し

本研究では，障害物との接触を伴った液滴の濡れた両対数グラフになっている．この結果，濡れ拡が拡がりにおいて見られた特殊な現象の物理的な機り挙動が前半と後半で変わることがわかった．グラ

構の解明を目的としている．そこで，将来の展望と京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 19

しては，OpenFOAM を用いて，この現象を再現し，この特異な現象の物理的な発現機構に関して数値解析の観点から，言及していきたいと考えている．

5

参考文献

[1] Tanner, L. H., Journal of Physics D: Applied
Physics, 12. 9 (1979), pp. 1473-1484.
[2] Lopez, J., Clarence, A. M. and Eli, R., Journal of
Colloid and Interface Science, 56. 3 (1976), pp. 460-468.
[3] Mitra, S. and Mitra, S. K., Langmuir, 32. 35 (2016), pp.
8843-8848.
謝辞：本研究を進めるにあたり，京都大学学術情報
メディアセンターの方々には多大なご協力を賜りました．ここに記し，感謝の意を表します．

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

無重力下での高プラントル数流体における温度差マランゴニ効果に起因する
液柱内対流場の二次不安定性
小笠原 亨
東京理科大学大学院理工学研究科機械工学専攻

1 緒言は行われていない．したがって，本研究では二次不安定性に関して，直接数値計算および固有直交

自由界面を有する気液界面において，温度分布

分解（以下，POD）16)を用いて解析を行った．

が存在する場合，表面張力の温度依存性により表私が所属している研究室は，国際宇宙ステーシ

面張力が不均一となり，対流が生じる．この対流ョン日本実験棟「きぼう」
によって行われていた，はマランゴニ対流と呼ばれ，微小重力環境下やマ

MEIS ( Marangoni Experiment in Space ) や Dynamic

イクロスケール下で顕在化する．これは自然対流

Surf. と称される流体物理実験に参画していた．その駆動力の一つとして着目されている．例えば，して，今後，JEREMI ( Japanese European Research

単結晶材料生成方法の一つである floating-zone（以Experiment on Marangoni Instabilities ) と称される，下，FZ）法は，純度の高い単結晶材料を生成する

流体実験を行う予定である． JEREMI では，気体ことができる．一方で，材料生産の高効率化のた側に強制対流を付与し，自由界面での熱伝達制御

めに，微小重力下で利用するとマランゴニ対流のを行う予定である．

不安定性によって生じる振動流が材料生成に悪影
1)

本研究では，マランゴニ対流の不安定化メカニ

響を及ぼすことが知られている ．したがって，ズムの解明による理学的貢献，及び産業分野へのマランゴニ対流の不安定性に関する理解や制御の波及効果を目指し，研究を実施している．

観点で研究が行われている．
図１は half-zone（以下，HZ）
モデルと呼ばれる，FZ 法の半分を模擬した，最も典型的なモデルである．上下壁面間に液柱を形成し，上部を加熱，下部を冷却することで液柱自由表面上に温度勾配を付与し，上から下へマランゴニ対流を発生させる．
この HZ モデルにおけるマランゴニ対流の不安定性にプラントル数（以下，Pr）が影響を与えることが一般的に知られている 2) - 4)．まず，上下壁間の温度差が小さい時，液柱内対流場は二次元定常流になる．低 Pr 流体の場合，温度差を大きくすると，対流場は３次元振動流になり（一次不安定性），
その後三次元振動流に遷移する（二次不安定性）．高 Pr 流体においては，二次元定常流から三2) - 11)次元振動流に一気に遷移する（一次不安定性）2), 7),10)．この対流後の高次不安定性およびカオス化過程においても研究されているが 12) - 15)，高 Pr 流体
における二次不安定性の詳細な解析および定義付

Fig. 1 Schematics of half-zone liquid bridge

2 計算方法液柱の高さを H，半径を R とし，アスペクト比Γ = 𝐻𝐻 ⁄𝑅𝑅と定義する（図１）．今回の計算では，
Γ = 0.7 で固定する．液柱は動的・静的な変形を考慮せず，常にストレートな形状を保持しているものとする．計算系全体は無重力環境下を想定している．上下壁間の温度差は ∆𝑇𝑇 = 𝑇𝑇+ − 𝑇𝑇- と定義

し，表面張力 𝜎𝜎(𝑇𝑇)は温度に依存し，線型的な変京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 21

3 結果

化をするものと仮定する．
𝜎𝜎(𝑇𝑇) = 𝜎𝜎(𝑇𝑇- ) + 𝜎𝜎2 (𝑇𝑇 − 𝑇𝑇- )(1)

ここで，𝜎𝜎2 = 𝜕𝜕𝜕𝜕 ⁄𝜕𝜕𝜕𝜕．流体は非圧縮性ニュート

周方向波数（以下，m）が m = 3 で周期的な回転振動流であった対流場は Re = 3250 において準ン流体であることを仮定し，Pr が 4 であるものと周期回転流へと遷移した．図２(a) - (c)は FFT によする．無次元化された流体の運動と熱輸送を記述って得られた回転周期を元に準周期回転流を回転する支配方程式は，Navier-Stokes 方程式，連続の座標系に変換し，温度変動の等温面を可視化した式，エネルギー方程式とする．表面張力を除く物

様子を時系列に示している．等温面の形および大性の温度依存性は無視する．

きさが時間によって変形していることが確認できる．したがって，これを二次不安定性後の対流場

1𝜕𝜕𝒖𝒖∗+ (𝒖𝒖∗ ∙ ∇)𝒖𝒖∗ = −∇𝑝𝑝∗ + ∇= 𝒖𝒖∗
𝜕𝜕𝜕𝜕
Re

(2)

1 = ∗
𝜕𝜕𝑇𝑇 ∗+ (𝒖𝒖∗ ∙ ∇)𝑇𝑇 ∗ =
∇ 𝑇𝑇
𝜕𝜕𝜕𝜕
Ma

(4)

∇ ∙ 𝒖𝒖∗ = 0(3)として，以下に結果を記述する．
二次不安定性前の温度場を POD によって分解した結果，エネルギーの大きい順に成分は m = 3, 6,
9 を示した．これは元の対流場の m および高調波成分である．図２(d) - (f)は二次不安定性後の温度

ここで，u*を無次元速度ベクトル，p を無次元圧場を POD によって分解した各成分の温度変動の力，𝑇𝑇 ∗ = (𝑇𝑇 − 𝑇𝑇C )⁄∆𝑇𝑇を無次元温度とする．円柱

等温面の瞬時場を可視化した様子である．一番エ

座標系(r, q, z)を用いている．これより以下においては*を省略する．無次元数 Re は対流の強さを測るパラメーターであり，以下のように定義する．
Re =|𝜎𝜎2 |∆𝑇𝑇𝑇𝑇𝜌𝜌𝜈𝜈 =

ネルギー量の大きい成分（k = 1）は二次不安定性
前同様 m = 3 であることが確認できる．しかし，次にエネルギーの大きい成分は m = 2, 4 の形を示している．したがって，これらの成分が生じるこ(5)とで二次不安定性が発生したことがわかる．
また，

ここで，r は密度，n は動粘度，Ma をマランゴニ

この対流場からは高調波成分も得られているが，数とすると，Ma = Re Pr という関係がある．

エネルギー値が新しく生じた成分のエネルギーよ

境界条件は上下壁上では，
速度はノンスリップ，りも小さいことが確認されている．

温度は高温壁・低温壁で Th および Tc でそれぞれ一定としている．液柱自由表面上での熱の授受はないものとし，断熱の境界条件を与え，速度に関しては，表面張力とせん断応力の釣り合いから，
マランゴニ効果による流体の駆動を行っている．

(b) t = t0 + ⌧ /2(c) t = t0 + ⌧(e) k = 2(m = 2)(f) k = 3(m = 4)

計算により得られたデータから時間平均に対する温度変動分 𝑇𝑇Hを算出し，一般的に知られている
POD16)を用いて解析を行う．POD からは固有値l, 固有ベクトル a, 固有関数 f，が得られる．周期関数を POD した際，同じ固有関数が２つずつ得られることが一般的に知られているため，各固有(d) k = 1(m = 3)関数の成分の固有値を足し合わせてエネルギーEk

Fig. 2 (a)-(c) 3D visualization of thermal-fields in rotating

を以下のように定義する．

frame of reference against fundamental frequency of

𝐸𝐸J = 𝜆𝜆=J + 𝜆𝜆=JLM(6)

ここで k は正の整数であり，k が小さいほど Ek は大きくなるように定義している．また，各成分の温度変動 P𝑇𝑇J を以下のように定義する．
𝑇𝑇HJ = 𝑎𝑎=J 𝜙𝜙=J + 𝑎𝑎=JLM 𝜙𝜙=JLM(7)22 − Academic Center for Computing and Media Studies, Kyoto University

hydrothermal wave after secondary instability (d)-(f) Largeenergy components obtained by POD
図３は二次不安定性後の温度場を POD によって得られたエネルギーの大きい７つの成分の各 m
に対するパワースペクトルを示している．この三次元グラフからも二次不安定性前には確認されて4 結言

いない m = 3 および高調波以外の m のピークがPr = 4 の流体で，無重力下でのストレートな液存在していることが明確である．

柱（Γ= 0.7）を想定して，二次不安定性の遷移過程を調査するため直接数値計算を行った．二次不安定性後の対流場は準周期回転流になることを確認した．この対流場を POD によって解析し，二次不安定性前には確認されなかった m = 3 および高調波以外の成分が得られた．この成分のエネル
ギー値の Re 数に対する発達過程から臨界レイノ
ルズ数 Rec(2)を求めた．

5 参考文献1)

A. Cröll et al.: J. Cryst. Growth, 191 (1998) 365.

2)

M. Wanschura et al.: Phys. Fluids, 7 (1995) 912.

3)

I. I. Ryzhkov: Phys. Fluids, 23 (2011) 082103.

4)

K. Fujimura: J. Phys. Soc. Jan., 82 (2013) 074401.

5)

R. Rupp et al.: J. Crystal Growth, 97 (1989) 34.

6)

M. Levenstam and G. Amberg: J. Fluid Mesh., 297(1995) 357.

のエネルギー値 Ek を Re の変化とともに示す．Re

7)

J. Leypoldt et al.: J. Fluid Mech., 414 (2000) 285.

が大きくなるにしたがって，それぞれのエネルギ

8)

M. Levenstam et al.: Phys. Fluids, 13 (2001) 807.

ーは大きくなることが確認できる．Re < 3200 では

9)

N. Imaihi et al.: J. Crystal Growth, 230 (2001) 164.

m = 2, 4 の成分の存在を確認することができなかった．このエネルギー値の分布から，二次不安定

10) C. Nienhüser and H. C. Kuhlmann: J. Fluid Mech.,
458 (2002) 35.

性における臨界レイノルズ数（Rec(2)）を求める．

11) K. Motegi et al.: Phys. Fluids, 29 (2017) 074104.

臨界値近傍の３点の値を線型近似し，Ek = 0 との

12) R. Velten et al.: Phys. Fluids A, 3 (1991) 267.

交点を Rec(2)と定義した．

13) V. M. Shevtosova et al.: Phys. Review E, 68 (2003)
066311.

Fig. 3 Full spectral decompositions of surface
temperature deviation at 𝑧𝑧 = Γ⁄2 for each component
after secondary instability

図４に二次不安定性後に生じた m = 2 および 4

14) I. Ueno et al.: Phys. Fluids, 15 (2003) 408.

3015) T. Matsugase et al.: Int. J. Heat Mass Trans., 89(2015) 903.

2016) K. Li et al.: J. Crystal Growth, 307 (2007) 155.

10

0

3200

3400

Re

3600

Fig. 4 Energy of Fig.2 (e) and (f) as a growth of Re

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

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

1 背景と問題して、風鈴および弦楽器のイラスト、風鈴および弦楽器の音を、それぞれ二つ用いた。

人間の脳においては、視覚、聴覚、触覚とい感覚間協応を持つ刺激を作成するにあたり、図った異なる感覚モダリティの情報は、基本的に異形の形状と効果音の音響特徴の対応に着目した。

なる場所で、異なる方法によって処理される。しまず、視覚刺激として、鋭角と直線で構成されたかし、我々は時に、異なる感覚モダリティから得図形（直線図形）二つと、曲線で構成された図形られる情報の間に共通性を感じることがある。例（曲線図形）二つを用意した。また、聴覚刺激とえば、動物の絵や写真（視覚情報）と鳴き声（聴して、ソフトウェアシンセサイザーによって音響覚情報）はいずれも、「犬」や「猫」といった共通

特徴を操作した効果音を多数作成し、その中からの表象に結びつく。また、このような経験に基づ

上述の四つの視覚刺激のそれぞれに最も合うと感く結びつきを持たないと考えられる視聴覚刺激にじられるもの四つを、147 名の学生に対して行っ

おいても、例えば高い音には明るい色、低い音にたアンケートの結果に基づいて選択した。

は暗い色を組み合わせると、逆の組み合わせよりも良く適合して感じられるように、非恣意的な対応が存在する。この現象は様々な感覚モダリティの様々な感覚特徴の間に観察されることが分かっており、感覚間協応と呼ばれる。
近年、経験に基づく結びつきを持つ異種情報間の共通性を表現する脳のメカニズムについては、
徐々に理解が進んできている。一方で、感覚間協応を支えるメカニズムや、経験に基づく結びつきと感覚間協応の違いについては、よく分かっていない。そこで、本研究では、fMRI を用いて成人の実験参加者が様々な視聴覚刺激を観察している
際の脳活動を計測し、このデータに MVPA

図 1 実験に用いた刺激（Multivariate Pattern Analysis）解析によるデ実験に用いたイラストおよび効果音のスペクトログ

コーディングを適用することで、視覚と聴覚の感ラム。スペクトログラムは横軸が時間、縦軸が周波数、覚間協応の神経基盤や、経験に基づく結びつきと色がパワーを表す。(a)風鈴のイラストと音、(b)絃楽器

の違いについて検討を行った。

のイラストと音、(c)直線図形と適合する音、(d)曲線図形と適合する音。

2 実験方法これらの視覚刺激、および聴覚刺激を一つずつ、モニタおよびイヤフォンによって呈示し、これを男女合わせて 16 名の成人が実験に参加した。

実験参加者が視聴している際の脳活動を fMRI に経験に基づく結びつきを持つ自然物体の刺激とよって計測した。全ての刺激を一回ずつランダム

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

順に呈示する系列をブロックとし、実験全体で 10

聴覚野においては聴覚刺激について高く、かつ想

ブロックを繰り返した。スーパーコンピュータを

像課題を行った場合や、自然物体を用いた場合に用いて、fMRI から得られる生の脳活動データに特に高くなる傾向が見られた。また、多感覚野に

様々な前処理を施した上で、以下の MVPA 解析

おいても、右上側頭回では刺激と課題の交互作用を行った。

が見られ、特に感覚間協応を持つ刺激を観察のみ行う場合に、視覚モダリティ内デコーディングの3 解析方法

正答率が高くなることが分かった。
モダリティ間デコーディングの正答率は全体的

ある特定のブロック（e.g. ブロック 1）を除いに低く、ほとんどの場合においてチャンスレベル

た残りのブロックで取得した脳画像データを入力を上回らなかった。しかし、左初期聴覚野においとして、二カテゴリの刺激（e.g. 直線図形と曲線て、視覚刺激で学習、聴覚刺激でテストを行った図形）に対応する脳活動を判別するように、サポ

場合には刺激と課題の相互作用が見られ、特に自ートベクターマシンを用いてデコーダを作成した。

然物体刺激を、逆感覚モダリティの刺激を想像しそして、学習時に除外したブロックで取得した脳ながら観察する場合には有意にチャンスレベルを

画像データについて、学習した二カテゴリ間の判超える正答率となることが分かった。

別をデコーダが正しく行えるかどうかテストした。
さらに、テスト時に用いるデータを入れ替えて、

5 考察

全てのブロックで取得したデータの判別の平均正答率を求めた。なお、視聴覚情報の結びつきにつ

感覚モダリティ間デコーディングの正答率が、いて検討するという目的に照らして、視覚刺激の各モダリティの初期感覚野において全体的に高かみに対して有意に賦活する初期視覚野、聴覚刺激ったことは、知覚情報が主にそれぞれの感覚器にのみに対して有意に賦活する初期聴覚野、および対応する初期感覚野で処理されることを反映して視覚刺激と聴覚刺激のどちらに対しても有意に賦いる。また、モダリティ間デコーディングの正答活する多感覚野（上側頭回）の三領域を左右それ率が全体的に低かったことから、本研究で用いたぞれの半球について定義し、これらの領域に含ま

刺激については、
視覚または聴覚から得た刺激が、れるボクセルから得られたデータを用いて上記の反対の感覚モダリティから得た刺激と類似したパ

MVPA 解析を行った。

ターンの脳活動を励起することは基本的に困難だったと考えられる。

4 結果

一方で、自然物体については、逆モダリティの刺激を想像した場合にモダリティ間デコーディン

自然物体については風鈴か弦楽器か、感覚間協グが一部可能であったことから、経験に基づく視応を生じさせる刺激については鋭角図形（および聴覚情報の結びつきは比較的強固であり、逆モダ

鋭角図形に適合する音）か曲線図形（および曲線

リティの刺激を想像することによって、モダリテ図形に適合する音）かを判別するように学習およィ非依存の活動パターンが一部の脳領域で生じるびテストを行った。また、学習とテストに同じ感可能性が示唆された。さらに、感覚間協応刺激を

覚モダリティの刺激を用いるモダリティ内デコー

受動的に観察する場合、本研究で多感覚野としてディングと、学習とテストで別の感覚モダリティ

定義された上側頭回で視覚モダリティ内デコーデの刺激を用いるモダリティ間デコーディングの両

ィングが一部可能であったことは、感覚間協応の方を行った。これらの正答率データに一般化線形特異性を反映していると考えられる。感覚間協応

混合モデルを適用し、刺激（自然物体・感覚間協による視聴覚情報の結びつきは、経験に基づく結応）および課題（観察のみ・想像）の効果についびつきとは異なり、逆モダリティの刺激の意識的て検討した。

な想像を行うことなく当該の刺激をただ観察する

まず、モダリティ内デコーティングの正答率は全体的に、初期視覚野においては視覚刺激、初期場合に生じる、自動的なモダリティ非依存の脳活動パターンに基づいている可能性がある。

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

随伴解析を用いた物体表面形状最適化による抵抗低減ŲèāµøăĦ¢¦¢¦Ť¯¦ĜĤŤ

1 諸言Ărd!ĺ 1ôřď,Ărd3Ġ 0œ
Ŗù#ďrľŬ$Ăr#ĮÅ!-0Ö×ÎÍ ėē0Ö×ÎÍ$7KY:^Õ£!Ěĵôřď#Ė,œŖù#Ďő!¼ŭ3c
0+õ Ö×ÎÍqĆËĽ Ôð1
0ďr»Đ$Łý10iõ!- 

⎞
⎛ ∂ψ
T
ρ ⎜ − u + ( ∇u) ψu ⎟ = ∇ ⋅ σ ψ
⎠
⎝ ∂t

{

T

σ ψ = −ψ p I + µ ∇ψu + ( ∇ψu )}∇ ⋅ ψu = 0
Ψ = (ψu ,ψ p )#ţ p+Ö×ÎÍqĆ#+#ďr» ψ u ψ p $11Ř¶ u] p !

Đæś$¯¦Ę!şŁ ĜĤōů/#«0Ŧp Ü/ ρ $©¶ t $äŢ I

+#ÊĀ Ü¡ Ôð10éĜĤ$

$ħßH\DYµ $ĮÅtÜ0ğŪķ,ěÿW6?^Ī#- _õĂd!ĺ 1Ŧpßġ¸$z# Navier-Stokes ßġ¸|îĐďr#ľâ»Đæś3«ŏÜõ ĥŢÂ]äŢÂ#»¸ľ10+

u Ă r ņ ì (Computational Fluid Dynamics,

#- vÂßġ¸İ3R;JY]H\DY

CFD)Ŧpņì(Adjoint Method)!-Ö×»¸#ĐÈņ   0 OpenFOAM[2]ÎÍqĆí æ¢ 0»Đ3Ń033Ĕņì3ļOpenFOAM $ĥŢũęĘ0Û3çţrģĀŘ¶-#8GPX\<!

2 解析手法-ßġ¸3ņ0Ă1$eĂĐÈ

SIMPLE Ā& PISO Ā3ĔĢťņĀ!
0+ňĩï¥-/¬B=^Y(Sub-gridņì«ŏ$§ŰŪķ(s%[1] )%1

Scale, SGS) # Ĉ 3 V I Y    ň ĩ  0

0²Ő#ŪķĿł»Đ3÷Ø|Ĩ»Đ/

Large-eddy Simulation (LES)3Ĕņì3ļ

Űħ#¹Űëmb#ŰÎÍ3æ¬0

SGS VIY!$öĊBU>X\B9^VIY§ŰŪķ#Ŀł»Đ$¡Ņ».

3Ĕňĩï¥#JT[A^$.4|î4ċ/7\SBħ Ôð13/* O ï¥(BM\ß!$ħŢŧ#Ï

0 ò|3Ò 0- .4»Đ3

~ï¥)Ă{xĕ!$_õĂëmĂ«ŏÞŬò|»ĐŬģ3§.4#~xĕ!$«ĂĂ~ëm3Ĕ_õĂ#¢NGF3 

ÞŬò|µĚ¾þ+.10Z6LYCŗ´# CFD ņì#Ùŝßġ¸$ūĹÅĂrÜ Re $ 3900 ($ 10000 Ărņì#Ùŝßġ¸0 Navier-Stokes ßġ¸&Ś$Ĳĭ |î»Đ 2 ŗ/#NGF»Đ3Ñ.

Ķ#¸0_ßŦpņì$Ŧp Ü Ψ

4|î# CFD ņì3$+!ļ#¿ĉ' lb#ßġ¸ z# Navier-Stokes

Ŧpņì3ļŦpņì#~$ďßġ¸&ŚĶ#¸ .ý+.10[3]

rľŬa#Ŧp Ü#u³ c.113Ö×ÎÍ!«0Ç¶»Đ å(éņì26 − Academic Center for Computing and Media Studies, Kyoto University

$.4#NGF)!Ĕ

参考文献

3 計算結果及び考察[1] ĄÚã, dé, dêw, ª®·g, ćŜąī, “qŰŪķ#šė”, Ūü¦oŎÝŋŷ, 120($+!ņìÊĀ#¤ºÅñŊ¨

± 12 , pp.1723-1729 (2000).

ű,j#ÜuņĀ!-0ĵí {ÊĻ qZ6[2] OpenFOAM® -The open source CFDLYCÜ Re=3900 #! CFD ņì3ļ

toolbox-, https://www.openfoam.com/

á½#ĜĤĵí#ûŒ3ļňĩĵí[3]-/ďr!y ÍtÜ3ý+0¨ű,j#formulateon for the computation of topologicalňĩó"u Á.1OpenFOAM !-0

and surface sensitivities of ducted flows”, Int. J.

ņì ¤º0 ĝŌ 

Numer. Meth. Fluids, 58, pp.861-877 (2008).

C.

Othmer,“A

continuous

adjointŘ¶I^E .Ĉ¶,QX@H53ĔeĂd#Ĉôř3Ð~ń3ļ 1 !
ğ-!Ĳĭ |î#¿ßŮėē
0`úzĘ ĸĈų[^YĈŴ .4!
å0!-Ć­($ÉŠ#Ì Ń.1ľŬ»Đ# å!-Ö×ÎÍ#ÔŇ(a) ďr»Đ¨Ũ!ļ0 ĝŌ 
Ŧpņì!-R^B0.4»Đ!
«Ç¶ì0aŕ#Ĉôř#ń!
Ĉ#ėēčkŔ! Ö×ÎÍ!«0Ç¶ ¢ /Ŧpņì!-eĂ

Ĉ#Ð~!-0ďĒĘņŞ3Ĕæś»Đ#÷
ĳõ Ç¶ņì3ļ0 2 #Ç¶ņì!ÀÇ¶ ¬(b) Ă1ßĈ¶

0ß!.

4#NGF3 å}¶ CFD ņì3¨à0-/ÍtÜ ¬ »Đ3Óĳ 0
 2 

4 結言(c) BM\ßĈ¶éĜĤōů3¨à0!-/ŗ´#ĂrņìŦpņì3Ĵ)20!-/|îĐďr!y Ö×ÎÍ3đĘ!æ¬0ÜuĘņìÊĀ#ĞĘ ñŉ3ļ  éņì$.4#NGF3 å0#ëmå ¨Ũ#ďr»Đ]ľŬ»Đ$nÆ!`úzĘ! åĻ/Ç¶ņì!À(d) QX@H5^©¶»Đ3Č. ! å0V^O5\<3Ĵ){10ÃŁ 0h¿$đĘ]Ųį¶ V^
O5\<ÊĀ3¨ŀ¡õ ôřď»Đ]ıŵ `úzĈôř#ńŶ(°)|î().ëm!«ÄĻ æśÊĀ3ôĬ0f§4

0

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

アンブレラサンプリングを利用した自由エネルギー反応経路探索法の開発
満田祐樹
大阪大学理学研究科化学専攻

1 背景

2 手法

自由エネルギー反応経路は反応物や生成物の構

自由エネルギー上の反応経路については、Bother

造、自由エネルギー差、自由エネルギー障壁といっらによって、アンブレラ積分法[3]を利用した多次元

た化学反応の性質を明らかにするために重要な指PMF 上の安定点および鞍点の、構造最適化および標である。反応座標に沿った自由エネルギー変化で最小自由エネルギー経路の計算をするアルゴリズ

ある平均力ポテンシャル(PMF)を分子動力学計算にムを開発されている[5]。しかし、この手法では鞍点よって求めることは、計算化学において重要な役割

付近の初期構造が必要であり反応経路の自動探索の一つである。PMF を計算する方法は様々だが、代は実現していなかった。その一方でポテンシャル面

表的な手法として、アンブレラサンプリング法

上における反応経路の自動探索は超球面探索法[1],[2]がある。この手法は、人為的なバイアスポテン[6],[7]によって実現されている。この方法では、ポテ

シャルを分子動力学計算に加えることで、指定したンシャル面が安定点もしくは解離チャネルがある

範囲でのサンプリングを可能とする方法である。こ方向に、非調和的にゆがむという点を利用する。そのバイアスポテンシャルをかけたサンプリングをこで基準振動数でスケールした超球面上で非調和

ウィンドウと呼ぶ。ウィンドウ 1 つは狭い範囲しか下方歪みが大きな方向に探索することで、鞍点を発

サンプリングされないが、バイアスポテンシャルを見し、反応経路ネットワークを自動で明らかにする

変えながらウィンドウを増やしていき、それぞれのことができる。この非調和的な下方歪みは PMF 上

サンプリングが重なり合うようにすることによっでも同じように起こっていると考えることができて、広範囲の PMF を計算することができる。

る。そこで我々はこの 2 つの手法を組み合わせるこ

我々はこれまでの研究で、アンブレラ積分法[3]をとによって、PMF 上で自動探索ができる手法を開利用することで、自動でパラメータを決定するアル発した。FERRMap 法では、初期構造として安定点

ゴリズムの開発を行った[4]。これによって、目的のを使用する。そこから、超球面探索法によって、周範囲に対してアンブレラサンプリングによる PMF

囲の鞍点を発見する。発見した鞍点からアンブレラ計算ができるようになった。しかしながら、多次元

積分法を使用して反応経路に沿ったサンプリング

空間の PMF において、どの範囲をサンプリングすを行い、新しい安定点へとつなぐ反応経路を発見す

れば自由エネルギー反応経路が発見出来るのかとる。見つかった安定点から超球面探索をする、といいう問題は残ったままだった。

う形で計算を繰り返すことで、自動的に反応経路ネ

そこで本研究では、アンブレラ積分法を利用し自ットワークを解き明かすことができる。反応経路に

由 エ ネ ル ギ ー 反 応 経 路 探 索 法 (Free Energy

沿ったアンブレラサンプリングによって自由エネ

Reaction Root Mapping; FERRMap)を開発した。

ルギー変化が分かるのはもちろんのこと、各鞍点、

これによって多次元自由エネルギー空間上で反応安定点のウィンドウを解析することによって、それ経路を全自動探索できるようになった。

ぞれの構造を解き明かすことも可能である。

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

3 結果と考察

4 まとめ

本研究では FERRMap 法の有用性を示すために、我々の開発した FERRMap 法によって、自由エ

水中のアラニントリペプチドの 4 次元 PMF につい

ネルギー上で反応経路を自動探索することが可能て、このアルゴリズムを適用し、反応経路ネットワとなった。FERRMap 法は以下の点において他の手ークを計算した。計算をする際には、βシート構造法に比べて有利である。

を初期構造として計算を開始し、63 点の安定点と

各案定点から、自動で鞍点を発見し、新しい

133 点の鞍点を発見した。図 1 にその反応経路ネッ

安定点を発見できる。このとき、事前知識なトワークを示している。我々の方法によって、このしに思いがけない反応経路を解き明かすこ

ネットワークの各点における構造、反応経路の自由とも可能である。

エネルギー差が解明される。これらを解析すること

アンブレラサンプリング法を使用するため、でアラニントリペプチドの結果からアミド鎖の代他の手法に比べて、MD 計算の安定性が高い。

表的な構造と反応経路を解き明かすことができた

例えば PMF 上に分子が解離するなどの原因

（図 2）
。

によって MD 計算が不安定となる範囲があったとする。メタダイナミクス法などではそこをサンプリングしてしまい MD 計算が止まってしまうが、我々の手法ではその範囲以外をサンプリングすることができる。


FERRMap 法による計算結果は、安定点及び鞍点の構造、自由エネルギー差、反応経路が同時に分かるため、反応機構の解明に有効である。



FERRMap 法はアンブレラサンプリング法を使用するため、並列効率が高く、計算リソースを有効に利用しながら反応機構解明が可能である。

5 参考文献[1] G. M. Torrie et al., Chem. Phys. Lett. 28.4,
図 1 アラニントリペプチドの初期構造と、反応経路ネットワーク。四角の点が安定点、円の点が鞍

578-581 (1974).
[2] G. M. Torrie et al., J. Com. Phys. 23.2, 187199 (1977).

点を意味する。

[3] Y. Mitsuta, T. Kawakami, M. Okumura, andS. Yamanaka, Int, J. Mol. Sci., 19, 937, (2018)[4] J. Kästner et al., The Journal of chemical
physics 123.14 144104 (2005).
[5] M. U. Bohner et al., J. Chem. Phys. 137.3, 0–
6 (2012)[6] K. Ohno et al., Chem. Phys. Lett. 384.4–6,
277–282 (2004).
図 2 安定点のラマチャンドランプロットと代表的な安定構造。

[7] S. Maeda et al., J. Phys. Chem. A, 109.25,
5742–5753 (2005).

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

タンパク質フォールディング駆動力の確率的解析
山守 優 ∗∗

1

産業技術総合研究所 · 人工知能センター · インテリジェンスバイオインフォマティクスチーム

はじめに

アンフォールド構造 Xu に対して安定であるとは,

球状タンパク質は, 室温で溶液中という条件で特Eintra (Xf ) + ∆µ(Xf ) =定のコンパクトで安定な構造 (フォールド構造) をと

Eintra (Xu ) + ∆µ(Xu )る。フォールド構造の安定性は, タンパク質内部での相互作用および環境 (溶媒分子) とタンパク質の相互作用のバランスによって実現していると考えることができる。本研究では分子動力学シミュレーションと(自由) エネルギー解析を用いて, タンパク質フォール
ディングの駆動力を同定することを試みた。 [1] そのために, フォールド構造のアンサンブルとアンフォールド構造のアンサンブルを (シミュレーションによって) 用意し, フォールディング駆動力の候補因子について確率的解析を行った。
フォールディング駆動力としては, タンパク質内およびタンパク質ー溶媒間の相互作用成分, 即ち静電相互作用成分 ·van der Waals 相互作用 (分散力) 成
分 · 排除体積効果成分が検討対象となることが多く本研究でも上記 3 成分をのうちどの成分がフォールド構造の安定化に有用かを検討した。排除体積効果とは, 溶質の挿入に伴う溶媒分子の除去に起因する
自由エネルギーのペナルティで, 疎水性相互作用の主要な要素と考えられている。
本研究では, 室温でサンプルしたフォールド構造のアンサンブルと高温でサンプルしたアンフォール
ド構造のアンサンブルに対して, 以下のようなエネ
ルギー論的解析を行った。タンパク質のある構造をX とし, X の純水中での出現確率を P wat とすると
−kB T log P

wat

(2)

が成立しているということである。そして, 静電相
互作用成分 ·van der Waals 相互作用成分 · 排除体積効果の 3 成分は, Eintra (X) + ∆µ(X) の項に含まれている。ここで, 例えば, V elec (X) をタンパク質の分子内エネルギーの静電相互作用成分とする。このとき静電相互作用成分がフォールディング駆動力であるとは, ほとんど全ての Xf と Xu に対して

V elec (Xf ) < V elec (Xu )(3)が成立するということである。確率的な解析を用いて記述すると, Vfelec (= Vfelec (X)) と Vuelec (= Vuelec (X))

を確率変数として P (Vfelec ) と P (Vuelec ) をその各々の確率として∫pelec =

Vfelec <Vuelec

dVfelec dVuelec P (Vfelec )P (Vuelec ) (4)としたときに, pelec ≃ 1 となることである。同様の

特定の相互作用成分がフォールディング駆動力であるための条件は, van der Waals 相互作用成分と排除

体積効果成分についても成立する。ただし, 全て式( 2) の成立するもとでの議論である。
本研究では, 式 ( 2) の成立を確認するために, 水和自由エネルギーの計算が重要となる。自由エネル
ギー摂動法や熱力学積分法などの標準的な手法を用いた場合, タンパク質そのものの溶媒和自由エネル

(X) =
(1)

ギーを計算するには, 莫大な計算資源を必要とする。

とかける。ただし, kB は Boltzmann 定数, T は温を必要とするためである。本研究ではエネルギー表

Eintra (X) + ∆µ(X) + constant,

度, Eintra はタンパク質の分子内エネルギー, ∆µ は

水和自由エネルギーを表す。フォールド構造 Xf が30 − Academic Center for Computing and Media Studies, Kyoto University

計算の際に溶媒系と溶質系をつなぐ多数の中間状態
示溶液理論を用いて, 溶媒和自由エネルギーの計算を行った。この方法は, 溶媒系と溶液系のシミュレー

ションのみで, 溶質ー溶媒相互作用エネルギーの分布関数に基づいて近似汎関数を用いて溶媒和自由エ
ネルギーの評価を行うために高速な近似計算手法で,

表 1: 各相互作用成分について式 ( 4) の p を計算した値

精度が高く応用例も多い。 [2]2手法

本研究では, 3 種のタンパク質, 即ち Trp-cage, protein G, ubiqutin を計算対象とした。各々の PDB
code は、1L2Y, 1GB1, 1UBQ で, 20 残基, 56 残基, 76 残基の球状タンパク質である。水モデルには TIP3P を, タンパク質と共溶媒の力場には AM-

BER99SB を採用した。分子動力学シミュレーショ
ンは, GROMACS 5.0.7 を用いて, 定温 · 定圧 N P T
アンサンブルで行った。全ての系において, 総溶媒分
子数は 20,000 とした。本稿の研究では 2 段階の MDを行った。即ち, タンパク質の構造に自由度を与え,フォールド構造のアンサンブルとアンフォールド構造のアンサンブルを採取するためのシミュレーショ
ンとタンパク質の構造を固定し, 水和自由エネルギーを計算するためのシミュレーションである。フォー

pelec

pvdW

pexcl

Trp-cageprotein G

0.68
0.20

0.90
0.99

0.96
1.00

ubiquitin

0.46

0.96

0.98

作用 (分散力) 成分 · 排除体積効果の 3 成分に対して行った。その結果を表 1 に示す。

本研究の結果は, タンパク質フォールディングの駆
動力 (フォールド構造のアンサンプルを室温で安定化する成分) は, Trp-cage, protein G, ubiqution の全てにおいて, van der Waals 相互作用成分と排除体積効果であることを示唆した。
実際には, 本研究のような確率的解析では, 因果関係を直接議論することはできないが, さまざまな相互作用が複雑に関与する系で現象の原因となる因子をこのような方法で抽出できることが分かったのは,有用であると考えている。

ルド構造のアンサンプルのためのシミュレーション
は 300 K でアンフォールド構造のサンプルのためのシミュレーションは, 400 K で行った。水和自由エ
ネルギーを計算するためのシミュレーションは全て

300 K で行った。圧力は, 1bar に設定した。水和自由エネルギーの計算には, ERMOD 3.1 を用いた。

3

参考文献[1] Y. Tokunaga, Y. Yamamori, N. Matubayasi, Jour.
Chem. Phys., 148, 125101, (2018).
[2] Y. Yamamori, N. Matubayasi, Jour. Chem. Phys.,146, 225103, (2018).

結果および考察
フォールド構造のアンサンブルとアンフォールド構造のアンサンブルの構造的特徴は溶媒接触表面積
によって確認した。全てのタンパク質について, アン
フォールド構造は比較して大きな溶媒接触表面積を持ち, 伸びた構造を多分に含んでいた。次に, 式 ( 2)の成立の確認を試みた。そのために, フォールド構造のアンサンブル Xf とアンフォールド構造のアン

サンブル Xu において, Eintra (X) と ∆µ(X) の相関をみた。式 ( 1) のように, Eintra (X) と ∆µ(X) は傾
き −1 で強く逆相関する。その結果, ２つのアンサ

ンブル中の全ての構造について式 ( 2) の成立を確認した。最後に以上の結果を踏まえ, 式 ( 4) のように,確率的解析を静電相互作用成分 ·van der Waals 相互京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 31

津波による底泥巻き上げ量の予測と海洋環境変化に関する数値解析村山 雅子 1， 中田 聡史 2， 林 美鶴 2， 小林 英一 2， 越村 俊一 3
1 富山高等専門学校国際ビジネス学科， 2 神戸大学大学院海事科学研究科， 3 東北大学大学院災害科学国際研究所

1 緒言
南海トラフ巨大地震が今後 30 年以内に高い確についてシミュレーション計算を行い検討した．

2 粒子追跡シミュレーション率で発生するとされており，この地震に伴い発生する津波流は，日本沿岸の広い範囲に来襲し海域により強い流れとなる．津波流によって海底で底質が巻き上げられ，その後に海中において底質の大規模な移動により，底質中の有害物質や有害プ
ランクトンが海水中に放出され，海洋環境に大きな影響を及ぼすことが考えられる．さらに海水中に放出された有害物質が潮流と津波流により近隣の海域に移動することにより，海水中の有害物質の濃度が変化することが考えられる．
そこで本研究では，津波流による底泥の巻き上げを予測計算し，海中の有機物の移動量と密度の変化について数値解析し，その解析結果から津波が海洋環境に及ぼす影響について示すことにより海洋環境の保全と減災策に貢献できる要素となることを目的とした．

2.1 津波流計算
津波流は，東北大学で開発された津波モデルを用い，中央防災会議が 2012 年に公表した M9 規
模の南海トラフの巨大地震を想定した検討結果から，新想定の東南海・南海地震の津波波源モデルに基づき算出した 1),2)．
格子間隔は 1350m，450m，150m，50m の組合せを用い，Fig.1 に示す対象海域においては 50m を採用した．海底摩擦係数の値はマニングの粗度係数に基づき，津波挙動は水平方向の流速分布が水深方向に一様とした非線形波理論により計算を行った．また，陸上への遡上計算は行わずに全反射とし，計算時間は地震発生
後から 10 時間とした．

2.2 粒子追跡シミュレーション

本研究の検討対象海域は大阪湾とした．大阪湾津波による底質の巻き上げが発生した際，粒子は後背地に工業地域を有している一方，多くの魚がどのように巻き上げられ，移流・拡散されるか場があり海洋環境保全が重要な課題である．大阪を粒子追跡シミュレーションプログラム 3),4)を用湾の一部の海域においては底質中に重金属が多くいて計算した．計算期間は 3 か月間とし，季節に含まれており，津波により底質の巻き上げが発生

よる移流量の変化を検討するため，計算開始日をすると重金属汚染により海洋環境が悪化する可能①2015 年 4 月 1 日，②5 月 1 日，③6 月 1 日，④性が懸念される．

7 月 1 日の 4 ケース設定した．

本研究では大阪湾を対象に，津波流による底泥の巻き上げ量を予測計算し，海中の有機物の移動

3 計算結果量と密度の変化について数値解析した．
計算には沿岸部を高解像度かつ三次元的に表現

粒子追跡シミュレーションを 4 ケースの計算期

可能な高精度陸海統合モデルを用いた．また，海間について実施した．計算結果から得られた粒子

洋環境は季節により変化することから，津波流に水平分布の一例を Fig.1 に示す．Fig.1 は 4 月 1

よる底泥の巻き上げ量および移流輸送の季節変動

日から 3 ヶ月間の計算期間で，4 月 1 日 0 時に地

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

震が発生し，津波流により巻き上げた粒子が拡散

算開始日から 3 ヶ月間の粒子の移流輸送についてし水流によって移送される時の計算開始から 8 週

シミュレーション計算を行った．

間後の 5 月 27 日 0 時の粒子の分布図を示す．

粒子数の変化を見ると，計算開始から 4 週間後

また，水中の粒子数の時刻変化を Fig.2 に示す．

には変化が緩やかになり，8 週間後には停滞する．

Fig.2 は計算期間ごとの大阪湾内の粒子数の時刻

計算期間で比較すると，計算開始から 2 週間後か歴を示す．

ら 5 週間後にかけて，計算期間により粒子数の変化にばらつきが大きいことが分かった(Fig.2)．今
回は沿岸海洋の成層構造が発達する春夏季を対象としたが，今後の課題として秋冬季を含めた季節変動について検討を行うことが必要である．

5 謝辞
本研究における粒子追跡シミュレーション計算
は平成 29 年度京都大学スーパーコンピュータ共同研究制度の支援により行われたものです．
また，
本研究は公益財団法人住友電工グループ社会貢献基金による研究助成の一環として実施しました．
津波シミュレーションデータの整備は神戸大学大Fig. 1 Snapshot of horizontal distributions of simulated particles

学院海事科学研究科の米田将太氏および谷口裕樹

during the period from initial on 1 April to the end of the氏にご協力いただきました．関係各位に心より感

simulation on 30 June, 2015.

謝申し上げます．

6 参 考 文 献
[1] 小林英一, 越村俊一, 久保雅義：津波による船舶の漂流に関する基礎研究, 関西造船協会論文集, 第 243 号, pp.49-56, 2005.
[2] 谷口裕樹, 小林英一, 二村泰憲, 越村俊一,
村山雅子：備讃瀬戸海域における新想定モデ
ル津波シミュレーション, 日本船舶海洋工Fig. 2 Time history of the number of the simulated surviving
particles in each computation period in Osaka Bay.

学会講演会論文集第 16 号, pp.37-38, 2013.
[3] Nakada, S., S. Suzuki, M. Hayashi, S.
Koshimura, E. Kobayashi : Simulation ofheavy metal transport induced by a giant
tsunami

4 考察とまとめ

based

on

Nankai-TroughEarthquake: Application to Osaka Bay,

大阪湾を対象に，津波流による底泥の巻き上げ量を予測計算し，海中の有機物の移動量と密度の変化について数値解析した．また，海洋環境は季節により変化することから，計算開始日を 4 月 1
日，5 月 1 日，6 月 1 日，7 月 1 日の 4 パターン
設定し，津波流による底泥の巻き上げ量および計

Proc. of the 26th Int. Ocean and Polar Eng.
Conference, ISOPE-I-16-385, 2016.
[4] 中田聡史, 林美鶴, 鈴木綜人, 越村俊一, 小林英一：底泥巻き上げ・粒子追跡シミュレー
ションから推定した大阪湾沿岸水環境へ及ぼす巨大津波のインパクト，第 25 回海洋工学シンポジウム, pp.530-535, 2015.

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

分子動力学シミュレーションを用いた小分子・高分子混合溶液の相溶性の計算手法の確立
山田 一雄 ∗∗

1

大阪大学基礎工学研究科

はじめに
物性の異なるポリマー同士を混合させることで望ましい物性・機能を発現する高分子材料を開発する
技術をポリマーアロイ技術と呼ぶ。ポリマー混合物の物性は、混ぜ合わせるポリマーのお互いへの溶けやすさである相溶性に依存することが知られている。
新規なポリマーアロイ材料を設計する安価な手段として、計算機実験を用いた相溶性の解析を行う手法が期待されている。本課題の目的は、複数の高分子

ネルギーを計算することを考える。高分子に属するモノマー単位の溶媒和自由エネルギーは incremental

chemical potential[1, 2, 3] と呼ばれ、同じ構造の繰
り返しモノマーの数が多い長鎖ポリマーほど一定になることが期待される。ポリマー鎖全体の自由エネ
ルギーは incremental chemical potential の総和となるため、ごく少数の incremental chemical potential
計算から全体の化学ポテンシャルを見積もることができれば 計算時間を大幅に短縮することができる。

種からなるポリマー混合系の相溶性を、計算機実験を用いて原子レベルの分子間相互作用から決定する

4ことである。

計算手法
本研究では、全原子 MD シミュレーションを用いてポリマー鎖の incremental chemical potential を計算

2

問題点
多成分系の相溶性は、混合前後の自由エネルギー変化である混合ギブス自由エネルギーで決まる。本
研究では、ポリマー混合系における混合ギブス自由エネルギー評価を用いて、ポリマー溶液系の相溶性した。MD 計算にはオープンソースソフトウェアである OpenMM を、自由エネルギー計算には ERmod(エネルギー 表示理論による溶液系の近似理論 [4, 5, 6,

7, 8, 9]) と、pymbar (multistate BAR 法) の 2 つを用いて、自由エネルギー計算手法の比較を行った。

を見積もる。分子間相互作用を考慮するために全原子分子動力 学 (MD) シミュレーションを行い、各ポリマー鎖の平 衡化学ポテンシャル (自由エネルギー)
計算から混合ギブス自由エネルギーを見積もる。しかしながら、非常に大きな内部自由度を持ち排除体積が大きいポリマー鎖は、自由エネルギー計算に膨
大な数のサンプリングが必要となり、大きな計算コストがかかる。

5

計算結果
本課題では、polyethylene (PE), polymethyl methacry-late (PMMA), polyvinyl alcohol (PVA), polystyrene(PS) について系の作成および数値計算を行った。incremental chemical potential の鎖長依存性を確かめるため、10 ユニット程度の短鎖から 100 ユニット程度の長鎖まで、異なる鎖長でそれぞれ平衡化および平衡サンプリング計算を行った。

3

解決方法
この問題を解決するために、高分子がモノマーの繰り返し構造からなる点に注目する。高分子全体ではなく、その構成要素であるモノマーの溶媒和自由エ

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

参考文献

-10
∆Eincr(i) (kcal/mol)

-10.5
-11

[1] Sanat K. Kumar, Igal Szleifer, and Athana--11.5

sios Z. Panagiotopoulos. Determination of thechemical potentials of polymeric systems from

-12
-12.5

Monte Carlo simulations.
66(22):2935–2938, 1991.

-13
-13.5
-14

Phys. Rev. Lett.,[2] Sanat K. Kumar. The chemical potentials of

-14.5
0

10

20

30

40

50

monomer index ipolymer systems from computer simulations.
Fluid Phase Equilib., 83(C):333–340, 1993.
[3] Christopher J. Rasmussen,

図 1: ∆E v.s. index.

Aleksey Vish-nyakov, and Alexander V. Neimark. Calculationの溶媒環境との平均相互作用エネルギーの計算結果

of chemical potentials of chain molecules by theincremental gauge cell method. J. Chem. Phys.,である。横軸は高分子に沿ったモノマーの番号であ135(21), 2011.

図 1 は PMMA, 鎖長 50 の i 番目のモノマーと周りり、i=0 および 50 が高分子の両端を表している。図2 より、PMMA の場合、鎖長 50 程度であれば、端から遠い中間領域付近のモノマーに関して、周りの溶媒分子との平均相互作用エネルギーが i によらず一定の値に収束する性質があることが分かる。図 2 は鎖長6081, 2000.
[5] Nobuyuki Matubayasi and Masaru Nakahara.
Theory of solutions in the energy representa-

N=50

tion. ii. functional for the chemical potential. J.
Chem. Phys., 117:3605–3616, 2002.

-4.8∆µincr(i) (kcal/mol)[4] Nobuyuki Matubayasi and Masaru Nakahara.
Theory of solutions in the energy representation. i. formulation. J. Chem. Phys., 113:6070–

-5
-5.2

[6] N. Matubayasi and M. Nakahara. Erratum:“theory of solutions in the energy representa-

-5.4
-5.6

tion. ii. functional for the chemical potential”[J. Chem. Phys. 117, 3605 (2002)]. J. Chem.

-5.8
-6
-6.2
0

10

20

30

40

50

monomer index i

Phys., 118:2446, 2003.
[7] Nobuyuki Matubayasi and Masaru Nakahara.
Theory of solutions in the energy representa-

図 2: ∆µ v.s. index.

tion. iii. treatment of the molecular ﬂexibility.

50 の PMMA の incremental chemical potential の計算結果である。横軸は図 1 と同じく高分子に沿った

J. Chem. Phys., 119(18):9686–9702, 2003.
[8] Nobuyuki Matubayasi, Kuo Kan Liang, and

モノマーの番号を表す。図 2 より、PMMA の場合、

Masaru Nakahara.

鎖長 50 程度であれば、端から遠い中間領域付近のモ

solubilization in micelle.
124:154908, 2006.

ノマーに関して、incremental chemical potential がFree-energy analysis of
J. Chem. Phys.,i によらず一定の値に収束する性質があることが分か[9] Nobuyuki Matubayasi, Wataru Shinoda, andる。前述の通り, ポリマー鎖全体の自由エネルギーはMasaru Nakahara. Free-energy analysis of theincremental chemical potential の総和であるため、この性質を利用して計算を効率化することが可能と

molecular binding into lipid membrane with themethod of energy representation. J. Chem.

なることが示唆された。

Phys., 128:195107, 2008.

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

撹拌時における混合評価に対する数値シミュレーション
加藤 賢也東北大学大学院環境科学研究科

1 緒言

2 実験方法、数値解析手法

アルミニウム製造プロセスにおいては、アルミ実験装置概略図を Fig. 1 に示す。一部に傾斜がニウムのインゴットを溶解炉において溶解し、溶ついた矩形容器に斜めから撹拌翼を導入する。撹

解炉、もしくは保持炉で成分調整や溶湯清浄化操

拌翼を十分長い時間回転させ流れを発達させた後作を行う。この成分調整や溶湯清浄化操作においに、KCl 溶液をトレーサーとして図の位置に導入ては、大量のアルミニウム溶湯を高速に処理するする。KCl 溶液の濃度によって電気伝導率が変化

必要性があるため、機械的に撹拌される。一般にするため、電気伝導率の時間変化をそれぞれ A~D

アルミニウム産業界で用いられる溶解炉の容量は点で計測する。その電気伝導率が最終的に落ち着規模にもよるが数トンから百トン程度のものであいた値の前後 5%以内に収束した時間を均一混合る。さらに、アルミニウム溶湯は８００oC 程度、時間であると定義し、混合性能を評価した。

かつ不透明である。このため、大規模な装置内の流動を把握するのは非常に難しく、経験的に撹拌性能が評価されている。
一方で、上記のプロセスの処理効率を向上させると、高温保持する燃料費や排出二酸化炭素を大幅に低減することができるため、プロセスの効率を向上することは非常に重要な課題である。そこ
で、アルミニウム溶湯の撹拌操作は数値解析や水を用いた実験と組み合わせて評価されてきた。
Bujalski ら[1]は水を用いたスケールダウンモデ
ルで流動状態を評価した。一方で、Bilodeau ら[2]は実機炉での不純物の物質移動と数値解析結果を評価した。他にも多種の形状、パラメータにおいて研究が行われているものの、実機炉中での物質移動と数値解析結果が良好に一致しているものはほとんどなく、炉内物質移動を正確に表現できているものはほとんどない。
本研究では、機械撹拌操作を用いた際の炉内物質移動を解明するため、水を用いたスケールダウ
ン実験とその場合に対応する数値シミュレーショ
ンを行い、スケールダウン水モデルにおいて炉内物質移動を解明することとした。

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

Fig. 1 Schematic representation of scale-downstirred vessel.
数値解析手法として、２段階で計算する方法を採用した。１段階目には、流動場を定常解析し、２段階目には、1 段階目で計算した流体速度と乱流渦粘性を用いて物質拡散の非定常解析を行った。
今回乱流 Schmidt 数は一定であると仮定した。1
段階目では Reynolds-Averaged Navier-Stokes(RANS)式を用いた乱流解析を行い、撹拌翼の回転は Multi Reference Frame (MRF)法を利用して、撹拌翼付近は回転座標、翼遠方では静止座標系で計算を行った。RANS 方程式を解くために、k- SST モデルを利用し、総格子点数はおおよそ125 万点であった。全てオープンソースである
OpenFOAM を用いて計算を行った。

3 結果と考察

Fig. 4 に示す。撹拌翼が 20 mm の高さにある場
合、均一混合時間が増大しており、効率が悪化す

無次元化したトレーサー濃度の時間変化を Fig.

る。一方で、10 mm, 30 mm の場合には均一混合

2 に示す。最初に各測定点で濃度のピークが発生

時間が減少し、効率が改善された。これは、撹拌する。これは、最初に投入したトレーサーの濃度

翼の位置を変更したことによって流動パターンがが高く、その濃度の高い部分が測定点を通過した

変化したことに起因する。詳細なメカニズム等にためである。最終的にはトレーサーの濃度は均一

関しては著者らの既報[3]を参照されたい。

になるため、
時間と共に一定値へと近づいていく。
数値シミュレーション結果と実験結果を比較すると、ピーク値等のズレはあるものの定性的には同様の分布となる。この実験ではトレーサーの入れ方に依存してばらつきが大きくなってしまうが、おおよそ同様の結果になったため、計算結果が妥当であると言える。

Fig. 4 Effect of impeller immersion depth on theperfect mixing time.

4 結言
本研究ではアルミニウム溶解炉における機械撹
拌効率の定量評価を行うため、水を用いたスケー
ルダウン実験とそれに対する数値シミュレーショ
Fig. 2 Time variation of tracer concentration.

ンを行った。数値シミュレーション結果と実験結 続いて、均一混合時間に対する撹拌翼回転速度

果は良好に一致し、水を用いたスケールダウンモの影響を Fig. 3 に示す。標準偏差のついたデータ

デルではシミュレーションの健全性が確認できた。

が実験値、実線が数値解析結果である。回転速度

今後、数トン-数十トンスケールのアルミニウム溶の大きい場合、数値解析結果は標準偏差外となる

解炉での物質移動を数値シミュレーションと実験

ものの、良好に結果が一致している。当たり前でで比較し、アルミニウム溶解炉に用いることのではあるが、回転速度が上昇するほど均一混合時間

きる数値モデルを確立する必要がある。

が減少する。

引用文献[1] W. Bujalski et al., Chem. Eng. Technol., 27 (2004)
310-314.
[2] J.-F. Bilodeau and Y. Kocaefe, Light Metal, 1(2001) 1009-1015.
[3] T. Yamamoto et al., J. Mater. Process. Technol., 259(2018) 409-415.

謝辞
Fig. 3 Effect of impeller rotation rate on the 本研究制度（若手・女性研究者奨励枠）を活用

perfect mixing time.

させて頂いたことを、この場を借りて厚く御礼申 均一混合時間に対する撹拌翼配置高さの影響をし上げる。

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

1 次元量子スピン鎖の非平衡定常状態における温度勾配の数値的研究
八角 繁男 ∗∗

1

京都大学大学院 理学研究科 物理学・宇宙物理学専攻

背景と目的する。また、zz–カップリングの強さ ∆ は µxx に比量子非平衡系の統計力学はいまだ発展途上にあり、統計力学における重要な問題の 1 つである。量子非平衡系は現象としてはありふれているが、その状態

べ十分小さい (∆/µxx ∼ 10−4 ) 場合のみを考える。

着目系の密度行列の時間発展は以下のリンドブラッ

ド型量子マスター方程式で与える。

を微視的に記述する方法は定常な場合でさえ確立されていない。そこで、非平衡定常状態 (NESS) の微視的な状態の記述法を得るため、まず具体的な系で非平衡定常状態の性質を調べることが必要である。
論文 [1–3] では、一様な 1 次元量子鎖における NESSの性質について論じられているが、非一様な 1 次元
量子鎖については [4] でローカル磁場が非一様な場合が調べられているのみである。そのため、非一様ρ̇ = L(ρ),L(ρ) = −i[H, ρ]
})
∑ ∑(
1{ †
Ljs Ljs , ρ .
+Ljs ρL†js −2
j=0,n−1 s=±(5)

式 (5) に表れるリンドブラッド演算子 {Ljs } は式 (6),(7) で与えられる。

な 1 次元量子鎖における NESS についての更なる知0
L0,± = γ exp(∓βL ω/2)σ±,見を得るため、[4] とは異なる非一様性を持つ系につn−1
Ln−1,± = γ exp(∓βR ω/2)σ±.

いて、NESS を求め、その温度分布を調べた。

2(7)はそれぞれ左と右の熱浴の逆温度を表す。このリンドブラッド演算子は、詳細釣り合いを満たしている。

着目系としてスピン間の相互作用が非一様な 1 次

NESS は式 (8) で与えられる。

元 XXZ モデルを考える。スピンの個数を n として、

L(ρNESS ) = 0.

着目系のハミルトニアン H は以下で与える。
n−1
∑

(6)

j
ここで、σ±:= (σxj ± iσyj )/2 であり、また、βL , βR

セットアップ

H=

(4)

(8)

今回は、ランダムに生成したハミルトニアンに対し

Hj + Hint ,(1)

Hj = ωσzj ,
∑ { (
)
}
vjk σxj σxk + σyj σyk + ∆σzj σzk .
Hint =て、式 (8) で与えられる NESS を数値的に求め、その温度分布についての性質を 3 節以降で述べるよう

j=0

(2)

に解析した。

(3)<jk>

ここで、σij は j 番目のスピンに作用する i–方向のパウリ演算子である。系の非一様性は、xx– or yy–カッ
2) で各プリングの強さ vjk をガウス分布 N (µxx , αxx

スピン-スピン間隔毎にランダムに与えることで実現

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

3

平均温度分布式 (1)-(3) にしたがってハミルトニアンをランダムに多数生成し、それらの NESS における各スピンの逆温度を求め、サンプル平均したものを表したのが

図 1: 各スピンの平均の逆温度を示したグラフ。色が濃くなるほど xx–カップリングの標準偏差 αxx は大きい。最も薄いものが αxx = 10−4 , もっとも濃い
ものが αxx = 1 の場合である。使用したパラメタ

セットは n = 10, ω = 1, ∆ = 10
−2

10

−5

, βL = 1, βR = 1.1 である。

−1

, µxx = 10,γ =図 2: 両端から 1 つずつ内側のスピンの温度差と xx–
カップリングの標準偏差 αxx の関係を表したグラフ。横軸は期待値 µxx によって規格化された標準偏差 αxx /µxx ，縦軸は熱浴の温度差によって規格化されたスピンの温度差 |β8 − β1 |/|βL − βR | である。4 つのパラメタセットを計算した。パラメタセット n =

図 1 である。横軸は各スピンのインデックスを表し、10, ω = 1, ∆ = 10−5 , µxx = 10−1 , γ = 10−2 , βL =1, βR = 1.1 を基準として、他のパラメタセットは凡

縦軸は各スピンの逆温度を表しており、また、相互例に載っている値のみ基準のパラメタセットとは異作用の分布の標準偏差 αxx が大きいほど濃い色の線

なる値を採用してる。

での統計平均である。図 1 から、αxx の増加に伴い、かる。なお、αxx /µxx ≲ 10−1 においては、ほぼ αxxで表されている。各 αxx について、約 1 万サンプル

平均的な温度勾配がはっきりと現れることがわかる。

の 2 乗、すなわち分散で温度差が与えられることが
わかった。

4

温度差と非一様性の関係
両端から 1 つずつ内側のスピン間の温度差を相互作用のゆらぎ αxx に対してプロットしたグラフが図2 である。横軸が相互作用の分布の期待値 µxx で規謝辞
本共同研究制度を活用させて頂きましたことに感謝致します。

格化されたゆらぎ αxx /µxx 、縦軸が両端から 1 つず

つ内側のスピン間の温度差を熱浴の温度差で規格化した量である。スピンの個数 n = 10 において、4 種類の異なるパラメタセットについて調べた。図 2 において、4 種類のパラメタセットにおける全てのデー
タ点は同一曲線上にあり、このことから、温度差とαxx の関係はその他のパラメタに依らない普遍的な性質であると言える。また、図 2 から、αxx が大きくなるにつれて温度差が大きくなるが、αxx /µxx ∼ 1

付近で 2 つの熱浴の温度差に saturate することがわ

参考文献
[1]

T. Ishida and A. Sugita, J. Phys. Soc. Jpn. 85,
074006 (2016).

[2]

T. Yuge and A. Sugita, J. Phys. Soc. Jpn. 84,
014001 (2015).

[3]

T. Prosen, Phys. Rev. Lett. 107, 137201 (2011).

[4]

W. Roeck, A. Dhar, F. Huveneers, and M. Schütz,
J. Stat. Phys. 167, 1143–1163 (2017).

京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 39非一様な都市構造物上における大気乱流の組織構造に関する数値解析
吉田
∗

1

敏哉 ∗

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

はじめに都市域において、運動量・スカラー輸送を担う乱流組織構造の理解はヒートアイランドや大気汚染といった都市特有の環境問題の解決に重要である。都市キャノピー近傍の roughness sub-layer では都市構造物が乱流組織構造に大きく影響する。実際、組織高さの標準偏差と平均高さ Have の比 Vh が Vh = 0.0(V00), Vh = 0.5 (V05), Vh = 1.0 (V10) となるようなブロック列を用いた。ブロックの密度を示す建ぺ
い率 λp は 0.25 としている。本研究では流入境界に
別計算領域で作成した乱流境界層流をタイムステッ
プ毎に与えている。

構造に伴い発生する運動量輸送はキャノピー外ではejection が卓越するが、キャノピー内では sweep が卓越することが知られている [1]。さらに、高さのばらつきを有するキャノピーでは sweep の割合がキャノピー近傍で増加することが示されている [2]。しかし、既往研究では基本的に ejection と sweep の空間
平均的な特徴を示しており、ejection や sweep を伴う組織構造の空間的な特徴と都市構造物の関係は十分明らかになっていない。特に実際の都市上の構造物は非一様な高さ分布をしていることから、建物高
さのばらつきによる乱流組織構造への影響を理解す
ることは重要と考える。
そこで本研究では高さのばらつきを有する粗度ブ
ロック列を使用した Large-eddy simulation を実行
し、建物高さのばらつきがおよぼす乱流組織構造への影響を調べた。特に ejection と sweep を伴う乱流組織構造と都市構造物との関係を調べた。

2

計算設定

Fig 1: block layouts in the unit area of (a) V00, (b)V05, and (c) V10.
Fig. 2a に sweep と ejection の差 ∆S の鉛直プロ
ファイルを示す。∆S > 0 の場合、sweep が卓越していることを示す。図から、キャノピー内では sweepが卓越し、キャノピー高度から上空に向かうに従い、

ejection が卓越することが分かる。また、Vh が大き
いほど sweep の寄与が増加している。これらの結果
は既往研究と整合的である。
続いて、ejection もしくは sweep を伴う乱流組織
構造の空間スケールを算出するために、条件付き平均を用いた。ここでは、(1) 運動量フラックスの瞬間場から平均レイノルズ応力の 2 倍以上の値を示す
ejection(sweep) の点を検出、(2) その点を中心とした風速変動成分 (u′ , v ′ , w′ ) の水平分布を抽出、(3) 同

計算領域は主流方向 4 km・スパン方向 2.4 km で、様の操作をスナップショット毎に繰り返し、得られ水平解像度は 2 m とした。計算領域内には Fig. 1 でた水平分布をアンサンブル平均、という手順で実行示されたような粗度ブロックが繰り返し並べられてした。

いる。各ブロックは水平方向に 10 格子で解像されてFig. 3 は V00 と V10 のケースにおける条件付き平均によって抽出された高度 2Have での ejection 及びいる。建物高さのばらつきの影響を調べるため、建物

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

Fig 3: Vector ﬁelds of conditionally averaged u′
Fig 2: Vertical proﬁles of ∆S in (a) V00, (b) V05,
and (c) V10. The black-solid, -dashed and -dottedlines indicate the height of blocks in V00, V05, andand v ′ at z = 2Have in (a) V00 and (b) V10 withejection and (c) V00 and (d) V10 with sweep. The
contours indicates distributions of u′ normalized by

V10, respectively.

the center values.

sweep 周りの乱流構造である。ここで、高度 2Haveを示した。

は V00 のケースの場合キャノピー上にあたり、V10のケースの場合キャノピー内にあたる。図から V00のケースで得られた ejecction 及び sweep 周りの乱流組織構造は対渦を形成しており、都市構造物群上で見られるヘアピン渦に対応する構造と考えられる。また、建物高さのばらつきを有する V10 のケースでも中心付近に同様な対渦構造が形成されている。ただ

参考文献[1] Coceal, O et al, ”Structure of turbulent ﬂow
over regular arrays of cubical roughness,” J
Fluid Mech 589 (2007), pp. 375-409

し、V10 の場合は下流側 (x′ = −6Have ∼ −2Have )[2] Kanda, M, ”Large-eddy simulations on the effects of surface geometry of building arrays

スケールの対渦はキャノピー上空の乱流構造に由来

on turbulent organized structures,” BoundaryLayer Meteorol 118 (2006), pp. 151-168に一回り大きな別の対渦を確認できる。この大きいするものと考えられる。以上から、建物高さのばらつきをもつキャノピーの場合、キャノピー内の運動量輸送は大スケールと小スケールの渦により行われていることが示唆された。

3[3] Inagaki, A, et al, ”Large-eddy simulationof coherent ﬂow structures within a cubical canopy,” Boundary-Layer Meteorol 142(2)(2012), pp. 207-222

まとめ
本研究では建物高さばらつきがおよぼす乱流組織

構造へ及ぼす影響を調べるため、高さのばらつきが異なる構造物群を使用した Large-eddy simulation を実行した。そのために、ejection と sweep に対応した乱流組織構造を条件付き平均により抽出した。その結果、建物高さのばらつきをもつキャノピーの場
合、キャノピー内の強い運動量輸送は大スケールと小スケールの対渦により引き起こされている可能性

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

半無限領域のスペクトル法による竜巻を模した渦の数値実験に向けた研究開発
松嶋 俊樹 ∗∗

1

理化学研究所計算科学研究センター

はじめに

ここで, u : 速度ベクトル, P : 圧力関数, F : 外力で回転する環境場の中に鉛直流を強制したときに生成される竜巻を模した鉛直渦は, スワール比に依存して多様な渦の形態 (1,2 セル構造, 螺旋渦, 多重渦)をもつことが知られている. 竜巻を模した渦の形態の変化に関する数値実験は, これまで室内実験の環境を模した設定である Ward 型の設定 (Ward, 1972)
や, 剛体回転する円柱容器の内部に浮力強制を与える Fiedler 型設定 (Fiedler, 1993) を採用した研究が行われてきた. どちらの設定も鉛直方向には流出境
界条件や固体壁を置いた有限領域の設定である. しかし, 竜巻を模した渦は, 上昇流を強制しようとしても, 角運動量が駆動する循環が卓越することにより下流側からの下降流が生成されるため, 有限領域の計算では流れは上部境界の影響を強く受ける問題がある. Matsushima and Ishioka (2017) は, 竜巻を模した渦などの大気中の微細渦の数値実験を想定し, 上
部境界を排除した半無限領域を計算領域とする新しいスペクトル法数値モデルを開発した.
本課題では, そこで開発した数値モデルの改良と二次元分割 MPI 並列化を行い, それを用いて竜巻を模した渦の生成に対する境界条件の影響を評価するある. Matsushima and Ishioka (2017) (以下 MI17と書く) は, ナビエ・ストークス方程式を平板の上方に広がる半無限領域で計算するスペクトル法数値モデルを提案した. そこでは, モデル座標として円柱座標系を採用し, 水平方向には立体写像と球面調和関数を用いたガラーキンスペクトル法を用いる. また,
鉛直方向には写像とチェビシェフ多項式を用いたペトロフ・ガラーキン法を用いて, 速度場の粘着境界条件と同等な条件である渦度の積分制約条件を時間発展に課す方法を用いる. 円柱座標系構造格子を採用した数値モデルでは, 中心軸付近に方位角方向の格子が密集するため, 差分法を用いると CFL 条件が厳しくなるという問題がある. 一方, この数値モデ
ルでは球面調和関数展開を用いるため, 波数切断によって CFL 条件を緩和できるという利点を持つ. また, 方位角方向に対して対称性を課すことができる
ため, 軸対称実験や高対称性実験などより理想的な数値実験が可能である. さらにスペクトル法である
ため差分法に比べ少ない自由度で精度の良い計算が行え, 渦の力学や境界層過程に対してより信頼性の高い議論を行うことができる.

数値実験を行った.

32数値モデルの概要

MI17 の数値モデルはスペクトル法モデルである

支配方程式は, 非圧縮ナビエ・ストークス方程式に対して外力を加えた以下の式である.

∇ · u = 0,Du
= −∇P + ν∇2 u + F .
Dt

並列化ため, ルジャンドル変換, フーリエ変換, およびコサ
イン変換からなるスペクトル変換が最も計算量が多く実行時間の多くを占める. スペクトル変換を行う

(1)
(2)

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

ためには, 各変換の際に変換を行う軸方向のデータが全て同じプロセスにある必要があるため, スペク
トル変換を MPI 並列化するには三軸のうち一次元

方向あるいは二次元方向にデータを分割しておき,

基底変換のために必要なテーブルを重複して確保すデータ分割しない軸方向にスペクトル変換を行ってる必要があるため, 必要なメモリ量が増える. チェビ

は通信を行いデータ分割した軸方向のデータを集め

シェフ次数を分割したデータに対して, 隣接するプる必要がある. この際, MPI Alltoall(v) による全対ロセスのデータを取得するスペクトルのハロ通信も全通信が必要になる. 二次元方向に分割する利点と必要となる. そのため, 図は省略するがこれらの計算して, 同じプロセス数による一次元 MPI 並列化に比は鉛直分割のプロセス数を大きくするほど計算時間

べ, MPI Alltoall(v) で一度に動作するプロセス数がが長くなる.

少なくなるため通信量が削減されること, 並列数が増大することがある.
これまでの MI17 の実装では一次元分割の MPI 並列化をしていたが, 竜巻を模した渦の三次元数値実験を行うための検討の結果, 一次元分割では Flat MPIの場合に十分な並列数がとれないことが分かった. これは, 竜巻を模した渦の数値実験では, 下端に粘着境界条件を課すと境界層に強い吹き込みジェットが生成され CFL 条件が厳しくなるため, 時間刻み幅を非常に小さくとらなければならず, 高解像度にできないことによる. そのため, 本課題では動径・鉛直二次元分割の MPI 並列化を行った. スペクトル変換の実装の難しさは, 解く問題サイズや並列数によって通信や転置計算が高速になる実装方法が異なるため,

図 1: スペクトル正変換 (forward) と逆変換 (back-ward) を総プロセス数 36 に固定して行ったときの実行時間に対する動径分割のプロセス数依存性.

複数の実装やチューニングが必要になる場合がある
ことにある. 本課題では軸対称計算用と三次元計算用のサブルーチンを分けて用意し, 軸対称計算用サブルーチンでは不要な計算やデータ移動を除外する
ように, 三次元計算用サブルーチンではループのブ
ロッキングを行い一度の通信で送受信するメッセー
ジサイズを調整できるようにした.
方位角・動径・鉛直格子点数をそれぞれ 512, 256,

512 として三角切断を用いたスペクトル正変換と逆変換をシステム B の一ノードを用いて行ったときの実行時間を図 1 に示した. 二次元分割の各プロセス数が同じになる動径プロセス数 6 の場合に最も計算時間が短いことが分かる. 動径分割のプロセス数 6に対して非対称な分布となっているのは, 鉛直分割したことによって, 各プロセス内でのルジャンドル変換の計算時間が不均等になるからである.
二次元分割によりスペクトル変換は高速化されるが, MI17 に特有な欠点もある. MI17 では, 水平方向の基底関数を時間発展を解くときにより都合の良い基底関数に変換する必要があるが, 鉛直分割によってその計算途中に MPI Alltoall(v) による通信が新たに必要となる. また, 一度に基底変換できる鉛直層数が減ることによって, 実効効率が減少する. さらに,

4

数値実験
開発した数値モデルで竜巻を模した渦が生成されるようにするため, Fidler 型設定に準じた設定を用いる. ナビエ・ストークス方程式の外力に, 浮力強制項と遠方の流れを剛体回転に近づける強制を与える:
[
)2 ]
(
1
ez
F =5.05 exp −10r2 − 20 z −2−Λ(u − Ωreθ ).
η(3)

ここで, Ω は剛体回転率, e は円柱座標系 (r, θ, z) における単位ベクトルである. また, 第二項の強制はvolume penalization 法による表記を用いており, ηは浸透率, Λ はマスク関数である. 初期は剛体回転に微小な方位角方向に依存する擾乱を与えた流れとし,
下端境界条件は粘着境界条件あるいはすべり境界条件とする. 以上の設定に対して, Ω = 0.2, ν = 4 × 10−4として, 時間刻み ∆t = 0.0002 で BDF2(粘性項) とBX2*(非線形項) による時間積分を t = 60 まで行った. なお, この数値実験はそれぞれシステム B の 4

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

ノード 144 プロセス (動径・鉛直分割数はそれぞれ

12) で約 3 日経過した.
粘着境界条件を課した場合における t = 40 から

t = 60 まで平均した方位角平均接線風およびストー
クス流線関数と時刻 t = 59.25 における渦度分布とをそれぞれ図 2, 図 3 に示した. 接線風は鉛直に伸びた構造をしており, (r, z) = (0.15, 0.05) 付近で最大値をとっている. ストークスの流線関数の 0 線は

z = 0.05 付近で r = 0 の境界と交わっており, この線の左側 (反時計回りの循環) と右側 (時計回りの循

図 2: 時間方位角平均接線風 (トーン) とストークス

環) のセルで流れが分かれている. ストークスの流線の流線関数 (コンター)

関数の 0 線が下端境界の近くで r = 0 と交わるこのような流れはいわゆる drowned vortex jump と呼ばれる跳水現象に類似した流れとなるレジームであると判別できる. この平均的な流れに対して, 下端境界付近の渦の内部コアでは動径から方位角方向に傾いた非軸対称擾乱が過渡的に発達する. 下端境界から少し離れた z > 0.15 では, 鉛直・方位角方向に傾いた方位角波数 1 型から波数 4 型の擾乱が過渡的に発達する. さらに z > 0.4 では渦崩壊によって水平に軸をもつ渦の発達する乱流領域になっている.
すべり境界条件を課した場合の, 図 2, 3 と同様の

図 3: 時刻 t = 59.25 における渦度分布

量を図 4, 5 に示した. 下端境界付近から z ∼ 0.3 まで, 接線風は鉛直に伸びた構造になっている. ストー
クスの流線関数の 0 線は r = 0.18 付近で z = 0 の境界と交わっており, いわゆる 2 セル構造をもつ流れのレジームであると判別される. 鉛直一様に近い下端境界付近から z ∼ 0.3 までの領域では, 非軸対称擾
乱の発達は粘着境界条件を課した場合より弱く, 渦の中心軸で渦度が小さくリング状の渦度分布をしている. 下端境界から離れた z > 0.4 では粘着境界条件を課した場合と同様に渦崩壊によって乱流領域となっている.

5

図 4: 時間方位角平均接線風 (トーン) とストークスの流線関数 (コンター)

まとめと今後の課題
半無限領域のスペクトル法の改良と並列化を行い,それを用いて竜巻を模した渦に関して境界条件の影響を調べる数値実験を行った. 本課題を利用して行われた二次元分割 MPI 並列化によって, 数値モデルの大幅な高速化に成功したため, 今後は動粘性係数,
剛体回転率, 境界条件のパラメータを変えたパラメー
タスイープ実験を行い渦の力学を調べる.
図 5: 時刻 t = 59.25 における渦度分布

44 − Academic Center for Computing and Media Studies, Kyoto University高周波分数冪ラプラシアン Navier-Stokes 方程式のエネルギースペクトルの考察
中井拳吾
東京大学 数理科学研究科

本稿は岸本展(京都大学数理解析研究所), 斉木吉隆(一橋大学商学研究科), 米田剛(東京大学数理∂tv + νPΛ(-Δ)A/2v + (v・∇)v + ∇π = F

科学研究科)らとの共同研究に基づくものである.

外力, PΛ(-Δ)A/2:=ℱ-1χ{|y|>Λ}|y|A ℱ, χ:定義関数

ただし, ν:動粘性係数, v:速度, π:圧力, F:とする. このように一部の周波数に対して粘性

1 はじめにの効果が入るように一般化した粘性を部分粘性と呼ぶこととする. また, A=2, Λ=0 とした(GNS)

Navier--Stokes 方程式は, 乱流を含めた流体運動を記述できることが知られている. 特に, 3は Navier--Stokes 方程式である.
直接数値計算では, 境界を 3 次元周期境界と次元乱流では外力によりエネルギーが注入されると, 流体要素がねじられたり引き伸ばされたし, (119×2+1)³ に細分化した格子点データによるフーリエスペクトル法を用いる. ただし, 離散りすることによってより小さなスケールの運動を引き起こす. つまり, 注入されたエネルギー

フーリエ変換は高速フーリエ変換(FFT) というアルゴリズムを利用するが, FFT により現れるる. このエネルギーが遷移する現象はエネルギ

エイリアスエラーは 2/3 ルールに基づいて除去する. また, 4 次精度の Runge-Kutta 法を基本

ーカスケード [2]と呼ばれている. しかし, エネ
ルギーがどのようにして遷移していくかのメカ

的には 1/800 秒刻みで時間発展をさせる. ただし, エネルギーの総和が発散する傾向にある場合

ニズムは複雑であり厳密に解明されていない.
さらには, エネルギーカスケードは任意の速は適宜時間刻みを細かくするアダプティブメッ
シュを用いる. 初期条件は適当ななめらかな速度

度場の時間発展の過程で観測されるわけではな
い. 例えば, Euler 方程式の定常解として知られ場とし, 定常乱流を得るために低周波領域のいく
つかの波数に外力を入れる.

ている Arnold--Beltrami--Childress Flow (ABCFlow) [1]はエネルギーカスケードが見られない.

3 計算結果

そこで粘性の効果に注目して流体のエネルギーの様子を考察することにした.

(1) 低周波領域の粘性の効果を除去した(GNS)のが小さなスケールに遷移することが知られてい考察
各波数 k でのエネルギーE(k)を次で定義する:

2 流体
Navier--Stokes 方程式を用いてエネルギーカ

E(k, t) = � |ℱ[�] (𝐾𝐾, 𝑡𝑡)|𝑑𝑑𝑑𝑑
��

スケードのメカニズムを直接考察することは任意の周波数に粘性の効果があるため困難である.

ただし, 積分範囲は Dk:={K∈ℝ³|k-0.5<|K|<k+0.5}と定義した. Λ=0 を固定し A=2,5,10 の

そこでまずは Navier--Stokes 方程式の粘性項を一般化した次の方程式(GNS)を考察する.

場合のエネルギースペクトル E(k)の違いを上図,
A=5 を固定しΛ=0,100 の場合の E(k)の違いを下京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 45

図に書き出した. 横軸が波数 k, 縦軸がエネルギーを表し, 参考のためコルモゴロフ則を表す k-5/3

き起こされていることがかる.

の直線を描いた.
いずれの図の E(k)も波数 k ∈[5,100]の領域ではコルモゴロフ則の線に沿っていることが見て取れる. つまり, 部分粘性や超粘性の場合についても通常の Navier--Stokes 方程式と同様にコルモゴロフ則に従っていることがわかる.

図 2: (A,Λ)=(2,0),(5,100) に対する時刻t = 0,10,35, 65 におけるエネルギースペクトルE(k, t) を両対数メモリでプロットした. ただし, 各パラ
メータにおいて十分に発達した乱流を初期値とした. また, 実線はコルモゴロフの5/3 則を表
す.

4 まとめ
これらの(1)(2)のことから, Navier--Stokes 方
程式の粘性項を変化させてもエネルギースペク
トルには影響しないことと, 低周波領域の粘性の効果の有無にかかわらず低周波部分はエネル
ギー減衰することが分かった. 通常の PC では十分大きな波数までは計算できないため E(k)とコルモゴロフ則の関係を見ることは困難であるが,スーパーコンピュータの利用によりコルモゴロ
図 1: 各パラメータA,Λに対するエネルギースペ

クトルE(k, t)を両対数メモリでプロットした.

フ則との比較を考察することができた.

ただし, 各エネルギースペクトルはt = 10 からt

5 まとめ

また, 実線はコルモゴロフの5/3 則を表す.

[1] T. Dombre, U. Frisch, J. M. Greene, M.

= 30 の値の時間平均をとって得たものである.

H_enon, A. Mehr, and A. M. Soward. Chaotic(2) 減衰乱流の考察
十分発達した乱流を初期値として外力をなくした場合の(A,Λ)=(2,0),(5,100)について時間発展をさせ, 時刻 t=0,10,35,65 における E(k)を書き出した.
この結果から低周波部分ではコルモゴロフ則の線と並行を保ちながら時間減衰していることが見て取れる. このことから, 高周波にのみ粘性を入れた GNS でも低周波部分のエネルギー減
衰は, 通常の Navier--Stokes 方程式と同様に引46 − Academic Center for Computing and Media Studies, Kyoto University

streamlines in the abc ows. Journal of Fluid

Mechanics, 167:353--391, 1986.

[2] L. F. Richardson. Weather prediction by

numerical process. Cambridge University
Press, 1922.

プログラム高度化共同研究報告

動的／静的水～土骨格連成有限変形解析コードの高度化～任意形状に対する領域分割法の適用の検討～
野田利弘

豊田智大

名古屋大学 減災連携研究センター

名古屋大学大学院 工学研究科

1 平成 29 年度事業の概要と目的

2 検討事項の概要

筆者らが開発した地盤の動的／静的変形解析コード（GEOASIA® [1], [2]）に対して、平成 24 年度（下図１は、計算の全体処理フローを示すともに、今回の検討での要点を以下に示す。

半期）[3], 25 年度[4], 26 年度[5], 28 年度[6]に続き、平成 29 年度も京都大学学術情報メディアセンターのプログラム高度化支援事業に採択していただいた。
この概要を報告する。
本解析コードは、粘土から中間土、砂までの広範な土の力学挙動を同一の理論的枠組みで記述する弾塑性構成式（All Soils）を搭載し、地盤の変形から破壊までの挙動（All States）を、動的と静的を問
わずあらゆる外力条件（All Round）の下で解析することを目指したものである。これに対し、本事業の目的は、OpenMP と MPI とのハイブリッド並列化、反復法ソルバである PETSc[7]の組み込みや、領域分割法による並列処理による計算手法などの導図１ 全体処理フロー

入により、大規模・高速化計算をすべく当該解析コードの高度化・高性能化を図ることである。
平成 29 年度の目的は、平成 26 年度（以下、「前(a)任意形状に対し計算可能となる処理方法・フローの検討版」）の領域分割法について、計算可能な有限要素形

領域割当：前版では、ddm_input 内で分割数、位状を拡大することである。つまり、前版では、対称置座標に基づき領域分割を行い、節点、要素を割り的な形状に対しては領域分割法が適用可能であった

当てた。これに対し、今回は ddm_input で行わず、

ものの、非対称形状などの任意の形状に対しては適別の手段（本検討では GUI ツール）で行い、割り用できなかったため、より任意な形状で計算ができ当てられた領域（各領域の要素番号）に基づきるようにした。
なお、
領域分割法による並列処理は、

ddm_input で付番変更などの改訂処理を行った。

計算対象領域を複数領域に分割して各分割領域をプ

節点番号の付番変更：領域毎に要素番号の順に行

ロセス並列化手法により並列処理し、領域境界間でい、それに基づきコネクティビティの更新も行った。

必要な値の送受信をプロセス間通信により行う計算(b) GUI ツール使用での手順検討手法であるが、前版から解析コード自体は変更せず、GUI ツールとして計算結果の可視化などに利用

前処理を中心に計算フローの変更を行った。本稿でしている Paraview を使用することでメッシュデー

は、この概要について紹介する。

タを可視化し、分割箇所を GUI で取捨選択して分京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 47

割割当を行った。
（具体の処理については省略。
）

意形状への適用を図るため、全体処理フローの中で、(c)前処理（ddm_input.f90）の改訂

GUI ツールによる領域割当を採用した。

要素を領域ごとに割当を行い、節点は ddm_input

GEOASIA 動作検証では、PETSc 収束傾向／ソにおいて領域順、要素順、コネクティビティ順に付

ルバ以外の送受信箇所などの動作確認として、変位番を行った。

や応力値についてオリジナルシリアルコードと掲載結果の比較を行ったところ、有意な差異は生じてい3 領域割当検討例と計算結果比較ないことを確認した。
前版の検証[5]では、ソルバに PARDISO を用いた

図２は領域割当検討例を示す。領域を対称的な形場合の比較として、PETSc の場合に、PARDISO 比状としないため、上下分割において、領域境界をあで約 200～1000 倍の高速化が図られている。本検討

えて凹凸にして 2 つの領域に割り当てた。

では簡単な形状に適用したが、この高速化は魅力的なため、今後は複雑な３次元要素形状、領域分割による動作検証や、前版において PETSc で動作できていないケースへの対処も含め、各種入力データでの動作検証などを是非、実施してゆきたい。
謝辞 京都大学学術情報メディアセンターの中島浩教授および牛島省教授ならびにクレイ・ジャパン・
インクの方々には、上記以外にもご尽力いただきました。ここに深く謝意を表します。

図２ 領域割当の例（上下 2 分割ケース）
図３はオリジナルシリアルコードと図２の 2 分割
ケースについて可視化による計算結果比較である。
変位や応力の数値による計算結果比較も行ったが、計算結果に実質的な差異は生じなかった。

参考文献[1] Asaoka, A. and Noda, T. (2007): All soils all states allround geo-analysis integration, International Workshop on
Constitutive Modelling - Development, Implementation,
Evaluation, and Application, Hong Kong, China, 11-27.
[2] Noda, T., Asaoka, A. and Nakano, M. (2008):
Soil-water coupled finite deformation analysis based on a
rate-type equation of motion incorporating the SYS
Cam-clay model, Soils and Foundations, 48(.6),
771-790.
[3] 野田利弘 (2013): 京都大学学術情報メディアセ

図３ 計算結果比較

ンター全国共同研究利用版広報，Vol.12, No.2, 33-35.
[4] 野田利弘 (2014):京都大学学術情報メディアセ

また、2p8t（2 プロセス 8 スレッド）での性能比

ンター全国共同研究利用版広報，Vol.13, No.2, 44-47.

較については、前版が 110.23 秒であったのに対し、[5] 野田利弘 (2015): 京都大学学術情報メディアセ

改訂版が 118.14 秒でほぼ同速であり、ddm_input

ンター全国共同研究利用版広報，Vol.14, No.2, 31-34.

の改訂に伴う性能への影響は、この比較では生じて[6] 野田利弘・豊田智大(2017): 京都大学学術情報メいなかった。

ディアセンター全国共同研究利用版広報，Vol.16,

4 おわりにNo.1, 62-65.
[6] Intel MKL（LAPACK）, GGE: https://software.intel.
com/en-us/mkl-developer-reference-c-ggev.

平成 29 年度は、領域分割法による並列処理の任

48 − Academic Center for Computing and Media Studies, Kyoto University[7] PETSc/Tao: https://www.mcs.anl.gov/petsc/.

高分子分離膜の大規模シミュレーション
矢野 恭平 1)

安武 遼 1)

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

1 はじめに

吉元 健治 1,2)

２）京都大学学際融合教育研究推進センター

2 モデル

下廃水処理に用いられる高分子分離膜は多数の本モデルでは、水、溶媒、ポリマーの三成分の細孔を有し、そのサイズや形状等で分離性能が決

局所体積分率𝜙𝜙� �𝒓𝒓�、𝜙𝜙� �𝒓𝒓�、𝜙𝜙� �𝒓𝒓�に加え、高分定する。代表的な分離膜作製方法に非溶媒誘起相

子鎖の絡み合い構造の局所歪み𝑾𝑾�𝒓𝒓�も変数とし分離(NIPS)法がある(図 1)。NIPS 法では、まず疎て取り扱う。𝒓𝒓𝒓は位置ベクトル、𝜙𝜙� �𝒓𝒓�はスカラー、

水性高分子を高分子にも水にも可溶な溶媒に溶かし、その溶液を水槽に浸漬させる[図 1(a)]。する

𝑾𝑾�𝒓𝒓�は二階テンソルである。以下では、𝜙𝜙� �𝒓𝒓�を
𝜙𝜙� 、𝑾𝑾�𝒓𝒓�を𝑾𝑾として表す。系の自由エネルギー𝐹𝐹と、高分子溶液中の溶媒が水槽へ、水槽中の水がには、従来の Flory-Huggins の混合エネルギー

高分子溶液へ流れ込み、溶媒置換が生じる[図𝐹𝐹��� と界面エネルギー𝐹𝐹��� に、式(1)で示される弾

1(b)]。水を含んだ高分子溶液は、高分子濃厚相と水濃厚相に分離し、水濃厚相は乾燥後に細孔となる[図 1(c)]。NIPS 法で形成される高分子膜の特徴
として、膜表面にはナノメートルスケールの細孔を、膜内部にはマイクロメートルケールの細孔を有する点が挙げられる。前者は分離性能を、後者は透過性を決定する。
近年、NIPS 法で形成される高分子膜内の細孔構造を、相分離シミュレーションで予測する研究が行われているが、実際の高分子膜内で見られるような非対称構造はまだ再現されていない[1,2]。
そこで、著者らは、粘弾性を考慮した三成分系での相分離モデルを新たに開発し、非対称構造の形成における粘弾性の重要性を示唆した[3]。ただし、粘弾性相分離モデルは計算負荷が大きく、大規模系でのシミュレーションが困難であった。
そこで、本共同研究では、特に負荷の高い速度場計算を高速化し、大規模系での数値計算を試みた。
溶

溶

溶

高分子
高分子 液
(a) (高分子+ 媒) (b) 水
液 水
溶

溶

水

水

媒

水

(c)

媒

図 1 NIPS プロセスの概略図

性エネルギー𝐹𝐹�� を加えた [3-5]。
𝐺𝐺
𝐹𝐹�� � � �𝒓𝒓 ��𝑾𝑾�� � 𝜹𝜹�� �� 𝒓𝒓𝒓𝒓𝒓�1�
4���

𝑉𝑉は系の体積、𝜹𝜹�� はクロネッカーのデルタである。
せん断弾性率𝐺𝐺は高分子の局所体積分率の関数として、𝐺𝐺 � 𝐺𝐺� 𝜙𝜙�� で表した[5]。なお、局所歪み𝑾𝑾によって生じる弾性応力𝜫𝜫��𝒓 は𝜫𝜫�� � 𝐺𝐺𝒓𝑾𝑾 𝛁 �𝑾𝑾 � 𝑰𝑰�𝒓𝒓𝒓𝒓𝒓���

より求められる[4]。
シミュレーションでは、局所体積分率 𝜙𝜙� (𝑖𝑖=1,2, 3)と局所歪み𝑾𝑾を，それぞれ連続の式と構成方
程式を用いて時間発展させた[4,5]。
𝜕𝜕𝜕𝜕�� �𝛁𝛁 𝛁 �𝒗𝒗� 𝜙𝜙� �𝒓𝒓𝒓𝒓𝒓���
𝜕𝜕𝜕𝜕

𝜕𝜕𝜕𝜕
1
� 𝑫𝑫� 𝛁 𝑾𝑾 � 𝑾𝑾 𝛁 𝑫𝑫�� � �𝒗𝒗𝟑𝟑𝟑 𝛁 𝛁𝛁�𝑾𝑾 � 𝒓�𝑾𝑾 � 𝑰𝑰��4�𝜕𝜕𝜕𝜕𝜏𝜏

𝑡𝑡は時間、 𝒗𝒗�𝒓 は各成分の速度場を表す。式(4)に含まれる𝑫𝑫� は𝒗𝒗� の勾配テンソル、𝑰𝑰𝑰は単位テンソル、𝜏𝜏𝜏は緩和時間である。各成分の速度場 𝒗𝒗�𝒓 は
𝒗𝒗� � 𝒗𝒗 � 𝐿𝐿� �𝛁𝛁𝛁𝛁� � � ��� 𝐿𝐿� �𝛁𝛁 𝛁 𝛁𝛁�� �𝜙𝜙� �𝒓𝒓𝒓𝒓𝒓���
より計算した[6] 。溶液の平均移流速度 𝒗𝒗 はNavier-Stokes 式より(次節参照)、化学ポテンシャ
ル 𝜇𝜇�𝒓 は系のエネルギーの汎関数微分δ𝐹𝐹𝐹𝐹𝐹𝐹� より求めた。𝐿𝐿� は成分𝑖𝑖の易動度で、体積分率の関数である[5]。

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

3 シミュレーション

3 に示すように、コア数を増やすことで計算時間は減少したが、どのコア数でも計算時間の五割近

3.1 計算のフロー

くは速度場計算で占められた。図 4 は、速度場の計算フローを図 2 に示す。基本的には、拡散方計算時間を１コアで要した計算時間で割った値と程式を逐次的に差分法で解く流れと同じである。

アから 64 コアまでは、計算時間がほぼ傾き-1 で素（もしくはセル）に三成分の体積分率𝜙𝜙� と歪みテンソル𝑾𝑾の初期値を割り当てる。次に、各セルにおける三成分の化学ポテンシャル𝜇𝜇�� (=δ𝐹𝐹𝐹𝐹𝐹𝐹� )と、応力テンソル𝜫𝜫�� を式(2)より計算する。その減少している。ただし 128 コア及び 256 コアでは、𝑥𝑥方向の厚みが２セル及び１セルになるため、コア間でのコミュニケーション時間の比率が増えてしまい、性能向上がほぼ飽和した。

計算結果を用いて、ストークス流れを仮定した

𝜂𝜂�𝒗𝒗 � 𝛁𝛁𝑝𝑝 � � 𝜙𝜙� 𝛁𝛁𝛁𝛁� � 𝛁𝛁 � 𝜫𝜫�� ��������
�

から溶液の平均速度場𝒗𝒗を解く[6]。𝜂𝜂は溶液の粘度、𝑝𝑝は非圧縮性を保つために生じる圧力である。平

均速度場𝒗𝒗を式(5)に代入して各成分の速度場𝒗𝒗� を求め、𝜙𝜙� を式(3)の連続式で、𝑾𝑾を式(4)の構成方
程式で逐次的に更新する。

60

600

50

500

計算時間[s]Navier-Stokes 式

700

40

400

30

300

20

200

10

100
0

1

2

4

8 16 32 64 128 256
コア数

計算時間の割合[%]

まず、系全体をグリッドに分割し，各微小体積要コア数との関係を対数プロットしたもので、２コ0

図 3 並列化による計算時間の向上(黄：速度場、
赤：コミュニケーション時間、緑：その他)
計算時間（１コア基準）

なお、式(5)から𝒗𝒗を計算する際、まず式(5)をフーリエ変換して波数空間での速度場を求めてから、逆変換により実空間での速度場𝒗𝒗を得ている。次節で示すように、フーリエ変換の計算には MPI
並列版 FFTW ライブラリ(http://fftw.org/)を用い

1

0.1

0.01
1

て、計算負荷の軽減に取り組んだ。

10

コア数

100

1000

図 4 コア数と速度場計算時間の関係

なお、上述の並列化計算では、系を一次元に分割しているが、二次元分割した場合についても検討した。具体的には、多次元分割に対応した
P3DFFT ( http://www.p3dfft.net/)を用いて、256×256×256 のバルク系を対象に、FFT のみの動
作検証と性能評価を行った。図 5 に示すように、
128 コアまでは、一次元分割と二次元分割による

図 2 計算フロー

3.2 計算の高速化
グリッドサイズが 256×256×256 のバルク系を対象に、MPI による並列化計算を行った。計算機には京都大学スパーコンピュータのシステム
A(1 ノードのコア数は 68)を用い、バルク系を𝑥𝑥方
向にコア数分分割した際の計算時間を調べた。図50 − Academic Center for Computing and Media Studies, Kyoto University

１ループの平均計算時間[s]

性能差は顕著には見られなかった。
0.04

一次元分割（合計時間）

一次元分割（コミュニケーション時間）

二次元分割（合計時間）

二次元分割（コミュニケーション時間）

0.02

0.00
コア数

64

128

256

512

64x1 64x1 64x1 64x1 128x1 128x1 128x1 256x1 256x1
２次元分割
8x8 8x8 8x8 8x8 8x16 8x16 8x16 16x16 16x16 16x32
ノード数(コア数/ノード) 1(64) 2(32) 4(16) 8(8) 2(64) 4(32) 8(16) 4(64) 8(32) 8(64)１次元分割

図 5 多次元分割による FFT の並列計算

3.3 数値計算の事例(a) t = 0(b) t = 500(c) t = 5000

まず、三次元バルク系(128×128×128 グリッド)で相分離シミュレーションを行い、相分離構造とせん断弾性率𝐺𝐺� の関係について調べた。高分子と溶媒と水の初期組成は、図 6(a)に示すようにス
ピノーダル領域内に設定した。図 6(b-d) は、相分離が十分に進んだ無次元時間� � ��00での高分子の局所体積分率𝜙𝜙� の空間分布を表し、赤橙色の部分が高分子ドメインを示す。せん断弾性率𝐺𝐺� �

1では、弾性が比較的小さいため、複数の高分子

ドメインが界面張力によって大きく変形し、一つのドメインを合一した[図 6(b)]。一方、𝐺𝐺� を大きくすると、微細かつ複雑に連結した高分子ドメインが形成された[図 6(c), (d)]。これは、界面張力により高分子ドメインが変形する際、弾性応力が大

図 7 NIPS 系での相分離構造の経時変化：(上段)全体像、(下段)上段(a)の赤線で囲まれた𝑧𝑧𝑧𝑧断面像。どちらも高分子の局所体積分率𝜙𝜙� の空間分布を表す。

4 おわりにきくなり、変形の進行が遅くなるためである。

(a)(b) G0 = 1

粘弾性を考慮した三成分系での相分離モデルを開発し、並列化で計算負荷を軽減させ、NIPS 法で作製される高分子膜の形成シミュレーションを行った。その結果、従来のシミュレーションでは考慮されていなかった高分子溶液の弾性が膜内の(c) G0 = 10(d) G0 = 100

細孔形状・分布に大きく影響を及ぼすことが明らかになった。今後は速度場計算及び並列化の最適化を更に進めて、実験データを再現できるような大規模シミュレーションを行う予定である。
本研究の一部はプログラム高度化共同研究によ

図６ 高分子/溶媒/水の三成分混合溶液のバルク系
での三次元シミュレーション：（a）三角相図上での初期組成、（b−d）高分ドメインの相分離構造。

次に、NIPS 系(64×64×128 グリッド)でシミるもので、共同研究者である京都大学学術情報メディアセンターの中島浩教授と牛島省教授、並びクレイ・ジャパン・インクの方々には謝意を表す。

ュレーションを行った。せん断弾性率𝐺𝐺� の値は

参考文献

0.01, 𝜙𝜙� � 0.��, 𝜙𝜙� � 0.�0に設定した。高分子溶[1] B. Zhou et al., J. Membr. Sci., 268, 150 (2006)

100 に固定し、高分子溶液の初期組成は𝜙𝜙� �

液が水と接触するとすぐに、接触界面から溶媒が水中へ、水が溶液中へ一気に流れ込んだ[図 7 (a)]。
その結果、界面付近の高分子が凝集し、層状の高分子濃厚相が瞬時に形成された[図 7(b)]。高分子
溶液の流れ込んだ水は、時間とともに溶液内部に拡散するため、溶液内部ではスピノーダル分解による相分離が徐々に進行していった。最終的には溶液全体で複雑に入り組んだ高分子ドメインの相分離構造が形成された[図 7(c)]。バルク系と同様に、高分子ドメインが弾性により変形し難くなっ[2] D. R. Tree et al., Soft Matter, 13, 3013 (2017)[3]土井正男,“ソフトマター物理学入門”，岩波書
店 (2010)[4] 吉元健治, 矢野恭平, 谷口貴志, 大嶋 正裕,
第 66 回高分子討論会要旨集 (2017)[5] T. Taniguchi et al., Phys. Rev. Lett., 77, 4910
(1996)[6] K. Yoshimoto and T. Taniguchi, in preparation(2018)たことが主な要因として挙げられる。

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

非均質異方性材料評価のための陽的ボクセル有限要素解析の高度化
中畑和之 ∗∗

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

1 はじめに近年，軽量で高強度な炭素繊維補強材料 (Carbon

Fiber Reinforced Plastics, CFRP) の産業利用が進んで

CFRP は樹脂と繊維のミクロ構造までもモデル化する場合は大規模な計算になるため，計算効率を上げるのは至要たる課題であった．

おり，確実な安全性が求められる航空機にも主要部材として使用されている．CFRP はプリプレグと呼ばれるシートを積層して硬化させることで作られる複合材料であり，積層方向や数を変えることで発現する強度が変化する．このとき，特定の方向に大きな強度が発現するように成形すれば，材料異方性が生じる．現状では，この異方性が CFRP の検査を困難にしている．異方性によって超音波の伝搬速度が方向依存性を示すため，これを正しく把握出来ない
場合は，超音波探傷試験 (Ultrasonic Testing, UT) の信頼性を著しく低下させることになる．ここで，数値シミュレーションによって超音波の伝搬がモデル

2

ボクセル有限要素法
陽的に更新する時間領域の FEM 解析について簡単に示す．x = (x1 , x2 , x3 ) 空間を超音波が伝搬するものとし，このときの変位を u(x, t) = (u1 , u2 , u3 )とする．異方性弾性体中の超音波の伝搬を支配する方程式は，以下のようになる．

∂τij (x, t)(i = 1, 2, 3)∂xj∂uα (x, t)(i, j = 1, 2, 3)
τij (x, t) = cijαβ∂xβ
ρ(x)üi (x, t) =

(1)
(2)

化できれば，CFRP の音速分布が分かるだけでなく，¨ は時刻 t に関する偏微分 (∂ 2 /∂t2 ) を表す．
ここで，()

欠陥による散乱波やその伝搬特性も把握できるため，

また，cijαβ は弾性スティフネス，ρ は密度である．

UT のアシストに繋がる．

上式を，重み付き残差法を用いて空間に対して離散過去のプログラム高度化共同研究で，差分グリッ化する．重み関数の任意性から，要素剛性方程式がドを用いる有限積分法 (Finite Integration Technique,

得られ，それらを全要素数分組み合わせると，以下FIT) を改良して頂き，UT のシミュレーションに大きの領域全体の剛性方程式を得ることが出来る．

く貢献できた [1]．しかし，FIT で一般音響異方性材

M d¨ + Kd = f

料を扱うにはグリッド値の補間が必要であり，計算(3)

精度および安定性に懸念が生じる．そこで，異方性

ここで K は剛性行列，M は質量行列である．式 (3)材料に対しては，一定の大きさを有するボクセル要

は時刻 t に関する 2 階の連立常微分方程式とみなす

素を使った有限要素法 (Finite Element Method, FEM)を採用する．ボクセルの角部からの数値振動の影響を少なくするには，最小波長の 1/10 程度の要素長を用いる必要がある [2]．また，この FEM は陽的に時間更新を行うため，時間刻みが陰解法に比べて小さくなり，結果として計算ステップ数が増える．特に，

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

ことができる．
次に，時間軸の離散化であるが，現在の時刻ステッ
プを k とし，その時刻から ∆t 前後の時刻に対応するステップを k − 1 と k + 1 とすれば，現在の時刻の加速度 d¨k は，1(dk+1 − 2dk + dk−1 )
d¨k ≈
∆t2

(4)

と近似できる．式 (4) を式 (3) に代入して展開する．

こで，全体剛性行列 K を作ると膨大なメモリを要す

ここで，質量行列を対角化し，集中質量行列 M̃ をるので，実際には各ノードの更新に必要な K と M

導入すると，変位ベクトル d の成分が独立となるただけを計算するようにしている．まず，材料毎の K

めに陽的な解法となり，単純な代入計算で変位が求と M のテーブルを用意しておき，初期条件を考慮しめられる．

て解の更新を始める．1 回の更新で，matrux-loop1,

dk+1 = (2−∆t2 M̃ −1 K)dk+∆t2 M̃ −1 f k −dk−1 (5)ノードの変位ベクトル d の計算，matrix-mpi は各プ

上式は，現在の時刻ステップ k における値と ∆t 前

ロセスに割り当てられた計算領域の袖部にある d ののステップ k − 1 の値を用いて，変位が更新できることを意味している．

通信，matirix-loop2 は d の更新，matrix-output は dの書き出しを行う．オリジナルプログラムでは，計算領域を複数に分割し，1 つの副領域の計算を複数のコアで実行できるようになっている．すなわち，マ

3 プログラムコードの高度化

ルチプロセスに分割後，副領域で Open MP による

3.1 オリジナルコードのフローと問題点
FEM を実行するためのフローを図 1 に示す．実質的には，図 1 中の matrix と書かれたサブルーチンで，式 (5) で表された変位ベクトルの更新計算を行う．こデータ入力

matrix-mpi, -loop2, -output を実行する．matrix-loop1 は各全セルの条件指定

input

スレッド並列を実施する (Hybrid 実行)．本コードでは，3 次元 MPI 通信を用いて，プロセス間のデータ通信を行っている．また，計算領域を複数に分割し，

1 つの副領域の計算を 1 つのコアで実行するプロセ
ス並列のみの実行も可能である (Flat MPI 実行)．
京都大学スーパーコンピュータ System C で，オ
リジナルコードの実行時間を調べたところ，大部分が matrix-loop1 の計算に要することがわかった．まvcell&src

た，データ数が増えると input に要する時間が大きな割合を占めることも問題であった．

配列確保・変数初期化 mat_pre
材料毎の剛性行列 K と質量行列 M の作成

mat_const

3.2 改良点
3.2.1

UT のシミュレーションでは，欠陥 (自由反射境界)

外力 f の計算
K の M の計算

コンパイラによる最適化

loop1

や探触子の位置などが問題によって大きく変わるため，この位置の判定をボクセル毎に行う必要がある．
これは，すべてのボクセルで判定を実行するため，

変位 d の計算

matrix-loop1 では IF 文を多用するコードになってい

袖部の通信

mpi

変位 d の更新

loop2

全変位 d の出力

outputた．ネストしたループはどれも小さく，アクセスス
トライドも 1 以上であったので，ベクトル化されていなかった．そこで，スカラ変数を連続する配列に変更したり，IF–ELSE 文に変更することで条件分岐数を減らした．こうした変更によって，コンパイラ
による最適化が行われるようになった．

指定位置の変位出力

post
3.2.2

全ノード計算時間更新

図 1: FEM のフローチャート

入力処理の最適化

CFRP は複合材料であるので，すべてのボクセルで異なる材料定数を定義できるようにする必要があ
る．これまでは，各ボクセルの材料定数を逐次的に京都大学 学術情報メディアセンター 全国共同利用版［広報］Vol.17 No.1（ Dec. 2018）− 53

入力ファイルから読み込み，一旦メモリにスプール変位を出力して終了する．1 ノード内で並列数を変した後に，材料定数を各プロセスに割り当てていた．

化させたときの実行時間を図 2 に示す．Hybrid 実行

そこで，MPI-IO により分散して読込みを行うこととよりも Flat MPI 実行の方が高速であったため，ここし，副領域分のデータだけを各プロセスに割り当て

では，Flat MPI 実行のみの結果を示している．最適た．こうした改良によって，入力ファイルの読み込

化版は「-xCORE-AVX2」オプションとブロッキング

み時間と一時メモリを大幅に削減できた．

サイズを 52 に設定した条件での結果である．オリジ
ナル版では，コア数を 7 倍にすれば 4.49 倍程度の実行速度の向上 (Speed up) が得られていたが，改良後3.2.3 MPI 通信の順序変更

は 7.07 倍の Speed up がみられ，非常に良いスケー

本コードでは，3 次元 MPI 通信を用いて，変位ベ
クトル d のプロセス間のデータ通信を行っている．

ラビリティを示すことがわかった．

次元配列に格納されているが，通信のために袖部のSystem C とは異なり，System A(Intel Xeon Phi
7250 KNL, 1.4GHz) のリソースを十分に活用するな
らば，KNL のハードウエアスレッドを使用した Hybrid 実行を行う必要がある．図 3 に，System A にお

データを連続した 1 次元配列に一旦パックして，通いて並列数 (プロセス数とスレッド数）を変化させた

信後にアンパックする方法をとっている．ここで，ノ場合の計算実行時間について示す．また，同図にはンブロッキング通信を採用していても，ブロッキン
グ通信に等しい処理となっていたので，処理の順番

System C の結果も比較のため掲載している．System
C での Flat MPI(72 並列) に比べて，System A の Flatを入れ替える事で通信の隠蔽が有効になるように変MPI(64 並列）では matrix-loop1 の計算時間が大幅更した．

に増加しているが，並列数が少ないことに加えて，副領域の袖部の 1 層について，x1 ，x2 ，x3 方向の変位を隣の副領域に通信する．本コードでは，変位は 3

CPU のクロックが 2.3GHz から 1.4GHz に下がったことが原因として考えられる．1 ノードを用いて評価

3.3 性能比較を行ったところ，Hybrid 実行における最速値は，64

京都大学スーパーコンピュータ System C(Intel Xeion

プロセス ×4 スレッドの 256 並列実行が最も良い性

E7-8880 v3, 2.3GHz) で性能検証を行った．均質・等能となった．スレッド数による性能向上はみられた

方性の数値モデル (要素数は 1500×500×1500) を用

ものの，16 スレッドで頭打ちとなっており，全体のいて検証し，時間ステップ数が 100 となったときに，傾向としてはプロセス数が多い方が良い性能を示す

org
mod

結果となった．

org
modElapsed �me[sec]
3200
2800

4.49
3.62

SpeedupElapsed �me[sec]

7.07

2400

other

2000

matrix post

1600

matrix output
matrix loop2

1200

matrix mpi

800

matrix loop2

400

matrix mat_const
matrix pre

0

2.95
1.85
1.00

9p
(3x1x3)

1.75

72p_1t64p_1t

72 並列

64 並列

SystemC

18p(3x2x3)36p(6x1x6)72p(6x2x6)プロセス数

図 2: プロセス数を変えた場合の System C での実行
時間54 − Academic Center for Computing and Media Studies, Kyoto University

64p_4t32p_8t16p_16t

256 並列8p_32tvcell&src
input

SystemA

図 3: スレッド・プロセス数を変えた場合の System

A と C の実行時間の比較

4 積層構造を考慮した CFRP 中の
超音波伝搬解析チューニングされた FEM コードを用いて，CFRP
中を伝搬する超音波のシミュレーションを行った．こ超音波探触子

50mm
x3

x1

20mm

擬似等方積層材oooo[(45 /0 /90 /-45 )18]s

x2

50mm

こでは，疑似等方積層材の数値モデルを作成した．疑似等方積層材とは，同じ直交異方性板を弾性対称軸

20mmが等角となるように積層接着したもので，面内で巨視的に等方性を持つようになる性質を有する．ここで模擬した疑似等方積層材は，プリプレグを 45◦ ，0◦ ，

90◦ ，−45◦ 方向に回転したものを規則的に積層して1.0 sおり，144 層からなる．数値解析では，層内は均質材

2.0 s

料とし，各層に炭素繊維の配向方向を考慮した弾性ス

20mm

ティフネスを設定した．モデルの寸法は一辺が 20mmとし，1 つのボクセルの一辺を 0.02 mm としたため，総ボクセル数は 10003 となった．直径が 6.3mm の探触子をモデルの上部に設置し，ここから中心周波数
が 2.25MHz の 1 波を送信した．時間刻み ∆t は 1.0

x3

ns とし，15000 ステップの時間更新を行った．可視

3.0 s化結果は x2 方向の変位に色をつけたものである．

x1

x1 –x3 面の超音波の伝搬を可視化した結果を図 4(a)4.0 s(a)に示す．x1 –x3 面は等方面であるが，表層付近のプ

50mm

リプレグの配向の影響を受けて，超音波は斜め方向に強い振幅で伝搬していることがわかった．層間の散乱によって，波頭に続いて波が尾引いて伝搬している様子が観察できる．また，疑似等方積層材 CFRPの試験片 (50mm×20mm×50mm) を作成し，それをレーザー超音波法によって可視化した結果を図 4(b)x3
6.0 sx1に示す．数値シミュレーションは，実際の現象を良好に模擬できている．

10.0 s(b)

50mm

図 4: (a) FEM による CFRP 中の超音波伝搬シミュレーション，(b) レーザー超音波法による可視化計測．

5 まとめ

ドチューニングについて直接ご指導頂いたクレイ・

今回の共同研究によって，コンパイラによる最適

ジャパン・インクの鈴木 幸朗氏に感謝申し上げます．

化が可能となるようなコードチューニングや，MPI
通信の見直し，入力処理の最適化等についてご指導を頂いた．これによって，スケーラビリティが大幅に改善され，大規模計算時の効率が向上した．特に，入力処理の最適化は，ファイルの読み込み時間の短
縮だけでなくメモリの大幅削減も可能となり，大規模計算の実行に大いに貢献するものである．
謝辞： 有益なご助言を頂いた「スーパーコンピュー
タシステム共同研究企画委員会」の皆様，並びにコー

参考文献[1] K. Nakahata，J. Chang，M. Takahashi，K. Ohira and Y.
Ogura，Finite integration technique for coupled acoustic
and elastic wave simulation and its application to noncontact ultrasonic testing, Acoustical Science and Technology, Vol.35, No.5, pp.260-268, 2014.
[2] 中畑和之, 徳永淳一, 廣瀬壮一, イメージベース波動伝
搬シミュレーションと超音波探傷法のモデル化への応用, 非破壊検査, 59 巻, 5 号, pp.231-238, 2010.

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

超小型マイクロ波放電式中和器の電子引き出しを対象とした
3 次元プラズマ粒子計算
鷹尾 祥典**横浜国立大学 大学院工学研究院 システムの創生部門

1 はじめにかある。磁場中において電子は磁力線に巻き付いて運動を行うため、磁力線を横切る方向への移動

近年、大学や中小企業など比較的小規模な組織は比較的少ない。しかし、MIPS の中和器からはでも運用が可能な超小型衛星の研究開発が盛んに

磁場を横切って電子が引き出されており、この機行われている。超小型衛星は短期間・低コストで構を解明することにより中和器の効率向上を見込製作できるため、教育や技術実証、商用利用などめると考えられる。

幅広いミッションに用いられている[1]。しかし、実験では測定が困難な情報を入手できるという超小型衛星は重量や容積、電力等に厳しい制約が点で数値計算は強力な解析手段となり得る。これあるため、推進機を搭載した例は未だ少なく、搭までに、プラズマ生成から電子引き出しまでを統載する推進機は小型である必要がある [2]。衛星

一的に解析できる 3 次元プラズマ粒子計算モデルに推進機を搭載できれば、衛星自身による軌道選を構築し [4]、中和器からの電子引き出し機構の択や深宇宙探査が可能になるなど、より高度なミ解析を行っている [5, 6]。

ッションへの利用が期待できる。

現状、１つの条件当たり十分な時間平均を取っ

このような要求に対応する小型推進機の 1 つとて統計的に信頼できる解析結果を得るためには、してイオンスラスタが挙げられる。東京大学と最低でも 1 ヵ月以上（～数ヵ月）計算を続ける必NExt generation Space system Technology Research

要があり、多くの条件に応じた結果を得て、工学

Association (NESTRA) によって開発された 50 kg

的に有益な解析を行うためには計算の高速化が不級の超小型衛星 HODOYOSHI-4 には小型イオン可欠である。

推進システム (MIPS: Miniature Ion Propulsion
System) が搭載され、2014 年 10 月に宇宙空間で2 数値計算手法の運用に成功した。この MIPS にはイオン源と中和器にマイクロ波アンテナとリング状磁石を用い

図 1 に本研究が対象とする小型マイクロ波放電

た ECR（Electron Cyclotron Resonance）放電を採用

式中和器の計算領域を示す [7]。金属面とプラズしている [3]。

マ源の境界においてその中心を原点とする直交座中和器は放電室に生成したプラズマのうち電子標系を設定した。プラズマ源を 20×20×4 mm3 のを放出するが、電子はイオンに比べて質量が桁違領域とし、その内部（z = 1.0 mm）にはリング状アいに小さいため推力にほとんど寄与しない。そのンテナが配置されている。また、プラズマ源の下ため中和器から効率的に電子を引き出すことが必流側（z = 4.0 mm）には 4 つの穴を持つ厚さ 0.6 mm

要であり、それによりイオンスラスタ全体の比推のオリフィス面を備えており、この面から z = 10

力向上につながる。また、MIPS の宇宙空間でのmm までを電子引き出しの解析領域とした。オリ作動は実証されたが、中和器からの電子引き出し

フィス面を含む金属面のポテンシャルはゼロとし、機構には未だ解明されていない物理現象がいくつ

電子引き出しのために z = 10 mm でのポテンシャ

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

図 2: 計算フローチャート
スウェル分布で、時間的・空間的に一様とする。
(iii) 荷電粒子と中性粒子の衝突は、電子の弾性・励起・電離の各衝突およびイオンの弾性・電荷交換の各衝突を考慮する。(iv) マイクロ波による磁場は、永久磁石による静磁場と比べて小さいため無視する。(v) 本研究ではマイクロ波の電力が小さいためにプラズマ電流の影響を無視する [8]。
図 2 に本計算のフローチャートを示す。はじめに初期条件を与え、定常状態におけるマイクロ波に よ る 電 磁 場 を 得 る た め に Finite-difference
Time-domain (FDTD) 法によりマクスウェル方程式を解く。このとき、時間ステップはtEM =1.49×10−13 s（4.2 GHz マイクロ波の 1/1600 周期）である。次に、PIC/MCC 法により粒子の運動を、
時間ステップte = 5.95×10−12 s（マイクロ波の 1/40
周期）の下で解く。
マイクロ波電場の時間変化は、先ほど求めた定常状態のマイクロ波電場 E0 を用
いて、EEM = E0 cos(t) で与えられる。また、永久図 1: 解析対象とする電子源の構造 [7](JSASS-2016-4687 Ⓒ日本航空宇宙学会)
ルは 20 V とした。なお、プログラム高度化前の計磁石による静磁場Bst はANSYS を用いて解析している（図 3）。マイクロ波周波数 4.2 GHz はイオン
10算においては、計算時間短縮のため図の対称性を|B|(T)
5

利用し xy 平面のうち第一象限（1/4 領域）だけ計

0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05

本研究で用いる粒子計算法は Particle-in-Cell(PIC)法であり、荷電粒子と中性粒子との衝突をMonte Carlo collision (MCC) 法により考慮している。モデルの詳細は過去の文献に委ね [4-6]、ここy(mm)算している。
ECR

0

layer

-5

では簡単に述べるに留める。計算セルは各方向とも 0.1 mm の等間隔格子であり、計算コスト削減のため、1/4 領域による対称計算を以下の仮定の下で行った。(i) 1 価のキセノンイオンと電子のみを粒子として扱う。(ii) 中性粒子は 300 K のマク

-10
0

5

z(mm)

10

図 3: 磁束密度分布 [7](JSASS-2016-4687 Ⓒ日本航空宇宙学会)

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

プラズマ周波数と比べて非常に大きいので、計算時間短縮のためイオンの運動では EEM を無視し、(a) 1/4 領域

静電場 EES はマイクロ波 1 周期で平均化した値を用いて、時間ステップti = 2.38×10−10 s（マイクロ波の 1 周期）で更新している。本研究ではプラズ
マ吸収電力 P0 を入力パラメータとしており、イオ
ンと電子の運動方程式を解く前後での粒子の運動
エネルギーの差から得られる電力Pabs がP0 を満たすように E0 を調節している。上記の計算を定常解(b) 全領域が得られるまで繰り返し行う。
計算条件は、中性ガス圧力 p0 = 1 mTorr、マイク
ロ波周波数 f = 4.2 GHz、吸収電力 P0 = 0.3 W であり、推進剤としてキセノンを用いている。また、
電子とイオンの初期温度をそれぞれ 2.0 eV、0.05
eV とし、いずれも初期密度 1×1016 m−3 で一様に分布させた。

3 プログラム高度化内容と性能評価
図 2 のフローチャートが示すように、本プログ
ラムで大部分の計算時間を占めるのは個々のプラ
ズマ粒子の運動変化を求めるループとなる。
プロファイル計測の結果によると、中でもポア

図 4: 高度化前後におけるプラズマ密度分布：(a) 1/4 領域（高度化前）、(b) 全領域（高度化後）

ソン方程式を解いてプラズマポテンシャルを求めるサブルーチンが実に全体の90%を占めている事

必要なく並列化することが可能になる。次に、ルが明らかになった。これは高度化前のプログラム

ープ内の依存関係を回避するために RedBlack 法

において、図 1 の計算領域が示すようにプラズマを採用することにした。また、現状扱っている粒が分布する領域にアンテナやオリフィスなどが存

子数から 1000 倍以上粒子を増やさない限り、
シス在するため、遅くても最も単純でバグが潜みにく

テム B の 1 ノードで実行可能なレベルであり、かいことを優先し SOR 法を採用した事に起因してつ、ポアソン方程式を解くサブルーチン内でのメいる。また、最内ループに IF 文を用いて計算領域

モリ使用量も 100 MB 程度と小さく、格子間隔や内のプラズマの有無を判断している点も要因とし粒子数に依存しないことからプロセス並列ではなて挙げられる。

くスレッド並列を採用した。

ポアソン方程式の高速化において、SOR 法ではこれらの対策の結果、ポアソン方程式を解くサ

配列の FLOW 依存（前方依存）、ANTI 依存（後ブルーチンにおいて、オリジナル（1n1p1t）に対方依存）が存在するためコンパイラによるベクト

してスレッド並列化（1n1p36t）の適用により、43.7

ル化およびスレッド並列が効かない状況である。

倍の高速化が実現できた。その他、細かいチュー

この状況を改善するため、まずは最内ループ内にニングを施すことで最終的には約 53 倍の高速化存在する IF 文において TRUE と FALSE のいずれを達成している。

の場合でも演算式を共通化するマスク方式を採用

なお、SOR 法とは演算順序が異なるため、計算して IF 文を除去した。これにより、本来ポアソン

結果の完全一致は得られないが、両解法における方程式を解かずに済むアンテナやオリフィスプレ

計算結果の差異は実用上問題ないレベルであったート内部についてもポアソン方程式を解く演算を

（図 4）
。

行う必要が生じるが、複雑な境界条件を判定する

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

ポアソン方程式を解くサブルーチンの高速化が実現出来た結果、全計算時間に占める割合は 50％[2] Micci, M. M. and Ketsdever, A. D.: Micropropulsion for程度まで削減できた。その結果、計算格子内のプ

Small Spacecraft, American Institute of Aeronautics andラズマ粒子の位置から計算格子の各コーナーのグ

Astronautics, Reston, 2000.

リッドへ電荷などを割り当てる allocation のサブ[3] Koizumi, H., Komurasaki, K., Aoyama, J., and Yamaguchi,

ルーチンが計算時間全体に占める割合が大幅に増

K.: Engineering Model of the Minature Ion Propulsion

加し 30%程度を占めるに至った。この allocation

System for the Nano-satellite: HODOYOSHI-4, Trans.

の計算時間削減のため、粒子番号によるデータ並

JSASS Aerospace Tech. Japan, 12 (2014), pp. Tb_19–

列化、omp atomic によるループ内排他処理、デー

Tb_24.

タ移動量削減、スレッド同期削減など施すことにより 15%まで低下した。

[4] Takao, Y., Koizumi, H., Komurasaki, K., Eriguchi, K., andOno, K.: Three-dimensional particle-in-cell simulation of a

その他のサブルーチンについても細かいチュー

miniature plasma source for a microwave discharge ion

ニングを行った結果、プログラム全体としては、

thruster, Plasma Sources Sci. Technol., 23 (2014), pp.

オリジナルのコードに対して約 36 倍まで高速化

064004-1–11.
[5] Takao, Y., Hiramoto, K., Nakagawa, Y., Kasagi, Y.,

することができた。
今後この高度化されたコードを利用して、電子Koizumi, H., and Komurasaki, K.: Electron extraction

引き出し効率の向上を目指した解析を進めていく

mechanisms of a micro-ECR neutralizer, Jpn. J. Appl.

予定である。

Phys., 55 (2016), pp. 07LD09-1–5.
[6] Hiramoto, K., Nakagawa, Y., Koizumi, H., and Takao, Y.:

謝辞

Effects of E × B drift on electron transport across theプログラム高度化の実施にあたって、クレイ・ジ

magnetic field in a miniature microwave dischargeャパン・インクの伊藤 泰善 様には大変ご尽力いneutralizer, Phys. Plasmas, 24 (2017), pp. 064504-1–5.

ただきました。また本研究は、JSPS 科研費[7] 平本謙太, 中川悠一, 小泉宏之, 小紫公也, 鷹尾祥

JP16H06370 の助成を受けて実施されました。こ典：小型マイクロ波放電式中和器のプラズマ源におこに感謝の意を表します。

ける E×B ドリフトの効果、第 58 回宇宙科学技術連
合講演会講演集, 2016, JSASS-2016-4687。
[8] Takao, Y., Eriguchi, K., Ono, K., Sugita, Y., Koizumi, H.,

参考文献

and Komurasaki, K.: A Validation Study of a 3D PIC[1] Richardson, G., Schmitt, K., Covert, M., and Rogers, C.:

Model for a Miniature Microwave Discharge Ion Thruster,

Small Satellite Trends 2009–2013, 29th Annual

50th AIAA/ASME/SAE/ASEE Joint Propulsion Conf.,

AIAA/USU Conference on Small Satellite, 2015,2014, AIAA-2014-3829.

SSC15-VII-3.

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

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

サービス再開日時

障害発生日時保守時間[h]2017/10/09

6:30

2017/10/11

9:30

51.00

2017/12/07

9:00

2017/12/08

9:30

24.50

2018/03/28

9:00

2018/04/01

0:00

87.00

2017/12/06

サービス再開日時
2:55

2017/12/06

ダウン時間[h]3:00２）サービス状況
ジョブ

サービ
ス時間
[h]

処理
件数

10 月

693.00

5,498

11 月

720.00

12 月

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

ノード
利用率

52,01824,241,191

21,717,609

1761.4

74

%

1,703

7,419

5,624,335

3,170,218

1798.578%719.42

22,380148,146

32,674,482

36,686,266

1793.677%

1月

744.0036,954

332,976

16,783,488

68,650,710

1799.478%

2月

672.0018,067

179,926

14,510,299

55,890,977

1799.574%

3月

657.0015,479130,732

13,414,253

55,378,805

1759.4

72

%

計

4,205

100,081851,127

107,248,048

241,494,585

1785.3

75

%

経過
時間[h]

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



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

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



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

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

0.08

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

サービス再開日時

システムダウン障害発生状況
障害発生日時保守時間[h]サービス再開日時

ダウン時間[h]2017/10/09

6:30

2017/10/11

9:30

51.00

2017/10/16

14:50

2017/10/16-

17:10

2.35-

2017/10/13

15:00

2017/10/13

16:10

1.17

2017/11/08

17:00

2017/11/08

18:40

1.67

2017/12/07

9:00

2017/12/08

9:30

24.50

2017/12/29

12:48

2017/12/31

8:55

44.12

2018/03/28

9:00

2018/04/01

0:00

87.00

2018/01/02

0:00

2018/01/04

8:55

56.92

2018/01/22

15:01

2018/01/22

15:48

0.78

2018/03/17

21:02

2018/03/19

10:30

37.47

２）サービス状況
ジョブ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

689.5080,717557,722

10,601,897

9,498,841

805.5

70

%

11 月

718.3351,629593,214

11,176,725

10,205,141

831.3

75

%

12 月

675.3852,790676,268

12,406,912

12,013,002

825.7

71

%

1月

686.30101,998

630,024

13,276,605

12,284,851

830.5

73

%

2月

672.0071,048498,405

11,932,318

10,970,167

831.9

70

%

3月

657.0039,016

351,428

9,364,962

8,582,865

821.0

58

%

計

4083.39

397,198

3,307,061

68,759,419

63,554,867

824.32

70

%

ノード
利用率

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



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

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



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

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

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

サービス再開日時

システムダウン障害発生状況
障害発生日時保守時間[h]サービス再開日時

ダウン時間[h]2017/10/09

6:30

2017/10/11

9:30

51.00

2017/10/16

14:50

2017/10/16-

17:10

2.35-

2017/10/13

15:00

2017/10/13

16:10

1.17

2017/11/08

17:00

2017/11/08

18:40

1.67

2017/12/07

9:00

2017/12/08

9:30

24.50

2018/01/22

15:01

2018/01/22

15:48

0.78

2018/03/28

9:00

2018/04/01

0:00

87.00

2018/03/17

21:02

2018/03/19

10:30

37.47

２）サービス状況
バッチ

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

689.58

1,961

10,564

156,433

104,929

15.3

33

%

11 月

718.33

2,549

12,812

136,407

67,505

16.0

43

%

12 月

719.502,980

19,153

277,784

135,491

15.8

57

%

1月

743.22

6,949

26,638379,283

175,912

16.0

61

%

2月

672.002,307

19,824

197,546

126,036

16.0

41

%

3月

657.001,490

31,051

139,207

116,192

15.8

32

%

計

4184.43

18,236

120,042

1,286,660

726,065

15.8

45

%

ノード
利用率

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



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

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



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

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

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

サービス再開日時

障害発生日時保守時間[h]2017/10/09

6:30

2017/10/11

9:30

51.00

2017/12/07

9:00

2017/12/08

9:30

24.50

2018/03/28

9:00

2018/04/01

0:00

87.00

サービス再開日時-なし

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

サービ
ス時間
[h]

処理
件数

経過
時間[h]

占有
時間[h]

CPU
時間[h]

平均稼動
ノード数

10 月

693.00284

1,540

661,982

1,752

482.0

27

%

11 月

720.00163

2,742

1,022,241

15,553

481.7

28

%

12 月

719.50167

1,222

707,103

5,191

481.9

21

%

1月

744.0040

2,423

1,263,255

3,681

481.6

34

%

2月

672.0086

1,211

442,478

3,072

482.0

16

%

3月

657.0057

3,636

1,069,615

3,815

482.0

28

%

計

4196.25797

12,774

5,166,674

33,064481.87

26

%

ノード
利用率

ノード利用率
100%
90%
80%
70%
60%
50%
40%
30%
20%
10%
0%

10月

11月

12月

1月

2月

3月



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

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



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

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

資料

別表１ スーパーコンピュータシステム
提供サービス
利用負担額
コース

タイプ

セット

エントリ

タイプA
タイプB
タイプC

基本
基本
基本
基本
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位
最小
追加単位

パーソナル

タイプA1
タイプA2
タイプA3
タイプB1
グループ
タイプB2
タイプB3
タイプC1
タイプC2
タイプA
大規模ジョブ タイプB
タイプC
専用クラスタ-ストレージ容量追加
ライセンスサービス

システム バッチ

12,600 円/年
100,000 円/年
100,000 円/年
100,000 円/年
200,000 円/年
100,000 円/年
240,000 円/年
60,000 円/年
600,000 円/年
300,000 円/年
210,000 円/年
105,000 円/年
252,000 円/年
63,000 円/年
630,000 円/年
315,000 円/年
130,000 円/年
130,000 円/年
156,000 円/年
78,000 円/年
20,000 円/週(7日）
10,000 円/週(7日）
21,000 円/週(7日）
10,500 円/週(7日）
13,000 円/週(7日）
6,500 円/週(7日）
630,000 円/年
315,000 円/年
10,000 円/年
20,000 円/年

システム資源

経過時間 ストレージ
無料
（時間）
（TB)
利用者数

最大1ノード相当((36コア、128GBメモリ)×1)
1
最大4ノード相当((68コア、16+96GBメモリ)×4)
168
最大4ノード相当((36コア、128GBメモリ)×4)
最大1ノード相当((72コア、3072GBメモリ)×1)
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)
336
8ノード((36コア、128GBメモリ)×8)
B
準優先
2ノード((36コア、128GBメモリ)×2)
8ノード((36コア、128GBメモリ)×8)
占有
4ノード((36コア、128GBメモリ)×4)
1ノード((72コア、3072GBメモリ)×1)
優先
1ノード((72コア、3072GBメモリ)×1)
C
2ノード((72コア、3072GBメモリ)×2)
準優先
1ノード((72コア、3072GBメモリ)×1)
8ノード((68コア、16+96GBメモリ)×8)A
占有
4ノード((68コア、16+96GBメモリ)×4)
8ノード((36コア、128GBメモリ)×8)B
占有
168
4ノード((36コア、128GBメモリ)×4)
2ノード((72コア、3072GBメモリ)×2)
C
占有
1ノード((72コア、3072GBメモリ)×1)
8ノード((36コア、128GBメモリ)×8)
B
4ノード((36コア、128GBメモリ)×4)
ストレージ容量10TBの追加につき可視化ソフト（AVS,ENVI/IDL)およびプリポストウェアの1ライセンスにつき
B
A
B
C

共有
共有
共有
共有

0.2
3.0
3.0
3.0
24.0
12.0
28.8
7.2
48.0
24.0
24.0
12.0
28.8
7.2
48.0
24.0
24.0
24.0
28.8
14.4

-

48.0
24.0

8
4
16
4
16
8
8
4
16
4
16
8
8
8
16
8

-

16
8

備考
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 ． スパコン連携サービス
学術情報メディアセンターのスーパーコンピュータシステムと密な連携により、学内における部局の組織が計算サーバ等を設置する
場合、下記の負担額を支払うものとする。
冷却方式
水冷
空冷

利用負担額
9,800 円/月
11,500 円/月

利用負担額算定単位
水冷冷却方式の計算サーバ等の定格電力 1kWにつき
空冷冷却方式の計算サーバ等の定格電力 1kWにつき

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

別表２ 汎用コンピュータシステム
区

分

利用負担額

仮想サーバホスティングサービス

単位

36,000円／年 1仮想サーバにつき

備考
1． 利用負担額は、総額表示である。
2． 上記表の仮想サーバホスティングサービスを利用するには、スーパーコンピュータシステムの利用者であること。
3． 1仮想サーバに割当てるシステム資源は、CPU：2コア、メモリ：4GB、ディスク：100GBである。
4． 仮想サーバホスティングサービスにおいて、下記の負担額を支払うことによりCPU、メモリ、ディスクを増量することができる。
区

分

利用負担額

単

CPU増量

3,000円／年 2コアにつき(最大8コアまで)

メモリ増量

3,000円／年 4GBにつき(最大64GBまで)

ディスク増量6,000円／年 100GBにつき(最大1,000GBまで)位

5． 利用負担額は、当該年度(4月から翌年3月まで)の利用に対して年額として算定するが、年度途中から利用を開始する場合には月数に応じて減額する。

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

※民間機関利用

別表３ スーパーコンピュータシステム
システム

A

B

C

システム資源

経過時間 ストレージ
無料
（時間）
（TB)
利用者数

利用負担額

8ノード(68コア、16+96GBメモリ)×8)

336

28.8

16

960,000 円/年

12ノード(68コア、16+96GBメモリ)×12)

336

43.2

24

1,440,000 円/年

16ノード(68コア、16+96GBメモリ)×16)

336

57.6

32

1,920,000 円/年

8ノード(36コア、128GBメモリ)×8)

336

28.8

16

1,008,000 円/年

12ノード(36コア、128GBメモリ)×12)

336

43.2

24

1,512,000 円/年

16ノード(36コア、128GBメモリ)×16)

336

57.6

32

2,016,000 円/年

2ノード(72コア、3072GBメモリ)×2)

336

28.8

16

624,000 円/年

3ノード(72コア、3072GBメモリ)×3)

336

43.2

24

936,000 円/年

4ノード(72コア、3072GBメモリ)×4)

336

57.6

32

1,248,000 円/年

備考
1． 利用負担額は、年度単位で算定している。また、総額表示である。パーソナルコース、グループコース又は専用クラスタコース
年度途中から利用を開始する場合及び年度途中で利用を終了する場合の利用負担額は、上記表中の利用負担額を12で除し利用月数を乗じて算出するものとし、100円未満に端数が出た場合は、10円単位を四捨五入するものとする。
なお、月途中から利用を開始する場合及び月途中で利用を終了する場合は、それぞれ1月の利用とする。
2． ストレージ容量はバックアップ領域（最大で総容量の1/2）を含む。

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

全国共同利用版広報・Vol.16(2017)総目次［巻頭言］
Vol.16，No.1 号の発刊にあたって ................................................................................................................1-1
Vol.16，No.2 号の発刊にあたって ................................................................................................................2-1［スーパーコンピュータ共同研究制度（若手研究者奨励枠）研究報告］
種々の波動問題に対する積分核に依存しない演算子積分時間領域高速多重極境界要素法の開発 .........1-2
光機能性有機分子材料の構造と電子構造の解明による有機系太陽電池の素子性能向上.........................1-4
乱流相分離制御による自己組織化構造の能動的選択 ..................................................................................1-6
Singlet Fission 相互作用の量子化学計算 .....................................................................................................1-9
高解像度大気海洋波浪結合モデルを用いた海面物理モデルの現地観測結果に基づく検証...................1-11
iPS 細胞懸濁培養の最適化を目的とした撹拌槽内粒子挙動解析..............................................................1-13
大規模津波による流出油の拡散挙動予測....................................................................................................1-15
視聴覚モダリティ間デコーディングによる感覚間協応のメカニズムの解明..........................................1-17
高次元準周期軌道における準周期軌道の解析 ............................................................................................1-19
メモリ一貫性モデルを考慮したモデル検査器の開発 ................................................................................1-21
固液混相流の数値粘性測定...........................................................................................................................1-23
高プラントル数流体を用いた HZ 液柱内温度差マランゴニ対流の数値解析―特異な対流場における熱的相互作用の影響― ........................................................................................1-25
カオス力学系におけるネットワーク解析....................................................................................................1-28
Numerical investigation of the effect of hydrogel surface shape on drag performancein turbulent flow ...........................................................................................................................................1-30
MuSTAR MD に基づくタンパク質構造変化の自由エネルギー経路計算手法の確立と応用 ................1-34
破壊力学に基づく損傷モデルによる鉄筋コンクリートの 3 次元破壊シミュレーション ......................1-36
分子性導電・磁性材料の基底状態および励起状態に関する理論的解析..................................................1-38
ポスト「京」100 ナノ電子状態計算むけの超並列構造保存波束ダイナミクスソルバー開発 ...............1-40
固相光化学反応の計算化学による解明........................................................................................................1-42
漂流物の衝突力に関するばらつきの基礎的検討 ........................................................................................1-44
気象条件に応じた都市の境界層乱流の組織構造に関する研究.................................................................1-47
海色衛星データを用いた沿岸および陸棚海域における河川由来物質の解析手法の開発.......................1-50
硫酸ランタン La_2(SO_4)_3 結晶中における水拡散チャネルの発見 .....................................................1-52
粒子との接触を伴う液体挙動の直接数値解析～自由界面―固体物体間相互作用を考慮した多相直接数値解析 .............................................................1-54
高精度視覚質感記憶の心理学的基盤と神経機構の解明.............................................................................1-57
臨界レイノルズ数における球の流体力に対する周囲の乱流場の影響 .....................................................1-59
3 次元ナビエストークス流の過渡的乱れ ....................................................................................................1-61［プログラム高度化支援事業研究報告］
動的/静的水～土骨格連成有限変形解析コードの高度化～固有値解析に関わる高速化、並列化の検討～ ........................................................................................1-62
津波による瓦礫浮遊物を考慮した構造流体連成解析 ................................................................................1-66

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

［バーストバッファによる高速ファイル I/O］

バーストバッファの利用について～ IME 編～...........................................................................................2-3
流体系シミュレーションでのプロセス分散書き出しにおける IME の効果 .............................................2-8［性能評価報告］
新スーパーコンピュータのベンチマーク評価報告 ....................................................................................2-11［サービスの記録・報告］
スーパーコンピュータシステムの稼働状況.................................................................................... 1-70，2-18
汎用コンピュータシステムのサービス状況................................................................................................1-75
センター利用による研究成果（平成 28 年度）..........................................................................................2-22［資料］
大型計算機システム利用負担金 別表...............................................................................................1-76,2-26
全国共同利用版広報・Vol.15(2016)総目次 .................................................................................................1-79
サービス利用のための資料一覧 ..........................................................................................................1-81,2-29［編集後記］
編集後記、奥付 .....................................................................................................................................1-82,2-30

68 − 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※ ホストへの接続は SSH(Secure SHell) 鍵認証のみ、パスワード認証は不可

2.

問い合わせ先 & リンク集


情報環境機構のホームページ
http://www.iimc.kyoto-u.ac.jp/



学術情報メディアセンターのホームページ
http://www.media.kyoto-u.ac.jp/



スーパーコンピュータシステムに関する問い合わせ先

利用申請などに関する問い合わせ先【情報環境支援センター】
E-mail : zenkoku-kyo@media.kyoto-u.ac.jp / Tel : 075-753-7424
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/

システムの利用など技術的な問い合わせ先【スーパーコンピューティング掛】
E-mail : consult@kudpc.kyoto-u.ac.jp / Tel : 075-753-7426
URL: http://www.iimc.kyoto-u.ac.jp/ja/services/comp/contact.html

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

奥付京都大学学術情報メディアセンター全国共同利用版広報 Vol．17，No．1
２
０１８年
１４日xx
発行
2018
年 １２月
12 月
日発行
広報編集部会
編集者 京都大学学術情報メディアセンター

深沢 圭一郎（部会長）

全国共同利用版広報編集部会

平石 拓（副部会長）

発行者 〒606-8501 京都市左京区吉田本町
京都大学学術情報メディアセンター

水谷 幸弘

小林 正信

熊谷 真由美尾形 幸亮

Academic Center for Computing and Media Studies
Kyoto University
Tel. 075-753-7407
http://www.media.kyoto-u.ac.jp/
印刷所 〒616-8102 京都市右京区太秦森ヶ東町 21-10
株式会社エヌジーピー

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

表紙デザイン：中山 豊（中山商店）