局所サンプリング型 MD／CFD ハイブリットシミュレーション
安田修悟、山本量一
京都大学工学研究科化学工学専攻
yasuda@cheme.kyoto-u.ac.jp
概要： ソフトマターの流動解析において有効な流体力学計算と分子動力学計算のハイブリット計算法について発表する。 ソフトマターの流動解析においてはその局所的な応力特性をモデル化する構成方程式を求めることは一般に困難である。我々のハイブリット法では流体力学計算において必要とされる局所的な応力に対して、構成方程式を用いることなく、各局所的な点に貼り付けられた分子動力学セルにおいて直接応力のサンプリングをする。発表では、局所サンプリング法を用いて行った高分子溶液の流動解析など幾つかの応用例について紹介する。

1 はじめにリングによって解決することを考える。

我々の身近にある製品の多くは、柔らかくて変形しやすい物質、ソフトマター（コロイド、高分子、液晶など）からできている。例えば、化粧品、日焼け止め、ヨーグルト、ゼリー、液晶テレビなどはこれらを使った代表的な製品である。また、最近のナノテクノロジーの出現によって、ソフトマターは特殊な機能を持つ新しい材料の宝庫としても注目されており、医療、エネルギー、電機、情報など様々な分野で今後益々ソフトマターを利用した機能性材料が開発されることが期待される。
計算機シミュレーションは、ソフトマター研究における重要なアプローチの一つである。ソフトマ
ターのシミュレーション研究の難しい点は、原子・分子レベルでの構造と巨視的なレベルでの運動との強い相関にある。例えば、実験で使う試験管（数センチ程度）に入ったコロイド溶液の流れをシミュレーションしようとすると、それより遥かに小さなコロイド粒子（数マイクロ程度）の配置や運動も考慮しなくてはいけなくなる。電荷を持ったコロイドを考える場合にはさらに小さなイオン（数ナノ程度）の運動まで考える必要が生じる。ここで問題になるのが、途方もない計算時間である。イオンの運動に目を向けた場合、試験管の中の流れはそれより 109 倍も大きなスケールになり、それを計算することは現在の世界最速の計算機もってしても我々が生きている間には終わらない。それだけでなく、この計算が現実的な計算時間で完了するようになるには、コンピューターの性能向上がこの先も全く鈍ることなく継続するという最も楽観的なシナリオでさえ、これから 100年は待たなければならない。
我々は、ソフトマターのシミュレーションが抱えるこのような難題に対して、マルチスケールモデ

