RISTニュース

No．68（2022）

分子シミュレーションによるダイラタント現象へのアプローチ

A Molecular Dynamics Simulation of theShear Thickening Phenomena
高度情報科学技術研究機構
城野 亮太、手島 正吾
概要
液体にせん断を加えたときに粘度が増し固体のように振る舞うことをダイラタンシーといい、身近なところでは水溶き片栗粉などに現れる流体力学の一現象である。ダイラタント現象の応用先としてスペースデブリ（宇宙ごみ）など高速に飛来する物体に対する柔軟かつ耐衝撃に優れた宇宙服材料への展開など、その用途に大きな可能性を秘めているが、粘度変化の起源、特に分子論的なメカニズムは未だわかっていない．我々は非平衡分子動力学法によって流体粒子中の剛体球のシミュレーションを用いたダイラタント現象の発現に成功した。剛体球周りおよび流体粒子周りの 2 次元対分布関数から、ダイラタント現象が発現するときには微細構造が異方性を持つことを示し、流体粒子と剛体球の大きさの違いに起因する速度の違いによってこの異方性が現れることを見出した。この速度差によって流体粒子と剛体球の衝突が引き起こされ系の粘度が上昇することを明らかにした［ 1 ］。
1 ．はじめにであるナビエ・ストークス方程式を解くこと

ダイラタント現象の学術研究では剛体球をによって得られる。液体のような非圧縮性流ニュートン流体に浸したものとみなせるシリ体に限定し、またレイノルズ数と呼ばれるカ粒子（アモルファスSiO2）とポリエチレン系の密度ρ、代表長さL、代表速さU、粘度グリコール（polyethylene glycol, PEG）水ηの比 Re＝ρLU /ηの値が十分に小さい時、溶液の組み合わせが実験・計算共に良く調べナビエ・ストークス方程式はストークス方程られている［2－5］。特に剛体に働く相互作式と呼ばれる形に簡略化できる。

用のうち流体からの寄与をストークス方程式に基づいて計算し、剛体間相互作用はニュートンの運動方程式を解くことによって時間発展させるストークスダイナミクスという手法

ここで流体の流速 v、圧力 p、外力であは、ダイラタント流体の起源に迫る成果を上る。この方程式を、溶液中を沈んでゆく球形げきた［6－9］。この方法では、ニュートン粒子のような系について適用すると解析解 F

流体に浸された剛体に作用する力は、剛体と＝6πμavが得られ、Stokesの抵抗法則とし流体の間の相互作用を特徴付ける流体力学的て知られている。しかし初期速度を設定した相互作用、剛体間の相互作用を表す非流体力、ときのナビエ・ストークス方程式の一般解および系の温度によるランダム力である。流

は懸賞問題［10］となるほどに難問で、残体力学的相互作用は、流体を支配する方程式

念ながらストークス方程式にも一般解は与え－20－

RISTニュース

