（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 ダイラタント系のシアシックニング機構の解明に向けた理論シミュレーション
英文： A theoretical study on the mechanisms of shear thickening
利用課題責任者
山田達矢
所属
高度情報科学技術研究機構・計算科学技術部
邦文抄録（300 字程度）： 我々の研究グループでは，分子シミュレーションの方法を使用して，高分子系や粒子分散溶液系のダイラタント性（あるいはシアシック二ング）に関する理論的検討を行っている．本研究では，ポリエチレングリコール(PEG)/シリカ粒子分散水溶液系の散逸粒子動力学法(DPD)シミュレーションを行った．この系は典型的なダイラタント系で，その振る舞いは古くからよく知られているが，対応する分子レベルの現象はわかっていない．適当なモデルを用いて一連のシミュレーションを行い，構造や粘弾性に関する解析を行ってみた．
英文抄録（100 words 程度）: Shear thickening/thinning behaviors of dilatant systems of polymer solutions andnano-particle-suspended solutions have been studied in our group with molecular simulations. We focused in thisstudy on a dilatant system of polyethylene glycol solution with silica particle suspension. We analyzed molecularlevel structures and viscoelastic properties of this system through a series of dissipative particle dynamics (DPD)
simulations.
Keywords（５つ程度）: dissipative particle dynamics (DPD), ダイラタンシー，ポリエチレングリコール/シリカ粒子分
散溶液，レオロジー
背景と目的ステップ間隔を大きく取れるため，ポリマーや生体分子遅い剪断刺激に対しては液体として振る舞う物質が，などの大きな分子系の構造を調べるのにしばしば用い速い剪断刺激に対して固体のような抵抗力を示す場合られる方法である．本研究でもこの理由で DPD 法を使があることが知られている（ダイラタント現象）．この種用したが，ポリマーダイラタント系への応用の前例は無の物質を衝撃緩衝材や制動装置などへ応用する研究く，ここでの調査は方法評価を兼ねる．

が行われているが，剪断速度に応じて粘度の増加する（シアシックニング）物理的メカニズムが十分に理解さ結果および考察れていないため，いまのところ材料設計の適当な指針

図１のように PEG 水溶液/シリカ粒子分散系を粗視化が確立していない．我々は，分子動力学および粗視化（DPD 粒子化）して，表１のような大小の構成のシミュレ動力学シミュレーションを研究方法として用い，シアシッーションモデルを用意した．これらモデルの構成にはダクニングの分子レベルの現象解明に取り組んでいる．

イラタント現象の観測されている実験系の物質構成[1]

本研究では，ポリエチレングリコール(PEG)/シリカ粒子を参考にしたが，実験系とは表１脚注のような相違があ分散水溶液系のダイラタント性に関する検討を行った．

るので注意されたい．粒子間相互作用についてはあまり情報がないため，暫定的に表 2 の二種のソフトポテン

概要
この系のダイラタント性は PEG とシリカのインタープレイによるものと考えられているが．構造に関する具体的シャルパラメータセットを用意した: 異種粒子が特に強い親和性や排他性を持たないもの（SP1）と，P 粒子と S粒子が親和するもの（SP2）である．

データは報告されていない．本研究では，適当なモデルを用意して散逸粒子動力学法（DPD）シミュレーション

A. 計算可能なシミュレーション時間に関する検討による解析を行ってみた．DPD 法は粗視化動力学計算

はじめに DPD 法によるダイラタント現象の直接シミュ法の一種で，通常の分子動力学計算に比較してタイムレーションの可能性について検討した：ダイラタント現象

（様式 20）成果報告書の時間スケールは 10-1-10-3sec の程度で（101-103/sec の解析という側面でも（A の節で述べた事柄に加えて）検剪断速度で観測される），これは分子シミュレーション討できる剪断速度領域が制約を受けることを象徴的にの時間スケール（タイムステップが fsec や psec オーダ示している．この問題は MD 法の計算については指摘ー）に比べるとずっと大きく，DPD 法でも実時間的に十されていたが，DPD 法でも同様であることが分かった．

