微粒子・溶媒・熱揺らぎの直接数値シミュレーション
岩下 拓哉 a,† , 山本 剛紀 a , 名嘉山 祥也 b , 金 鋼 c , 山本 量一 a,d ,a

京都大学大学院工学研究科 化学工学専攻，615-8510 京都市西京区京都大学桂bc

九州大学大学院 工学研究院 化学工学部門，819-0395 福岡市西区元岡 744

分子科学研究所 計算分子科学研究系，444-8585 愛知県岡崎市明大寺町字西郷中 38 番地
d

科学技術振興機構 CREST, 332-0012 埼玉県川口市本町 4-1-8†

iwashita@cheme.kyoto-u.ac.jp
概 要

コロイドなどの微粒子の分散系の流動は，微粒子と溶媒流動の相互作用によって複雑な様相を示す．微粒子分散系の流動挙動を理解・予測するためには，微粒子と溶媒流動の連成問題，さらに溶媒分子による熱ゆらぎの効果を同時に考慮したシミュレーションが必要である．我々は，微粒子・溶媒・熱揺らぎを連成する直接数値シミュレーションの方法 (Smoothed Proﬁle 法) の開発を行い，また大規模計算のために Smoothed
Proﬁle 法の並列化についても取り組んでいる．本研究では，微粒子分散系の特異なレオロジーおよび沈降過程について検討した．これらの系では微粒子の集団的なゆらぎとシステムサイズの競合が重要となる場合があるため大規模な直接計算が必要である．本研究により微粒子・溶媒・熱揺らぎおよび微粒子集団ゆらぎと多スケールにわたる階層構造と流動の様相が明らかにされると考える．

1

背景

微粒子分散系は，微粒子が液体に分散している系である．これは食品・塗料・顔料・化粧品・スラリーなどに見出すことができ，その用途および製造・加工プロセスにおいて流動挙動の把握が重要である．微粒子のサイズがミクロンオーダー以下程度の場合はコロイドと呼ばれる．コロイド分散系では熱揺ぎの効果が重要となる．
微粒子分散系の流動挙動は粒子間相互作用や熱揺らぎなどに影響され大変複雑な様相を示す．とりわけ溶媒の流動により媒介される流体力学相互作用は、長距離的かつ多体的であるため多粒子の時空間構造に決定的な役割を果たす．したがって微粒子分散系の流動の理解のためには，粒子の運動と溶媒の流動を連成して解くことが本質的に重要である．またコロイド系における熱ゆらぎは溶媒分子に起因するものであるから，溶媒を陽に考慮した連成計算が必要である．我々は流体力学相互作用を正確にかつ効率的に解ける Smoothed Proﬁle 法 (SPM) の開発を行っている.
現在使っている計算機能力（100GFlops 級）では，粒子の直径の 10〜20 倍程度の系のシミュレーションが限度である. これで解決できる問題は，均一で等方的な系を仮定したものに限られる．粒子沈降や凝集など，粒子の分布が空間的に不均一になり，その効果を考慮することが本質的である問題については，粒子の直径の 50 倍以上の系のシミュレーションが求められ，その実現にはペタ Flops級の計算機が必要になる．また，これぐらい大きな系のシミュレーションを実現することができれば，現実的な外部境界条件を設定するとこが可能になり，微粒子分散系に関する理学的・工学的な種々の問題に対して，初めて直接シミュレーションで答えを出すことができるようになる.

2計算手法

微粒子多体系の挙動を解くには，溶媒分子のスケールでなく微粒子スケールの階層の自由度を解くことが必要である．SPM では，溶媒の自由度について連続体的な取扱いを行う．有限の体積をもった微粒子は陽に扱う．微粒子多体系の問題に本質的となる微粒子と溶媒の相互作用
Rは，F HdS i · ﬀ, を解くことで考慮される．ここでi =R
ﬀ は分散媒の応力であり， dS i は i 番目の粒子について表面積分を表す．有限の格子を用いる数値計算では，時々刻々この表面積分を正しく計算することは至難であったが， SPM では微粒子-溶媒界面に滑らかなプロファイルを陽に導入して，(i) 物理的に正確で，(ii) 効率的な，相互作用の計算法を実現している．これにより多数の微粒子と溶媒の流動を直接シミュレーションすることが可能となっている．SPM では，溶媒の内部自由度 (イオンによる電荷分布，配向自由度など) を適切な密度汎関数を通して考慮することも可能である [1-3] が，本研究では溶媒がニュートン流体の場合について多粒子系の流動を検討した．
SPM は多粒子問題の効率的な解法ではあるが，溶媒流動および微粒子の直接計算であるため，対象とする問題によっては，大規模系の計算が必要となる．そこで， 共有メモリ環境において openMP を用いて並列化を行った．
Fig.1 に速度向上の様子を示す．Fig.1 の場合では並列化率 95%を達成している．

Fig. 1: 共有メモリ環境における openMP を用いた並列化による速度向上.

3成果