2 局所サンプリング法
本稿では最近の成果である、高速で振動する平板間における高分子溶液の運動をマルチスケールシミュレーションによって解析した結果[1]を紹介し我々のマルチスケールモデリングの概要を手短に説明する。[2]
図 1 に本研究でシミュレーションを行った系の略図を示す。平板間の距離 2H は高分子を構成する一つの粒子の直径σの数千倍程度である。振動板は高速に振動しており、粒子のサイズσをナノメートル程度と考えた場合、振動板の一周期はおよそがナノ秒程度である。高速に振動する平板上では粘性拡散層が生じることが知られている。粘性拡散層は流体力学のスケールでは薄層として取扱われるが、分子の直径に比べると通常はるかに大きなスケールを持つ。従って、境界層内部での流体の振舞を解析することは通常の分子動力学（MD）では困難な問題の一つである。
高分子溶液の巨視的な流れは、連続体モデルの運動量保存式によって計算される。ここで、各計算点での x 方向の流速 vx はその点での局所的なせん断応力σxy の y 軸方向の勾配によって変化する。通常のニュートン流体では、応力と流速の空間勾配の間に比例関係、σxy=η∂vx/∂y, があることが知られているので、その比例定数、即ち粘性係数ηさえあらかじめ決めておけば、各点で求めた流速から各点での応力を求め、その瞬間の流速変化を計算することができる。しかし、高分子溶液にはそのような単純な比例関係は存在せず、その一般な構成方程式は求められていない。
我々のマルチスケールモデリングでは、巨視的な流れを計算する各計算点に分子レベルでの高分子の運動を計算するための小さな MD セルを貼り付けておき、そこで局所的な流速場と過去の履歴に基づく局所的な応力を直接分子動力学を用いて計算する。
（局所応力サンプリング）MD セルには流体計算で得られた局所的なひずみ速度の情報を与え、Leeds-Edwards の周期境界条件を用いて MDセルにせん断流れを作り、その応力をサンプリングする。その際、MD 計算は流体計算で用いられる時間ステップに対応する時間だけ行われ、各MD セルで得られる高分子の配置・構造は、そのまま次の流体計算の時間ステップにおける MD 計算の初期配置としてメモリーされる。この方法を用いることで、高分子溶液に特徴的な粘弾性や履歴効果を正確にシミュレーションすることができるようになる。
アルゴリズム上の特徴としては、流体の計算に用いるメッシュ間隔に比べて小さな MD セルを用いることで、系の全領域を全て MD で計算した場合に比べて圧倒的に少ない計算量で巨視的な流れを計算することができる。
（具体的には、各ＭＤセルには 10 粒子を繋げた高分子鎖が 100 本入っており、流体計算の各時間ステップにおいてＭＤ計算は 1000 ステップ実行する。）さらに、特筆すべき点には、計算時間の大部分を占める MD 計算部はそれぞれのセルで独立に計算を行うことができる
ため、各 MD セルに対してそれぞれ一つのＣＰＵを割り当てることによって、数千以上のＣＰＵを用いる超並列計算においてもほぼ理想的な並列化効率を容易に達成することができる点にある。

4 次世代スパコンへの展開
我々のマルチスケールシミュレーション法は、通常の分子動力学法とは異なり、次世代スパコンのように膨大な CPU を持つコンピューターとの相性が極めて高いというユニークな特徴も持っている。これは、計算負荷の最も高い分子動力学法による計算部分が各計算点で独立であるため、計算の並列実行が容易であることに起因する。この方法を次世代スパコンで用いることにより、分子数が 100 億個にも達する世界最大規模の分子動力学シミュレーションを実現することができると考えられる。これまで有効なシミュレーションを行うことが難しかった、大規模な流動・変形現象を伴う機能性材料の製造プロセスの解析といった、理論や実験では全く歯が立たない未開拓の問題へのアプローチが可能となることが期待される。
参考文献[1] S. Yasuda and R. Yamamoto, EPL 86, 18002(2009). [2] S. Yasuda and R. Yamamoto, Phys.
Fluids 20, 113101 (2008).

図１：問題と計算系の概略図。

3 計算結果と並列性能
平板間距離を 3200σとしそれを64スリットに分割して流体計算を行った。図２に高速振動平板上での高分子溶液の流速分布と対応するニュートン流体の流速分布の比較を示す。ニュートン流体は高分子溶液の線形領域での運動粘性係数を用いて計算した結果である。流速分布が明らかに異なるのがわかる。これは高分子溶液に特有なシア・シニングによる。即ち、振動板近傍の激しく動いている速度勾配の大きな領域では流体の粘性が小さくなり、振動板の運動は流体内部へ浸透しにくくスリップしているような状況になっているのである。
局所的なレオロジー特性についても大変興味深い現象があることが明らかにされた。これについては、ポスターで発表する。
並列性能のベンチマークとして、平板間距離を12800σとし平板間を 256 スリットに分割してハイブリット計算を 1000 ステップ行った結果を表 1 に示す。128CPU まではほぼ理想的な並列計算が実行されるのが確認される。

図２：高速振動平板上での流速分布。各色は異なる時間での流速分布を表す。
表 1：並列計算ベンチマーク。並列化効率は 32CPU を用いた並列計算の結果に対する効率を表す。
CPU

計算時間[分]並列化効率

32277.1―64128139.1
69.6

99.6 %
99.5 %