分なシミュレーションができない懸念があったからであ調べた剪断速度領域では shear thinning という結果る．になった．また viscosity の剪断速度ゼロの極限を平衡計算を実行したところ，表 2 の大小のモデル系の計系の計算から Green-Kubo 法により求めたところ，5-6h0

算速度は，スパコン 4 ノード使用でそれぞれ~2M との 値 が 得 ら れ た の で ， 0-10-6t0-1 の ど こ か で shear~20M ステップ/日であった．したがって，それぞれのモthickening の剪断速度領域があることが示唆される．

デルを使って実行できる計算総量は，現実的にはしかしながら計算で得られた viscosity の増減は実験系~100M と~1G ステップの程度で，シミュレーション実時に比べると小さく（実験系では値が数桁変化する），ダ間は~1Mt0 （数十 µsec）と~10Mt0 （数百 µsec）の程度まイラタント現象に対応する変化かどうかは疑問である．

でと試算される（１ステップは 0.01t0, 1t0 = 25 psec, DPD

ここで得られている viscosity の数値の由来を検討す

単位については図 1 脚注参照）．非平衡動力学計算でるため，2 成分や 1 成分系の計算を行ってみた．（表 3．

は定常状態到達までにある程度大きな剪断変形を与簡易な検討のために小モデルを使用した．また統計解える必要があるが，これを仮にg=1 とすると，定常状態析に要する計算ステップ数を少なく抑えるため，あまりを 調 べ ら れ る の は ， 小 さ く て せ い ぜ い ~104/sec や小さくない剪断速度ということで 10-4/t0 の系を調べた）．

~103/sec の剪断速度までということになる．斯くして，ダ計算モデルからシリカを取り去ると viscosity が若干上イラタント現象の剪断速度領域（低粘度-高粘度の遷移昇する，また PEG を取り去ると viscosity が極端に減少すの領域）を調べるのは困難であるということが分かっる（溶媒レベルの値まで）という結果が得られていて，こた．

れは系の viscosity の殆どが PEG 由来であることを意味している．また，viscosity は PEG を寸断すると極端に減B. 相互作用モデル SP1 を用いた DPD 計算
次に，モデル系の不具合や DPD 法の計算手法に絡少するため，PEG の紐状に長いという特徴と関係したものであることを示している．

んだ問題などについて知見を得るため，10-2/t0 から各剪断速度での構造を図 3 に，PEG の構造パラメー

10-6/t0 （4•108/sec から 4•104/sec）の範囲の剪断速度のタを表 4 に示した．これらのデータを全体的に見ると，シミュレーションを行い，viscosity や構造に関する一連剪断速度の遅い非平衡系(10-6/t0)と平衡系の構造にはの解析を行ってみた．ここでは表 2 の SP1 の相互作用違いが殆ど無く，ある程度剪断速度が大きくモデルを使用した．

（10-4-10-3/t0 の辺り）なると違いが明確になるようである．

大モデル系の viscosity 解析の結果を図 2 に示した

見た目の構造については（図 3），剪断速度が 10-3/t0 以（小モデル系の結果は大モデル系と同様なので割愛）．

上で PEG 鎖が剪断速度方向に並ぶという特徴（平衡系大きい揺らぎのために viscosity の数値をあまり正確にでは見られない）を持つことがわかる． PEG の広がり具計算できない問題のあることが明確になった．この問題

合を表す両端距離や慣性半径（表 4）は剪断速度に応はそもそも応力計算値の揺らぎが大きいことに由来しじて比較的にはっきり変化する：慣性半径のデータを見ていて，剪断速度が小さくなるに従って viscosity（=応力ると，剪断速度が大きくなるのに伴って PEG が剪断速度÷剪断速度）の揺らぎが深刻になっているためである方向（x 軸方向）に引き伸ばされ，これに垂直な方向（y（応力計算値の揺らぎが剪断速度によらずにほぼ一定および z 軸方向）には広がりを失うことがわかる．PEG のになっていることにも注意されたい）．剪断速度の小さ結合長や結合角といった局所的な構造は，基本的に剪

