京都大学(京大)大学院工学研究科の山本量一教授と村島隆浩特定助教、谷口貴志准教授らの研究チームは、計算流体力学法と高分子動力学法の組み合わせによる「高分子流体のためのマルチスケールシミュレーション法」を用いて、円筒の周りを流れる高分子溶融体の複雑な流動挙動を再現することに成功した。同成果は近日中に欧州の物理学雑誌「Europhysics Letters」のオンライン速報版で公開される予定。

プラスチックに代表される高分子系材料は分子鎖が長く(分子量が大きく)、それらが複雑に絡み合っているため製品に十分な強度を持たせることが可能になる。プラスチック製品は、高温で溶融した高分子材料を金型に流し込む「射出成形」という方法で作られる。金型のような複雑な流路の中で高分子が流体として、どのような挙動を示すかを理解することは、材料や製造プロセスの設計を行う際に重要だが、これまでは工学的に重要な高分子は分子量が大きく、絡み合いがあるために高分子鎖の計算流体力学シミュレーションの計算負荷が大きいこと、ならびにマクロスケールとミクロスケールを組み合わせたマルチスケールシミュレーションは、高分子鎖の状態が流れに沿って同じである場合に限られていたこと、の2点から現実の高分子材料の問題に広く適用できていなかった。

同研究チームは今回、個々の流体粒子の計算負荷を軽減するモデルで絡み合い高分子鎖をシミュレーションし、マクロな計算流体力学シミュレーションに埋め込むことで、流体とともに流れるシミュレーション法を開発した。

具体的には、ソフトマターのシミュレーションが抱える難題をマクロスケールとミクロスケールを組み合わせた「マルチスケールモデリング」によって解決した。

図1はシミュレーションを行った系(30cm×30cm)の略図だが、中央に配置された円筒状の物体の直径(約6cm)は、高分子を構成する1つの原子のサイズの1024倍ある。高分子溶融体の巨視的な流れは、連続体モデルの運動量保存式によって計算され、ここでの各計算点での流速はその点での局所的な応力によって変化する。通常の粘性流体(ニュートン流体)の場合、応力と流速の空間勾配の間に比例関係があるので、その比例定数、すなわち粘性係数さえ実験で求めておけば、各点で求めた流速から各点での応力を求め、その瞬間の流速変化を予測することができる。

問題の概略図。流体は流体粒子として扱い、流体粒子の運動方程式を解く。また各流体粒子の応力は、流体粒子に関連づけられた粗視化高分子モデルにおいて計算される。図の中央の黒円は円筒を真上から見ているところを表す。黄色の矢印は流体粒子の流れの方向を示している

しかし、今回の研究で扱うような高分子溶融体にはそのような単純な比例関係は存在しないため、各計算点での応力はその点での巨視的な流速や変位だけでは表すことができず、原子・分子レベルでの高分子の配置、形状、絡み合いの状態に依存しているため、研究グループでは、マルチスケールモデリング法を用いて、まず巨視的な流動を表現するために流体を多数の流体粒子で表した。この流体粒子の運動方程式を解くためには応力場を知る必要があり、その応力場を求めるために、前のステップで求めた各流体粒子の速度からその位置における局所的な変形速度(速度の勾配)を求めて、それを外場として、各流体粒子に内在する膨大な数の高分子が絡み合った状態の時間変化を粗視化高分子動力学法により計算した。

この計算により求めたある流体粒子中の多数の高分子鎖の微視的情報を用いて、その流体粒子の位置における応力を統計的手法により求め、すべての流体粒子の位置で求めた応力により応力場を構成し、その応力場を用いて各流体粒子の運動方程式を解き、微小時間後のすべての流体粒子の位置と速度を決定する。この計算を繰り返し行うことで高分子流体の時間変化が求められ、この結果、高分子溶融体に特徴的な粘弾性や履歴効果を正確にシミュレーションすることが可能となったという。

特に、微視的な状態の移流を正しく取り扱うことが可能となり、任意の流動条件に対して適用可能なマルチスケール法の開発に成功したという。

今回の研究では、重要な2つの結果が得られたという。図2は、円筒状の障害物の周りを流れる高分子溶融体の速度分布、せん断速度分布、せん断応力分布とナビエ・ストークス方程式によって計算されたニュートン流体のものとの比較だが、両者を比較すると、速度分布、せん断速度分布に関しては上流側(図中の黒円に対して左側)と下流側(図中の黒円に対して右側)でほぼ対称な振る舞いをしているにも関わらず、せん断応力分布は両者で大きく異なる振る舞いをしており、高分子溶融体は上流側と下流側で非対称になっていることが分かる。

円筒状の障害物の周りの流れ。(a)高分子溶融体の挙動。(b)ニュートン流体の挙動。上から流速分布、せん断速度分布、せん断応力分布。図の値は適当な量で無次元化された絶対値を表す。(a)図中の赤色の線で囲われた内側の領域は、局所的なワイゼンベルグ数Wiの値が1より大きい領域。黄色の矢印は流れの方向を示している

これは高分子の持つ粘弾性に起因する履歴効果により応力の応答が遅れるためで、その結果、非対称な分布が生じている。また、円筒状の障害物の近くではせん断速度の値に比べて、せん断応力の値が小さな値を示している。これはシアシニングと呼ばれる現象で、局所的に粘度が低下していることを示している。これも高分子流体に特有の現象であり、流動による高分子鎖の配向が生じた結果起こる現象である。単純なせん断流れの中ではワイゼンベルグ数Wiが1より大きい場合にシアシニングが生じるが、今回のような複雑な流れの中では局所的なワイゼンベルグ数Wiが1より小さい領域でもシアシニングが生じていることが分かるが、これは履歴効果のためにシアシニングが遅れて出現しているためである。

また、図3は各場所における平均的な高分子同士の絡み合い量の空間分布を示したもので、円筒状の物体の周囲では高分子の絡み合いが少しほぐれて、絡み合いが減った領域が下流側に広がっていることが分かる。このような巨視的な流動中の微視的な物理量の空間分布は、従来の計算手法や現在の最新の計測技術でも測定することが困難な情報で、同シミュレーション法の計算によって、世界で初めて定量的に示すことができたという。

流動中の絡み合い量の空間分布。高分子同士の絡み合いが流動中にどのように変化するかを調べた図。円筒の周囲とその下流側(円筒の右側)で絡み合いが弱まっているのが分かる。黄色の矢印は流れの方向を示している

今回開発されたマルチスケールシミュレーション法により、目に見える巨視的な流動と原子・分子レベルの微視的な配置や形状と、絡み合いの状態を直接関係づけて計算することが可能となった結果、これまで計算手法では知ることができなかった高分子同士の絡み合い量の空間分布を知ることができるようになった。

また、任意の流動条件に対して適用可能であり、実験では調べることが困難な微視的な情報を解析できるため、新機能性材料の開発や物質科学のさまざまな分野に応用できるものと期待されるという。

なお、同シミュレーション法は、通常の分子動力学法とは異なり、膨大なCPUを持つコンピュータとの相性が高いという。これは、計算負荷の高い分子動力学法による計算部分が各計算点で独立であるため、計算の並列実行が容易であることに起因しており、そのため日本の次世代スーパーコンピュータ「京」のような膨大なCPUを持つ計算機システムを用いることで、分子数が100億個に達する世界最大規模の分子動力学シミュレーションも実現できる可能性があるとしており、それらを活用することで、今後、新機能性材料の開発や物質科学のさまざまな分野への応用を図ることが期待できると研究チームでは説明している。