（様式 20）成果報告書

TSUBAME 共同利用 令和 2 年度 産業利用 成果報告書

利用課題名 フィラー分散ポリマー複合材料の相分離構造シミュレーション
英文：Simulation of phase separated structures of polymeric materials with fillers
本田 隆
Takashi Honda
先端素材高速開発技術研究組合
Research Association of High-Throughput Design and Development for
Advanced Functional Materials (ADMAT)
http://www.admat.or.jp/
粒子フィラーと高分子混合系における相分離構造を計算するために、高分子の SCF 理論において、粒子フィラーを球形状に形状を固定された溶媒である DIP(Diffuse Interface Particle）として取り扱い MPI+GPGPU 大規模計算が可能なプログラム SOBA(Soft Blend Analyzer)の開発を継続した。新たな改良点として、フィラーの配置がローカルミニマムな構造に陥らないように計算途中で移動させる手法をいくつか導入した。
SOBA (Soft Blend Analyzer) is developed to calculate the phase separation structure of particle filler
and polymer mixing system and use MPI + GPU computers. The SCF theory of the polymeric materials
is introduced, and the particle filler is treated as DIP (Diffuse Interface Particle), which is a solvent
whose shape is fixed in a spherical shape. To avoid local minimum structure, we introduced several
methods of moving the filler in the SCF calculation.

Keywords: SCF theory, polymer composite, particle filler, diffuse 𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 particle, phaseseparation
高分子材料は、充填材とブレンドされ利用される場Interface Particle: DIP と記す）を利用した高分子の合が多く、タイヤ用途ではエラストマーにカーボンブラッSCF 計算を用いた。DIP の形状は次のモデル式で記クやシリカが添加され強度や低燃費性の改良がされる。

述される。

また、電気伝導性・熱伝導性を付与する高分子材料に

Φ𝑝 (ｒ) =tanh{𝛼(𝑅−|𝐫−𝐫0 |)}+12(1)は、CNT を添加することも行われている。これらの材料

ここで、Φ𝑝 (ｒ)は位置ｒにおける DIP 粒子の濃度分布、物性は、フィラーの分散性に影響される。しかし、フィラ𝑅はフィラー粒子の半径、𝐫0 は粒子中止の位置、αはーと高分子複合材料の構造解析や評価は難しく、材料界面の幅を設定するパラメータである。の開発期間は長時間となる場合が多い。
本課題では、高性能な高分子複合材料の開発期間

フィラーを充填した系の高分子の SCF 法の説明は、先の報告書等[4、5]に記載したので詳細は割愛する。

短縮のため、高分子の SCF 法（密度汎関数法の一つ：

２． フィラー間の非圧縮条件の導入

以降では簡単のため SCF 法と記載する）[1,2]を利用し、DIP が重ならない効果は、SCF 計算における非圧縮フィラーの凝集構造と高分子マトリックスの相分離構造条件にて自動的に取り入れるとしてきたが、DIP 界面を同時にシミュレートできるプログラム Soft Blends

幅が広い場合、非圧縮条件のみではある程度 DIP がAnalyzer(SOBA)[3]の改良を行ってきた。本年度は、重なる場合が発生した。そこで、フィラー間の相互作用フィラー凝集構造を計算する方法を改良し、かつ並列エネルギーを導入することにした。これは次式で定義さ化効率についてまとめた。

れる。

概要

𝐹𝑝𝑝 (ｒ) = 𝜒𝑃𝑃 Φ𝑝 (ｒ)Φ𝑃 (ｒ)(2)

１． 理論：フィラーを充填した系の高分子の SCF 法

ここで、𝜒𝑃𝑃 はフィラー用の相互作用パラメータである。

フィラーの凝集構造をシミュレートするためにフィラー

３． 粒子形状の導入方法の改良の境界に濃度勾配がある溶液粒子（以降では Diffuse

ローカルミニマムな相分離構造にトラップされるのを（様式 20）成果報告書防ぐ方法として、計算初期に DIP の粒子形状を導入せず相分離を先行させ、あとから粒子形状を導入する方法が 有効であることがわかった。
４． フィラーの移動方法の改良
DIP を高分子マトリックス内で移動させるのに、以下のような簡易な動力学モデルを利用していた。ている、ｆ(𝑟)は SCF 法で得られた力であり、この式は粒Figure 1. Phase separated structure of DIPs in an
A/B polymer blend
System size is 323, 𝜒𝑃𝑃 =106, two-dimensional
system and periodic boundary conditions wereintroduced. (a) initial structure, (b) after 50 step
without the shape of DIPs, (c) after 100 step with theshape of DIPs.

子を働く力ｆ(𝑟)の方向に微小移動させることを意味する。

うに 2 次元の小さい系での 4 FlatMPI 計算の例であこのモデルは高分子マトリックスが高粘度で慣性力がる。(a)：計算初期の無秩序状態、(b)：DIP 粒子の形状無視できる場合に導入できる。しかし、このモデルではを導入しない相分離の先行による相分離構造、SCF 計算の収束性がよい場合、粒子はあまり移動せ(c):DIP の形状を導入した後の相分離構造を示す。(b)ずローカルミニマムな相分離構造に収束する場合が多の状態では相分離構造は粒子形状の影響を受けないいことが判明した。そこで、次式による疑似的な粒子のので滑らかな界面形状を有している。その後粒子形状慣性を導入した。

を導入するので最終的な構造(c)では相分離界面に粒

𝑑𝐫0＝𝐷𝑒 ｆ(ｒ)𝑑𝑡(3)

ここで、𝐷𝑒 は疑似的な易動度で時間の効果も考慮され

𝑑𝐫0 𝑖+1𝑖+1＝𝑣 𝑖+1 (ｒ)＝𝑣 𝑖 (ｒ) + 𝐷𝑒 ｆ (ｒ) − 𝜂𝑣 𝑖 (ｒ)𝑑𝑡(4)子が存在し、界面は粒子形状の影響を受け、粒子配置も密な六方配置に変化している。

ここで、上添え字 i はスキームの番号を意味し𝑣 𝑖 は i 番目のスキームにおける疑似的な速度であり、右辺の式
において第 1 項の i 番目の速度は第 2 項の i+1 番目の力の効果が加えられ i+1 番目の速度に更新される。
しかし、それだけでは力の効果が積分されて粒子の移動が止まらなくなるので、第 3 項による疑似的な粘性係数η（>0）を導入して移動を抑制する。このηは初期値0 とし、スキーム毎に増加させ、最終的に 1 とし、式（4）を最終的には式（3）と同等とする。
結果および考察１．モデル系
鎖長 N=10 のポリマーA/B 溶融系(組成 0.5:0.5)に

R=４、α=3、𝐷𝑒 =0.0015～0.015 の DIP フィラーを添加した場合をベンチマーク用のモデルとして評価した。
DIP は 100SCF 毎に移動させた。大規模系でのフィラーの初期配置は、最初に面心立方格子状の位置座標を発生し、モンテカルロ法により空いている位置座標に

Figure 2. Phase separated structure of DIPs in an
A/B polymer blend
System size 1283, mesh size 2563, 𝜒𝑃𝑃 =106, periodic
boundary conditions were introduced. The numberof DIPs is 782 and the volume fraction of the DIPs is0.1. χAB=χAP=0.4.

３．フィラーの移動方法の改良

粒子を配置する方法をとった。

図 2 に示すのは、粒子形状を無視する方法と疑似的

２．粒子形状の導入方法の改良な慣性を導入した方法を組み合わせて求めた相分離

図１に示すのは粒子形状の導入方法を変化させた構造の一例であり、A/B ポリマーの相分離界面の等値場合の相分離構造の変化であり、状態変化がわかるよ面と DIP 粒子の配置を図示した。最初に、粒子形状を（様式 20）成果報告書無視する方法＋疑似的な慣性を導入した方法による計算を 50 ステップ行い(𝜂 = 1 − 0.999𝑖 )、次に粒子形状を導入する計算のみを 50 ステップ行った結果である。
構造はマクロ相分離であり、DIP 粒子との相互作用が無い B セグメントのドメイン内に DIP 粒子が移動していることがわかる。図 3 に示すのは、図 2 と同じ条件でA/B 間の相互作用パラメータχAB と A セグメント/DIP間の相互作用パラメータχAP を同等なχ＝χAB＝χAPとして変化させた場合のフィラーの 2 体相関関数の比較である。 面心立方格子状の位置は充填構造として安定位置に近いが、体積分率が 0.1 としてあるので、χが大きくなると DIP は凝集し、第一ピークが低 r 側にわずかに移動することがわかる。

Figure 4. Comparison of calculation times [sec] forvarious system sizesする）、および８（MPI＋GPU）計算での計算時間を図４に示す。CPU 数は横軸の値（基本ユニットセルの倍数）と同じである。メッシュサイズが 64 倍未満では１GPU の計算が早く、その後は MPI 計算が高速となる。メッシュサイズが 512 倍では、512CPU のFlatMPI 計算が高速であるが、8GPU を利用しても3 倍程度の計算時間の増加で済む。
次に、GPU を利用する場合、できるだけ GPU のメモリを利用するジョブの方がコストパフォーマンスに優れるので、１GPU のみでの並列化効率βを測Figure 3. χ parameter dependence of pair
distribution function of DIPs
The condition is almost same as Figure2 but χ=χ
AB=χAP are changed.

４．並列化効率定した結果を図 5 に示す。βは次の式で定義される。
𝛽=𝑇１𝑇𝑛

ここで、𝑇１ は基本ユニットでの計算時間、𝑇𝑛 はシステム

本課題で利用した SOBA の SCF 計算の並列化効率を計測したのでそれも報告する。標準とした問題は、高分子のミクロ相分離構造である Gyroid を計算する SCF 計算であり、この問題では、経路積分（場の exp 演算や四則演算を伴うステンシル計算）等を 100×2×100 回行う演算を行う。システムサイ
ズと並列数を変化させ計算効率を計測した。システ
ムは全て周期境界条件とした。基本ユニットセルのメッシュサイズは 323 とした。利用した GPU は P100,
CPU は Xeon E5-2697A V4 である。
1GPU と８GPU での計算効率システムサイズを変化させた場合の１GPU、CPUのみの FlatMPI（ここでの CPU 数はコア数を意味

Figure 5. Calculation efficiency using 1 GPU“b” means the block size of CUDA calculation.

（様式 20）成果報告書サイズを n 倍した場合の計算時間である。システムサ

システムサイズを大きくしても計算効率は 50％程度

イズを大きくすると計算効率は低下するが、216 倍で停滞するということで、MPI＋GPU 並列の効果が（メッシュサイズ 1923）のシステムサイズにしても効大規模系でも発揮できているといえる。

率は 2.6%程度を維持しており、理論的な推算値

まとめ、今後の課題

1/216=0.46%より 5 倍以上効率がよく、CUDA コアフィラーが分散した高分子溶融系の相分離構造のロの並列計算により高速化されていることがわかる。

ーカルミニマムな構造を避ける手法を開発することがでまた、CUDA 計算におけるブロック数を変化させてきた。この手法は、今後高分子複合材料の相構造を予みたが、その影響は無視できる程度であることも分測する上で有用な技術となると考える。また、DIP 粒子かった（点はほとんど重なっている）。

間に粗視化ポテンシャルを与えることにより繊維状のフMPI＋GPU での並列化効率

ィラーを含む系の相分離構造の予測も可能となる。

GPU のメモリを多く利用する計算を行わないと計開発したプログラムは、 http://octa.jp/ より SOBA

算のコストパフォーマンスはよくならいので、基本ユVer.3.や SUSH11.1 として公開予定である。

ニットセルのメッシュサイズを 1923 とし、n(MPI＋

今後は、様々な大きさや形状のフィラーを含む複合材GPU)計算における並列化効率βを計測した結果を料や、個体壁表面への密着性等の様々な高分子複合

図 6 に示す。この場合βはウィーク・スケーリングの材料のアプリケーションの例を増やしてゆきたい。さら計算効率となる。図 6 中に示すキーの右脇の数字に、粒子の動力学における GPU を利用した並列計算

xyz は各方向へ拡張した領域の数でシステムのトポの改良や、SCF 法の高速化にも取り組んでゆきたい。

ロジーの違いを意味する。xyz の積が MPI プロセス数 n になる。最大のシステムサイズは xyz=666 で基謝辞

本ユニットの 216 倍である。この大きさでも計算効率

この成果は、国⽴研究開発法⼈新エネルギー・産業技は 基 本 ユ ニ ッ ト の 53 ％ 程 度 で あ る 。

術総合開発機構（NEDO）の委託業（JPNP16010）

xyz=122,212,221 の場合効率が低下するが、これの結果得られたものです。

は１ノード内（４GPU）の計算になる。何等かの計算

参考文献

効率の低下が発生しているようであるが理由は不明[1] 公益社団法人 新化学技術推進協会 (編集)、“増である。図 6 の Y 軸は各点の差異がわかるようにス

補版 高分子材料シミュレーション: OCTA 活用事例集”、ケールを拡大してあるが、この計測でわかったことは、

化学工業日報社 (2017).
[2] T. Honda and T. Kawakatsu, “Applications ofSUSHI in OCTA system” in Nanostructured Soft
Matter: Experiments, Theory and Perspectives, A.
V. Zvelindovsky, ed., Springer-Verlag (2007).
[3] http://octa.jp/ にて公開中（アカウント取得と loginが必要であるが利用は無料） (2020).
[4] 本田 隆、“フィラー分散ポリマー複合材料の相分離構造シミュレーション”、TSUBAME 共同利用 平成
30 年度 産業利用 成果報告書[5] H. Morita, M. Toda and T. Honda, Polymer

Journal, 48, 451 (2016).

Figure 6. Calculation efficiency of weak scaling