い 10-6/t0 については， 100M 超のステップ数を費やし断速度に応じて値を変えない（表 4）ので，上記の変形ても大きなエラーバーが残っており， viscosity の統計は PEG の 2 面角捩れを通じた折り畳みの具合によるも（様式 20）成果報告書のである．表 4 のデータを総合的に見ると，系の構造変強いられていることがわかる．この間の viscosity の数値化はどこかの剪断速度で急激に起こるというタイプのもはシリカの無い場合の値（~9.9 h0，表 3 参照，シリカののではないようである．

無い系では SP1 と SP2 の区別はないことに注意）より有以上，パラメータセット SP1 を用いた計算の結果を詳意に大きい．斯くして，速度の異なる剪断流を跨ぐよう細に調べたが，次のように実験と相容れない部分があな PEG•シリカネットワーク（図 4 でいえば縦方向の隣接る． i) viscosity がほぼ PEG 由来でシリカとの interplay

ユニットセルに跨るもの）を持てば高 viscosity 状態を与による寄与が無い，ii) 平衡系と非平衡系で構造があえることが分かった．

まり違わず（10-6/t0 ~ 4•104/sec は実験的にはかなり大一方で ~20M ステップを過ぎると，viscosity が極端にきな速度である），iii) 違いがあっても剪断速度に応じ低下するという結果になった．これは，上記の引っ張りて連続的に変化するものである（ダイラタント転移の前のために，シリカに張り付いた PEG が徐々にはがれ，終後では急激な構造変化が起こるともの考えられている）．には上述の剪断流を跨ぐ架橋が失われてしまったため

以上の事柄はここで使用したモデルに問題があり修正である．一旦こうなると，剪断流を跨ぐ架橋が再生するの必要があることを意味している．

ことはない（剪断流方向(横軸方向)の架橋は保持されるが）．この構造では PEG がシリカに緊密に張り付いて，C. 相互作用モデル SP2 を用いた DPD 計算

またシリカ球が剪断流方向に整列するため，PEG の縦上記の i) の点に関係したモデルの変更として，P 粒軸方向（剪断速度勾配方向）の空間的な広がりが 1 個子-S 粒子間の親和性を強くしたモデル SP2（表 2）を用いのシリカ球程度しかない（図 4，表 5 の rgy_y）．このためた計算を行ってみた．ここでは前節のような大規模な検に viscosity が溶媒並みに小さい値になるというわけで討はしないで，剪断変形に対する系の応答を簡易に検ある．斯くして，1-20M ステップの高 viscosity 構造は過討するということにして，1x1x1 モデルを使用し 10-4/t0渡的な短寿命のもので， むしろ 24M step 以降の低の剪断速度の非平衡系の計算だけを行った．

viscosity 構造が定常状態である．

計算結果を図 4 と表 5 に示した（これに対応した SP1

こうして考えると，図 4 の 1-20M ステップの構造が高による計算の結果を図 5 と表 6 に参考のために掲載しviscosity であるからと言って PEG 水溶液シリカ分散系のた）．相互作用モデル SP2 を用いた計算では，SP1 の場ダイラタント状態に単純に対応するとは言えない．ここ

合と対照的に， PEG のほとんどのセグメントがシリカ表で行ったシミュレーションが実験に忠実なものではない面に張り付き（図 4 と 5 の initial structure 参照）, これに（剪断速度（10-4/t0 ~ 4•106/sec）が実験系でダイラタンシ付随してシリカ球が PEG を介して会合するのが特徴でーの観測される剪断速度（101-103 /sec）よりずっと大きある． また，PEG のセグメントの一部は，このシリカ球クいなどの違いがある）ことには注意する必要があるが，ラスターの間を架橋する．このために，ポリマー・シリカ

上述のような架橋には剪断速度の大小に関係なく寿命のネットワークが実効的に空間的に無限に広がった格があるはずなので，高 viscosity が維持されるメカニズム好になっている．この構造は，PEG 水溶液シリカ分散系が何か無くてはならない．この点で現状モデルには欠の高 viscosity 状態として ref.1 において想像されている陥があることは明らかであり，今後の検討のためにはものに近い．