工業的に重要なプロセスである多粒子分散系のレオロジーおよびコロイドの沈降過程の解明について検討を行った. 微粒子分散系は，コロイド懸濁液，スラリー，塗料などに見出すことができ，その粘度挙動を内部流動挙動から理解することが重要である. 熱揺らぎを考慮した SPM[4] を用いて，単純せん断流下のせん断応力を求め，微粒子分散系の粘度挙動についてを検討した．Fig.2 は高濃度分散系の粘度を計算したものである．
低レイノルズ数の計算データは実験相関式とよく符合しており， SPM のレオロジー予測能力を示している．
20Einstein relation
Kriger-Dougherty

18

Shear viscoity

161412108642000.1

0.20.30.4
Volume fraction

0.50.6

Fig. 2: ゼロせん断粘度 vs. 体積分率.

Fig. 4: SPM による沈降粒子のスナップショット．重力は垂直下向きに作用している．矢印は溶媒の速度場を示している．

4Fig. 3: 体積分率 0.51 におけるシアシックニングと構造因子 S(k). 横軸はせん断速度で，縦軸は粘度．
微粒子分散系では，溶媒にニュートン流体を用いていても，粒子の構造化により粘度はせん断速度に依存し，非ニュートン的な流動を示す．したがって粘度挙動を理解す
るには，せん断速度および体積分率による粘度変化を粒子構造とともに理解することが必要である． Fig.3 に高濃度(体積分率 ϕ = 0.51) の場合の粘度をせん断速度に対して示す．低せん断時のほぼ一定粘度から急激に粘度が増加する様子がわかる．高濃度・高せん断速度における粘度の増加 (シアシックニング) は実験においても報告がなされている．Fig.3 には，各せん断速度における構造因子 S(k)が示されている．これらの結果から，Fig.3 のシアシックニングは，粒子構造の秩序-無秩序転移に伴って生じていることが示唆される．微粒子分散系のレオロジーは未解明の部分が多く，直接数値シミュレーションによりミクロな粒子構造とマクロなレオロジーを同時に検討することにより明らかにされていくものと考えている．
Fig.4 は、沈降粒子のシミュレーションのスナップショットを示している．流動場の空間相関が粒子直径の数倍程度の長距離相関をもつことがわかる．本シミュレーションでは，微粒子の密度揺らぎがシステムサイズ程度の大きさを持つことが観察されており，計算の大規模化により意味のある結果を生み出すチャレンジングな問題である．
これらの微粒子分散系では微粒子の集団的なゆらぎとシステムサイズの競合が重要となる場合があるため大規模な直接計算が必要である．本研究により微粒子・溶媒・熱揺らぎおよび微粒子集団ゆらぎと多スケールにわたる階層構造と流動の様相が明らかにされると考える．

今後の課題

現状では、粒子直径の 10 倍程度のシステムサイズに対して意味のあるシミュレーションを行うために， 10GFlops級の計算機を用いて 0.5 日程度が必要である．このような小さな系のシミュレーションは周期境界条件の下で行うのが現実的であり，解決できる問題は外部境界条件によらない等方的な問題のみに限られる．粒子直径の 100 倍程度のシステムサイズに対しては，10PFlops 級の計算機で 5日程度である．これぐらい大きな系のシミュレーションであれば，現実的な外部境界条件を設定することが可能になり，粒子沈降や凝集を伴う微粒子分散系の問題に対して，初めて直接シミュレーションで答えを出すことができるようになる．ここまで系のサイズが大きくなくとも，例えば粒子サイズの 50 倍程度のシステムサイズであれば，1 ランあたり 1PFlops 級の計算機で 3 日程度の計算量でよい．
現段階では，OpenMP による並列化のみの実装であるが，今年度後期より京都大学学術情報メディアセンターと協力して，MPI ＋ OpenMP による大規模化に向けたプログラム開発に取り組む計画である．また，高速フーリエ変換 (FFT) を用いたスペクトル法を使用しているが，FFTの大規模化は困難である．MPI の使用を念頭に置いた場
合，この点を改善することが期待される．具体的な手法としては，波数空間のスペクトル法から，領域分割で並列化が比較的容易な実空間の緩和法にリプレースすることを考えている．
参考文献[1] Y. Nakayama, K. Kim, and R. Yamamoto: Simulating (electro)hydrodynamic eﬀects in colloidal dispersions: smoothed proﬁle method, Eur. Phys. J. E.
(2008) in print.
[2] K. Kim, Y. Nakayama and R. Yamamoto:Direct numerical simulations of electrophoresis,Phys.Rev.Lett., 96, 208302 (2006).
[3]R. Yamamoto: Simulating particle dispersions innematic liquid-crystal solvents , Phys. Rev. Lett. 87,075502(4) (2001)[4] T. Iwashita, Y. Nakayama and R. Yamamoto: A
numerical model for Brownian particles ﬂuctuating ,J.Phys.Soc.Jpn, 77, 074007 (2008).