No．68（2022）られていない。そこで解析解が与えられてい・ )＝－・ のように表される。一てη(γ〈Pxy〉/γる二体問題［11、12］をベースとして、数般に分子シミュレーションにおける応力テン値的に多体問題にアプローチする方法が取ら

ソルの各要素は揺らぎが大きく収束しにくいれてきた。最近の研究では固体粒子間の摩擦ことが知られている。実際本研究で用いたがダイラタンシーを再現する上で重要である

3000粒子系で応力テンソルのアンサンブルことがわかっている。また、粒子周りの微細平均を求めると標準偏差が6.5×108 Paもの構造解析により粘度変化が起き始める臨界せ揺らぎとなり、例えば水の粘度は10－3 Pa s

ん断速度前後での特徴的な構造変化が実験であるようにmPa s（ミリパスカル秒）のオー［13、14］・計算［ 6 ］の両面で観測されてダーで議論される物理量を求めるには心もといる。しかし、流体力学的相互作用の起源とない。応力テンソルをせん断速度で除した粘その振る舞いは、分子レベルでの相互作用と度の標準偏差が十分小さな値になるためにはいう観点からは完全には明らかではない。さ少なくともせん断速度が108 s－1よりも速くらにこれらのマクロスケールシミュレーショなくてはならない。したがってこのような速ン手法では流体が連続体モデルとして扱われいせん断速度でダイラタンシーが起きる条件るため、流体の微細構造についてはこれまでを探す必要がある。

議論されてこなかった。

我々はペクレ数と呼ばれる移流（せん断速我々は、系の支配方程式をニュートンの運度に特徴づけられる流れ）と拡散（分子の熱動方程式F＝maに求め、流体を構成する粒揺らぎに特徴づけられる量）の比が、100程子についても剛体球と同じように取り扱うこ度になるとダイラタンシーが起きる傾向があとで、ダイラタント現象発現の分子論的メカるという報告［16、17］に着目した。ペクニズムを明らかにしようと考えた。

レ数Peは拡散係数Dと代表長さa、せん断速・ を用いてPe＝γ・ a2/D のように表される。
度γ2 ．計算方法と系の構築

ペクレ数Peが100になるためには、系の代表系は同種の3000粒子からなるとした。粒長さを正十二面体による剛体球の粒径 1 nm、径約 1 nmの剛体球は正十二面体の頂点20点

せん断速度を108 s－1とすると、拡散係数Dはに粒子を配置し、その粒子間に掛かる力が 2～10－11 m2 s－1でなくてはならない。用いた×10－3 kcal・mol－1 Å－2以下になるまで構造力場での拡散係数が～10－9 m2 s－1 であるこ最適化することで得た（図 1 ）。シミュレーとから質量を106倍することで分子運動を抑

ション中は結合長および結合角を固定し剛体として取り扱った。この剛体球を0、10、15、20、25個配置し残りを同種の粒子で流体粒子として埋めることで総粒子数が等しい系を構築した。例えば剛体球が20個の系では、400粒子（20粒子×20個）が剛体球を、残り2600個の粒子が流体粒子を構成している。粒子間に働く力はStillinger-Weber 型の力場［15］を用いた。
・ における粘度η・ ）は応力
せん断速度γ（γ

図1テンソルPxyのアンサンブル平均〈Pxy〉を使っ－21－正十二面体の頂点に粒子を配置することによって作成した剛体球

RISTニュース

No．68（2022）

制し、拡散係数が 5×10－12 m2 s－1の系を作

ニュートン流体領域ではPeが増加するに連成した。このような質量改変により得られたれより急峻になっていく。しかし臨界せん断剛体球の重さ（6×10－16 g）と流体粒子の示速度Pe＝50.0を超えてダイラタント流体領域す拡散係数（5×10－12 m2 s－1）は、ダイラに入るとこれらの溶媒和構造は壊れ、特に第タンシーの研究に広く用いられるポリエチレ二溶媒和圏については失われた。これはせん

ングリコール中のシリカ球の系における粒断によって流体粒子間の距離がある特定の値

径40nmのシリカ球の重さとポリエチレングに制限され構造化していくが、あるせん断速リコールの拡散係数と似通っていた。この度を境にそのような構造を維持できなくなる様にして得られた計算条件で温度圧力一定ことを示す。このような変化は図 3（b）に示の平衡化を行い密度が 1 kg mm－3 の系を得した剛体球重心周りの流体粒子の動径分布関・ の非平衡分子動力学計算はた。せん断速度γ数

LAMMPS［18］ に 実 装 さ れ て い るSLLOD第一溶媒和圏に相当する6.5Åのピークの振法［19］を用いて行った。系の温度は能勢

る舞いは流体粒子の動径分布関数-Hoover 熱浴により300 Kになるよう制御し同様の傾向を示したが、臨界せん断速度を超た。時間発展は10 fs刻みで行い、定常状態をえると新たに5.6Åにピークが現れた。この得るための準備計算の後、99 nsの計算を行っピークは剛体球を構成する正五角形面のホロた。

ウサイト（中空位）に流体粒子が配位した構（r）ではより顕著である。剛体球のRp（r）とpp造に対応し（図 3（b））、ペクレ数の増加とと3 ．結果もに成長した。したがってこの構造がダイラ3.1

タンシーの起源である可能性がある。一方、レオロジープロット

図 2 はペクレ数に対する粘度の値をプロッ

図 3（c）に示したように剛体球重心間の動径トしたものである。参照系とした剛体球が含分布関数

まれていない系を除いた全ての曲線はダイラ

Åのピークは、せん断をかけることによってタンシーの挙動を示した。剛体球の体積分率

ニュートン流体・ダイラタント流体問わず構あるいは濃度が上昇するに連れ、ダイラタン造が壊れ、臨界ペクレ数からダイラタント流（r）は平衡状態における10.3RR

シーが始まるペクレ数の値である臨界ペクレ数は減少し、これまでの一般的な実験の報告と一致した。このダイラタンシー発現機構を明らかにするため剛体球20個を流体粒子2600個に浸した系（図 2 の緑線）について詳細に解析した。以降この臨界ペクレ数50を境に低ペクレ数側をニュートン流体領域、高ペクレ数領域をダイラタント流体領域とする。
3.2

動径分布関数

図 3 は系の動径分布関数(r)である。流体粒子周りの流体粒子の動径分布関数

pp（r）（図 3（a））ではせん断なしの平衡状態（Pe＝0）で見られた第一溶媒和圏および第二溶媒和圏を示す3 Åおよび4.5 Åのピークは、－22－

図2剛体球が0、10、15、20、25個入った系のレオロジープロット

RISTニュース

No．68（2022）

体領域に移行すると剛体球同士が接近してい粒子が剛体球へ追いつけず、分布として疎なた。したがって、ダイラタンシーがおきてい領域となって現れている。このようにして点るときには剛体球と流体粒子の間あるいは剛対称異方的な対分布関数が形成された。本研

体球間の接近、すなわち衝突があり、これがダイラタンシーの起源である可能性がある。
3.3

微細構造

先の動径分布関数を、せん断流れに沿った方向（x軸）とそれに垂直なせん断速度の勾配方向（y軸）に分けて対分布関数 h（x, y）として表示したものが図 4 である。 9 つの分布図のうち、左列はニュートン流体領域、中列は臨界せん断速度、右列はダイラタント流体領域のときの分布であり、上段は流体粒子周りの流体粒子の対分布 hpp（x, y）、中段
は剛体球周りの流体粒子の対分布 hRp（x, y）、下段は剛体球周りの剛体球の対分布 hRR（x,y）である。流体粒子間の対分布関数 hpp（x,y）は内側に第一溶媒和圏に相当するほぼ真円の分布があり、ニュートン流体領域からダイラタント流体領域にかけて大きな変化は無い。一方で剛体球周りの流体粒子の対分布
hRp（x, y）については臨界せん断速度までとダイラタント流体領域とで顕著な違いが見られた。ニュートン流体領域では等方的であった hRp（x, y）は、ダイラタント流体領域で異方的であった。
分布図の原点を剛体球重心に取った場合、せん断流れに沿った方向（x軸）とそれに垂直なせん断速度の勾配方向（y軸）について
y > 0, x < 0および y < 0, x > 0の領域では、流体粒子は剛体球重心よりも速いため流体粒子が図3

剛体球へ衝突する。分子運動に起因する速度に比べてせん断流れが速い場合、粒子にとってせん断流れによって生じるせん断流れ層を行き来することは難しいため、小さな流体粒子は大きな剛体球に堰き止められ大きな力がかかり、粘度上昇の原因となっている。一方y > 0, x > 0および y < 0, x < 0の領域では、剛体球の速度に比べ流体粒子の速度が遅く、流体－23－20個の剛体球を2600個の流体粒子中に浸した系の動径分布関数。（a）流体粒子-流体粒子間、（b）剛体球の重心と流体粒子間、（c）剛体球の重心間。いずれのグラフもペクレ数Peが0（せん断をかけない平衡状態）、42.9（ニュートン流体として振る舞っている非平衡定常状態）、50.0（臨界せん断速度）、51.4（ダイラタンシーが発現している非平衡定常状態）での動径分布関数を示した。

RISTニュース

No．68（2022）

究では粒子間に働く相互作用は分子間力のみいない。このようなマイクロスケールで初めを仮定し、マクロスケールで導入されているて定義される量について分子レベルから知見潤滑力や摩擦力といった項を直接取り入れてが与えられると考える。

図4

流体粒子周りの流体粒子の対分布関数（a）－（c）、剛体球周りの流体粒子の対分布関数（d）－（f）、剛体球周りの剛体球の対分布関数（g）
－（i）。
（a）、
（d）、
（g）はニュートン流体領域（Pe＝42.9）、（b）、（e）、（h）は臨界せん断速度（Pe＝50.0）、（c）、（f）、（i）はダイラタント流体領域での分布を示す。

－24－

RISTニュース

No．68（2022）

剛体球周りの剛体球の対分布hRR（x, y）はJ-I., Microstructure of the Fluid

ニュートン流体領域ではことなるせん断流れ

Particles around the Rigid Body

層にある剛体球同士の衝突を示していたもの

at the Shear-thickening Stateの系の粘度としては大きく変化せず、ペクレ

toward Understanding of the Fluid

数が増加するに連れ同じせん断流れ層にある

Mechanics, Sci. Rep.,11, 24204 (2021)

剛体球同士でしか接近できなくなることがわ［ 2 ］ Brown, E., Jaeger, H. M., Shear

かった。したがって衝突には剛体球と流体粒

thickening

inconcentrated子のような異なる大きさの粒子が存在するこ

suspensions:とが重要であると考えられる。

mechanisms and relations tophenomenology,

jamming, Rep. Prog. Phys.77, 046602(2014)

4 ．おわりに
本稿ではダイラタンシーという現象につい［ 3 ］ Hasanzadeh, M., Mottaghitalab,て非平衡分子動力学によるアプローチを行

V., The Role of Shear-Thickening

い、粒子間の衝突によって粘度変化が引きFluids (STFs) in Ballistic and Stab-

起こされていることを示す結果を紹介した。

Resistance Improvement of Flexible

ニュートンの運動方程式から流体力学の支配

Armor, J. Mater. Eng. Perform., 23,

方程式であるナビエ・ストークス方程式にア1182 (2014)

プローチする方法は、マクロスケールで初め［ 4 ］ Kamibayashi, M., Ogura, H.,て定義される摩擦等の語句をミクロスケール

Otsubo, Y., Shear-thickening ﬂow ofで説明できる可能性があり現象の理解・制御

nanoparticle suspensions ﬂocculatedに大きく役立つと考えている。今後は各階層

by polymer bridging, J.Colloid.Inter.の支配方程式を統一的に理解するための枠組

Sci., 321, 294 (2008)

みを構築することでマクロ現象とミクロ現象［ 5 ］ Lee, J., Jiang, Z., Wang, J., Sandy, A.
R., Narayanan, S., Unraveling theを理解していきたい。

Role of Order-to-Disorder Transitionin Shear Thickening Suspensions,

謝辞

Phys.Rev.Lett., 120, 028002 (2018)

本研究は、防衛装備庁安全保障技術研究推進制度JPJ004596の支援を受けて実施しま［ 6 ］ Brady, J. F., Bossis, G., Stokesian

した。また、本研究成果は筑波大学計算科

Dynamics, Annu. Rev. Fluid Mech.

学研究センターの学際共同利用プログラム

20, 111 (1988).

（cygnus）を利用して得られたものです。プ［ 7 ］ Chen, K., Wang, Y., Xuan, S., Gong,

ロジェクトリーダーの筑波大学藤田淳一教

X., A hybrid molecular dynamics

授、メンバーの筑波大学赤田圭史助教、住友

study

電気工業株式会社新領域技術研究所小野木伯

rheological behaviors of shear

薫様、大久保総一郎様、藤森利彦様には細部

thickening fluid, J.Colloid Interf.

にわたる議論をしていただきました。この場

Sci., 497, 378 (2017)onthenon-Newtonian［ 8 ］ Nakayama, Y., Yamamoto, R.,

を借りて感謝申し上げます。

Simulation Method to Resolve
Hydrodynamic Interactions in

参考文献［ 1 ］ Jono, R.,

Tejima, T.,

Fujita,Colloidal Dispersions, Phys.Rev. E－25－

RISTニュース

No．68（2022）

71, 036707 (2005)［15］ C h a n , H . ,［ 9 ］ Seto, R., Mari, R., Morris, J. F.,

Cherukara, M. J.,

Narayanan,B.,

Loeffler,T.

Denn, M. M., Discontinuous ShearD., Benmore, C., Gray, S. K.,

Thickening of Frictional Hard-Sankaranarayanan, S. K. R. S.,

Sphere Suspensions, Phys.Rev Lett.,

Machine Learning Coarse Grained

111, 218301 (2013)

Models for Water, Nat.Commun., 10,［10］ h t t p : / / w w w . c l a y m a t h . o r g /millennium-problems/navier-stokes-

379 (2019)［16］ N a k a m u r a ,equation

H.

Ishii,M.,

Makino, S., Dominant Model of［11］ Brenner, H., O’Neill, M. E., On the

Shear Thickening Behavior of

Stokes Resistance of Multiparticle

Concentrated

Systems in a Linear Shear Field,

Monodispersed Colloidal Particles,

Chem.Eng.Sci., 27, 1421 (1972)

J.Soc.PowderTechnol.Jpn., 56, 438［12］ Mazur, P., Saarloos, W.V., Many-Suspensionsof(2019)

Sphere Hydrodynamic Interactions［17］ Bender, J., Wagner, N. J., Reversible

and Mobilities in a Suspension,

Shear Thickening in Monodisperse

Physica, 115A, 21 (1982)

and［13］ Parsi,F.,

Gadala-Maria,F.,

Bidisperse

Dispersions, J.Rheology, 40, 899

Fore-and-Aft Asymmetry in a(1996)

Concentrated Suspension of Solid［18］ P l i m p t o n ,Spheres, J.Rheol., 31, 825 (1987)［14］ C h e n g , X . ,McCoy,

Colloidal

AlgorithmsJ. H.,

S.,
for

Fast

Parallel

Short-Range

Molecular Dynamics, J Comp Phys,Israelachvili, J. N., Cohen, I.,

117, 1 (1995)

Imaging the Microscopic Structure［19］ E v a n s , D . J . , M o r r i s s , G . P . ,of Shear Thinning and Thickening

Nonlinear-Response Theory for

Colloidal Suspensions, Science, 333,

Steady Planar Couette Flow, Phys

1276 (2011)

Rev A, 30, 1528 (1984)－26－