より適切なものを構築する必要がある．

この系に剪断変形を与えると，上記の架橋部が引っ張られる（架橋の上部と下部が異なる速度で横方向に動くため）ことにより, viscosity の大きい状態（~20M ステップまで, h ~10-20 h0）が現れる．このタイムステップ

まとめ、今後の課題
PEG/シリカ粒子分散水溶液系の DPD シミュレーションを行った．検討の結果，以下のような知見を得た．の間，平衡系に比べて PEG の結合角が大きい，結合偏i) DPD 法で定常化できる剪断速度系は小さくてもせ角のエネルギーが大きい（表 5）などの特徴が認められ，

いぜい~104/sec や~103/sec の剪断速度までで，ダイラタ剪断変形が与えられたことにより PEG が伸長の変形をント性の現れる剪断速度領域のシミュレーションは困難

（様式 20）成果報告書
である．ii) DPD 法による計算でも，MD 法の計算の場合と同様で，大きい揺らぎのために viscosity の数値を正確に計算できない問題が深刻である． iii) SP1 を SP2 の相互作用モデルを用い DPD 計算を行ってみたが，それぞれ，シリカの有無が viscosity にあまり寄与しない，定常状態が溶媒並みの viscosity しか示さない，など実験と相容れない結果を与えるものであり，今後の検討のためにはより適切なものを構築する必要がある．
参考文献:[1] M. Kamibayashi, H. Ogra, and Y. Otsubo, Journal of Colloidand Interface Science, 321, 294–301, (2008). 表 2 のモデル構
成は，シリカ 8 vol% (9 nm 径) / PEG 0.5w%(分子量 50 万)の実験系を参考にした．
[2] R. D. Groot and K. L. Rabone, Biophysical Journal, 81,
725-736 (2001). W 粒子の拡散係数を独自に計算し，この論文の方法に倣いt0 を算出した．
[3] Z. Shen, H. Ye, M. Kroger, and Y. Li, Phys. Chem. Chem. Phys.,
19,13294 (2017). 実際には「シリカ粒子の ダイラタント特性に関する ソフトウェア OCTA を使った 粗視化シミュレーション」 報告書の様に値を調整した．

（様式 20）成果報告書

図 1. PEG 水溶液/シリカ粒子分散系の DPD 粒子(P, S,W )への粗視化．
この粗視化により，エネルギー, 距離，質量，時間, viscosity の DPD
単位はそれぞれ E0 = kBT, L0 = 0.445 nm, m0 =18 amu., t0 = 24.7 psec
(1step = 10-2t0), h0 =1.1 mPa∙s になる[2].

表 1. シミュレーションモデルの構成
1x1x1
2x2x2(-P-)1140a
432(W)
87,402699,216(S1139)b
756a) P 粒子 1140 個を直鎖状に結合して PEG 分子をモデル化した(先行研究[3]の結合伸縮と偏角のポテンシャルを適用)．b) S 粒子の面心立方格子から球状(直径 ~ 9 L0)に領域を切りだしシリカクラスタを作成した(S-S 結合に調和ポテンシャル r0 = 0.78 L0, k=16.5 E0 を適用). これ
らのモデルの構成にはダイラタント現象の観測されている実験系の物質構成[1]を参考にしたが，次のような実験系との相違がある：i)PEG 分子１本の鎖長が実験系の 1/10； ii) シリカ球の径が実験系の1/2; iii) PEG の濃度が 10 倍大きい. i)と ii)は計算系のサイズを限定的にするためのもので，iii)は，PEG 濃度が実験系のままだと，性質の違いが顕著に観測できない（純溶媒系とあまり違わない）ので，その対処である. DPD 単位については図 1 脚注参照．

表 2.ソフトポテンシャル反発パラメ-タ(所謂 aij, E0 = kBT 単位)
SP1SP2a(P-P) a(P-W) a(P-S)
25a 26.3a25caa
2526.320ca(W-W)
25a25aa(W-S)
25c30ca(S-S)
35b35ba) 先行研究モデルの値[3]を適用. b) シリカクラスタが会合しないよう大きめの値に設定． c) SP1 では異種 DPD 粒子が特に強い親和性や排他性を特に持たないよう設定；一方 SP1 では P 粒子と S 粒子が親和的になるよう設定．

図 2. 2x2x2 モデル系の各剪断速度での viscosity 解析の結果： 記載の数値は平均値,標準誤差(角括弧)と計算ステップ数（丸括弧)．h(+0)
は Green-Kubo 法に基づき平衡系の計算データから解析．DPD 単位については図 1 脚注参照．

表 3. 構成の異なるモデル系の viscosity (h0 単位)*
P+W+S P+W
W+S
W
P’+W+S
7.70
9.94
0.88
0.55
0.98
[0.46] [0.63] [0.15] [0.17] [0.18]*1x1x1 モデルを使用し剪断速度 10-4/t0 で計算. [ ]は標準誤差.
P’: ポリマ-結合を寸断し２粒子分子にした．h0 =1.1 mPa∙s.

図 3. 2x2x2 モデルの各剪断速度での計算のスナップショット: P 粒子(シアン)，S 粒子(赤), W 粒子は見やすさのためプロット省略.剪断速度の方向は横軸(x 軸)方向で，速度勾配の方向は縦軸(y 軸)方向である．
DPD 単位については図 1 脚注参照．

表 4. 各剪断速度でのモデル系の構造パラメータ*
剪断速度 0/t0 10-6/t0 10-5/t0 10-4/t0 10-3/t0 10-2/t0
0.380 0.380
0.380 0.380 0.380 0.380l(L0)[0.016] [0.016] [0.016] [0.016] [0.016] [0.016]
140.15 140.14 140.15 140.16 140.27 140.61
q(deg.)[0.009] [0.070] [0.066] [0.066] [0.010] [0.072]
34.6
34.837.259.5
123.8 154.2
de-e(L0)[0.03] [1.37] [2.07] [2.33] [5.09] [1.08]
15.3
15.3
16.0
23.3
45.7
59.7
rgy(L0)[0.05] [0.36] [0.48] [0.62] [1.49] [0.14]
7.78
8.03
9.10
19.99 44.92 60.52
rgy_x(L0)[0.17] [0.63] [0.20] [0.82] [1.62] [2.17]
8.81
8.74
8.38
7.43
5.07
3.14
rgy_y(L0)[0.28] [0.59] [0.28] [0.37] [0.10] [0.08]
8.61
8.42
8.94
7.78
6.72
4.18
rgy_z(L0)[0.12] [0.42] [0.68] [0.22] [0.28] [0.20]*各エントリ: PEG 結合長と結合角（l,q）, PEG 両端距離（de-e）, PEG 慣性
半径とその x-z 成分（rgy, rgy_x など，剪断速度は x 軸方向,速度勾配は y軸方向）. [ ]内数値は標準誤差．数値は DPD 単位で表現(DPD 単位については図 1 脚注参照)．

図 4. SP2 相互作用モデルを用いた計算における各ステップでの構造．
P 粒子(シアン，紫，緑. ピンクで 4 個の PEG 分子を区別)，S 粒子(赤),W 粒子は見やすさのためプロット省略. 剪断速度（10-4/t0 で計算）の方向は横(x 軸)方向で，速度勾配の方向は縦(y 軸)方向．初期構造は平衡系のスナップショットから取った．数値は viscosity の 4M ステップ
区間平均と標準誤差（DPD 単位については図 1 脚注参照）．

（様式 20）成果報告書
表 5. SP2 相互作用モデルによる計算結果（10-4/t0 で計算）

表 6. 前段の続き

steph(h0)
ref.*
10.1
0-4M
17.5
4M-8M 21.4
8M-12M 11.3
12M-16M 9.11
16M-20M 13.8
20M-24M 1.31
24M-28M -0.29

stepde-e(L0)
ref.*
42.1
0-4M
61.5
4M-8M
83.9
8M-12M
89.4
12M-16M 66.6
16M-20M 74.5
20M-24M 80.8
24M-28M 58.5l(L0)0.38000.38000.38000.38000.38000.38000.38000.3800q(deg.) Ebond(E0) Eangle(E0) Enb(E0)140.18 15862.5
140.93 15856.1
141.16 15855.8
141.10 15853.8
141.12 15854.2
141.15 15853.5141.04 15855.0
141.03 15853.1

2322.29 468867
2325.71 468885
2332.54 468911
2331.95 468886
2331.94 468852
2332.72 468852
2329.27 468830
2328.92 468786

rgy_x(L0)
10.7
17.0
27.4
28.7
22.2
25.1
25.5
19.8

rgy_y(L0)
9.1
12.5
9.9
7.9
7.1
8.3
8.3
6.6

rgy_z(L0)
8.6
6.9
6.6
8.6
8.3
7.9
8.5
6.7

*ref.: 平衡系(剪断速度 0)の結果．各エントリ: viscosity(h), PEG 結合

表 5. 前段の続き
stepde-e(L0)ref.*28.80-4M
52.3
4M-8M
97.1
8M-12M
90.612M-16M 88.8
16M-20M 88.3
20M-24M 68.7
24M-28M 51.9

rgy(L0)
17.4
23.2
30.6
31.5
25.1
28.1
28.9
22.5

rgy(L0)
13.5
21.4
34.0
31.5
31.0
30.5
25.9
22.2

rgy_x(L0)
4.2
16.0
30.5
28.5
27.4
27.3
22.7
18.9

rgy_y(L0)
9.6
8.7
7.4
5.8
5.4
5.4
3.9
4.1

rgy_z(L0)
7.0
7.6
6.8
6.3
6.9
6.3
6.7
6.8

*ref.: 平衡系(剪断速度 0)の結果．各エントリ: viscosity(h), PEG 結合長と結合角(l,q), ポテンシャルエネルギー成分(Exxx,結合伸縮・偏角と非結合相互作用), PEG 両端距離(de-e), PEG 慣性半径とその成分(rgy,剪断速度は x 軸方向,速度勾配は y 軸方向). 数値は DPD 単位（図 1脚注参照）で表現.

図 5. SP1 相互作用モデルを用いた計算における各ステップでの構造．
P 粒子(シアン，紫，緑. ピンクで 4 個の PEG 分子を区別)，S 粒子(赤),W 粒子は見やすさのためプロット省略. 剪断速度（10-4/t0 で計算）の方向は横(x 軸)方向で，速度勾配の方向は縦(y 軸)方向．初期構造は平衡系のスナップショットから取った．数値は viscosity の 4M ステップ
区間平均と標準誤差（DPD 単位については図 1 脚注参照）．

表 6. SP1 相互作用モデルによる計算結果（10-4/t0 で計算）
steph(h0) l(L0) q(deg.)ref.*~2 0.3800 140.120-4M
10.8 0.3800 140.18
4M-8M 17.6 0.3800 140.218M-12M 10.2 0.3800 140.16
12M-16M 6.7 0.3800 140.11
16M-20M 9.5 0.3800 140.15
20M-24M 8.5 0.3800 140.18
24M-28M 7.4 0.3800 140.17

Ebond(E0) Eangle(E0) Enb(E0)16080.1 2318.38 469798
16079.1 2318.78 469797
16080.0 2322.17 469785
16079.4 2321.52 469793
16082.0 2319.05 469802
16080.2 2320.31 469798
16079.9 2320.86 469794
16080.6 2319.96 469796

長と結合角(l,q), ポテンシャルエネルギー成分(Exxx,結合伸縮・偏角と非結合相互作用), PEG 両端距離(de-e), PEG 慣性半径とその成分(rgy,剪断速度は x 軸方向,速度勾配は y 軸方向). 数値は DPD 単位（図 1
脚注参照）で表現.