（様式 20）成果報告書

TSUBAME 共同利用 平成 30 年度 産業利用 成果報告書

利用課題名 フィラー分散ポリマー複合材料の相分離構造シミュレーション
英文：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/
3 種類（粒子、ファイバー、板）の形状のフィラーを存在させた高分子溶融体における高分子マトリックスの相分離構造を高分子の SCF 法で MPI+GPU 計算可能とするプログラムを開発した。また、粒子形状のフィラーであれば、粒子形状を界面に濃度勾配をもった形状を固定した液体粒子（Diffuse Interface Particle: DIP）として扱い、粒子位置と相分離構造の同時最適化を可能とした。さらに DIP 粒子を粗視化ポテンシャルで結合し、フレキシブルなファイバー状のフィラーとして扱う実装もした。開発したプログラムは SOBA(Soft BlendsAnalyzer)として公開を開始した。
We developed a program that enables MPI + GPU calculation of the phase separated structure ofpolymer matrix in which three type fillers (particles, fibers, and plates) exist with the SCF methodof polymers. We also developed a method to simulate simultaneous optimization of filler’s positions
and phase separated structure, for particle fillers of which the interfaces have concentration
gradients (Diffuse Interface Particle: DIP). In addition, the DIP particles can be combined by usingcoarse graining potential and treated as a flexible fibrous filler. The developed program has beenreleased as Soft Blends Analyzer (SOBA).

Keywords: SCF theory, polymer composite, particle filler, fiber filler, plate filler, diffuse interface
particle

背景と目的
高分子材料は、様々なフィラーが添加され、その物Blends Analyzer(SOBA)[3]の改良をすることを目的とする。

性が改良されて用いられることが多い、例えばカーボンSCF では高分子のコンフォメーションのエントロピーブラックやシリカが添加されたタイヤ材料がその代表例の効果を精度よく計算に取り込めることができる。よっである。最近は CNT を添加し電気伝導性・熱伝導性にて、フィラー周囲の高分子マトリックスの相分離構造を優れた材料が開発されており、今後も様々な高性能材高精度に計算することができるが、計算するスケールと料の創生が期待されている。しかし、フィラーが添加さしては粗視化した高分子鎖が認識できるようなミクロ・メれた高分子材料は、フィラー高分子間の相互作用や高ゾスケールとなる。しかし、このスケールは複合材料に分子の絡み合い・分子量・自己組織化の影響等によりするために添加されるフィラーに比べると小さい。よっ様々な相分離構造をとる。よって、フィラー・高分子複合て、フィラーが存在する系の計算を SCF 法のスケール材料の構造解析には時間を要し、開発期間は長時間とに合わせると、領域分割ができる MPI ライブラリを利なる場合が多い。

用した大規模並列計算が必要となる。さらに、SCF 計

本課題においては、高性能な高分子複合材料の開発期間短縮のため、高分子の SCF 法（密度汎関数法算の高速性を実現するためには GPU の利用に対応する必要もある。の一つ：以降では簡単のため SCF 法と記載する）[1,2]

本課題においては、高分子の SCF 法の MPI+GPUを利用し、フィラーの凝集構造と高分子マトリックスの相並列計算における多くのバグを取り計算の安定化をは分離構造を同時にシミュレートできるプログラム Softかり、さらに MPI+GPU 並列計算において、領域を跨（様式 20）成果報告書
ぐフィラー存在下においても問題なく計算が実行できる

開発中の SOBA ではこの計算を全て有限差分法

ように境界問題とフィラー取り扱いにおけるマルチスケ（FDM）で行う。フィラーが存在する場合は、経路積分ール問題を解決することに注力し、高分子複合材料系の初期値（高分子末端の存在確率）が、フィラーがメッでの高分子マトリックスの SCF 計算を実行可能とする

シュに占める体積分率で影響を受ける。それを SOBA成果を得た。

で利用するスキームで表すと次のようになる[4]。
Q(0, ｒ𝐼 ; 0, ｒ𝐼 ) =概要
１．理論：フィラーを充填した系の高分子の SCF 法

1𝑣𝑓(3)

ここで、𝑣𝑓 は高分子末端が存在するメッシュ点での自由高分子を粗視化し、いくつかのモノマーをひとつのセ体積である。この自由体積を求めるため、SOBA ではグメントという粗視化単位として扱うと、セグメント間のフィラー周囲のみメッシュを細分化し、数値積分して𝑣𝑓排除体積効果が無視できるスケールとなり、高分子鎖を求める。の一端を空間のある点に固定すると、それにつながる

２．理論：フィラーの位置と相分離構造の自動最適化

高分子鎖を形成するセグメントの分布は、いわゆる酔前述の SCF 法において、𝑣𝑓 を求める方法は、メッシュ歩の軌跡と同等となり Gauss 分布となる。高分子のの細分化に計算コストがかかるので、フィラー位置は固SCF 法はこのような Gauss 鎖近似が成り立つスケー定とするのが現実的な制約となる。そこで、フィラーがルでの理論である。これを物理的に記述するために、存在する系の SCF 法の計算コストを削減するための鎖上のセグメントの位置を s、空間でのセグメントの位手法として導入したのが、フィラーの境界に濃度勾配が置を r で表し、高分子内の部分鎖の両端をそれぞれ s、ある溶液粒子（以降では Diffuse Interface Particle:r と s’、r’に固定したときのセグメントの統計的な重みDIP）とする方法である。SOBA ではフィラーの濃度分を 𝑄(𝑠, 𝐫; 𝑠 ′ 𝐫 ′ ) として表す。これは、以下の酔歩に相布∅ｐ (𝐫)を以下の式にてモデルして計算するようにした。

当する拡散方程式（Edwards 方程式と呼ぶ）に利用できる。
𝜕𝜕ｓ𝑄(𝑠, 𝐫; 𝑠 ′ 𝐫 ′ ) = [𝑏2𝑉(𝐫)−] 𝑄(𝑠, 𝐫; 𝑠 ′ 𝐫 ′ )
6𝑘𝐵 𝑇

(1)

∅ｐ (𝐫) =tanh{𝛼(𝑅−|𝐫−𝐫0 |)}+12(4)ここで、𝑅はフィラー粒子の半径、𝐫0 は粒子中止の位置、αは界面の幅を設定するパラメータである。この方法でここで、b はセグメントの有効長、T は温度、𝑉(𝐫)はセグ

は、フィラー界面に濃度勾配があるので、界面でも複数メント間相互作用と非圧縮条件からセグメントに働く外のメッシュ点が存在し、界面での物理量も計算可能とな部ポテンシャルエネルギーとする。以降では、この統計る。残念ながら高分子のコンフォメーションのエントロピ的な重みを経路積分と呼ぶこととする。経路積分を利ー効果を精密に計算するのを諦めることにはなるが、フ用すれば、ある点 r に集まるセグメントの濃度分布φ(r)ィラー位置を固定したまま複数回 SCF を回し、高分子は次の式で求まる。

フィラーへ働く化学ポテンシャルからフィラーへ働く力を計算し、力の方向に微小量フィラーを移動させた後、再𝑁∅(𝐫) = 𝐶 ∫ 𝑑𝑠 ∫ 𝑑𝒓0 ∫ 𝑑𝒓𝑁 𝑄(0, 𝒓0 ; 𝑠, 𝒓)𝑄(𝑠, 𝒓; 𝑁, 𝑟𝑁 )0(2)
ここで、C は計算で選択するアンサンブルから決まる規度 SCF を回す疑似的動力学をすることにより、DIP の位置と高分子マトリックスの相分離構造の同時最適化が可能となる。

格化定数、外側の積分は鎖に沿って全てのセグメントを積分する意味、内側の二つの積分は鎖の両末端位置を空間的に全て積分するという意味である。∅(𝐫)が結果および考察
１．SCF 法の MPI+GPU 計算の安定化

求まらないと𝑉(𝐫)を求めることができないので、この計最初に、流体力学効果を導入した動力学計算におい算 は 非 線 形 な 計 算 で あ り 自 己 無 撞 着 場 （ Self-て、問題によっては計算がアボートする不具合が発見

Consistent Field: SCF）計算を必要とする。

されたのでデバッグし、これまでのコードに存在したデ（様式 20）成果報告書
バイス上でのメモリの確保におけるバグを取り除いた。

２．3 種 類の 形状のフ ィラ ー存在系 の相分 離計算 のその結果、安定した計算が行えるようになった。さらに、

MPI+GPU 計算

セグが相互作用パラメータχを位置により可変とし温図 2 に示すのは境界を固体壁した球、ファイバー、板度勾配のある系を疑似的に設定できるような改良を施の 3 種類の形状のフィラーが存在する系における対称した。つまり、温度が高い領域はχが小さく溶融状態とジブロック・コポリマーの SCF 計算の結果の一例であなり、温度が低い領域はχが大きく相分離するという条る。このような計算を MPI+GPU 計算で実行すること件である。図１に示すのはジブロック。コポリマーのミクができるようにコードを改良することができた。フィラー

ロ相分離構造を 64（Core+GPU）MPI 並列で計算した形状は任意の大きさを設定することができるが、それに例である。図 1(a)に示すのは、温度勾配に相当するχ

は高分子マトリックスの SCF 計算とオブスタクルスケーを球対称で設定した場合のプロファイルであり、球の中ルの計算のマルチスケール計算を同時に実行するス

心部で相分離し、外周は高温で溶融状態を保ち、周期キームが必要である。当初、フィラーは大きいので、全境界条件の影響を排除した自由界面を存在させた計

MPI プロセスにて、全フィラーを扱う同等なスキームを算とした例である。図 1(b)に示すように自由界面を伴う実行する設計としていたが、256 プロセスを利用する大球状の相分離構造を形成させることに成功した。

規模 MPI＋GPU 並列計算とした場合、フィラーの取り

図 1 χパラメータに位置依存性を導入して自由界面を図 2 3 種類の形状のフィラーが存在する系にお形成させたジブロック・コポリマーの相分離構造

ける対称ジブロック・コポリマーの相分離構造

64（Core＋GPU）MPI 並列計算、Gyroid 構造の8FlatMPI 計算、システムサイズ 323(643 メッシュ）、鎖4096unit_cell(5123 メッシュ）、鎖長 N=20、ブロック比

長 N=20、χN=20、ブロック比 f=0.5、静的学計算、

f=0.35、動力学計算、(a)χパラメータに分布を導入、(a)セグメント濃度分布、(b)等値面表示、Φ=0.5、球、フ(b)等値面表示、流体力学を導入：粘度μ=0.001、等値

ァイバー、板の 3 種類のフィラーが存在する周囲で相面表示：Φ=0.75、時間 t =76

分離している。

（様式 20）成果報告書
扱いの効率の悪さが露見し計算速度が低下したのでスキームを改め、各分割小領域に存在するフィラーのみを取り扱うようスキームを変更した。このような分割スキ

1

𝑉𝑏 = 𝑘𝑏 （R 0 − |𝐫0 − 𝐫1 |）

2

2

1

𝑉𝑎 = 𝑘𝑎 (𝜃 − 𝜃0 )22

(5)
(6)

ームは、粒子を扱う MD の MPI 計算と同等と認識しがここで、𝑘𝑏 、𝑘𝑎 はそれぞれポテンシャルの強度を設定すちであるが、大規模粒子、長いファイバー、広面積板等る任意のパラメータ、R 0 は RIP 粒子の半径、𝐫0 と𝐫1 はボを導入した計算では一つのフィラーが複数の小領域をンドを形成する 2 つの DIP 粒子の中心位置、𝜃は２つ跨ぐ条件が発生するので一般的な粒子を扱う MD ののボンドがなす角、𝜃0 は任意の基準となる角度である。

MPI 計算以上にスキームに工夫が必要であり、プログ

図 4 に示すような大規模 MPI+GPU 計算により、多くラムを改良し各フィラーの存在領域の認識を効率よく行の統計量を得ることができ、パラメータの差異が構造にうようにした。

与える結果を評価することができた。その結果、フィラ

３．ポリマーブレンド中の DIP 粒子

ー周囲はポリマーマトリックスであり、真空中に置く粗図 3 に示すのは DIP 粒子を A/B ポリマーブレンドに視化 MD とは異なり、非圧縮条件（ポリマー溶融体なの投入した例である。粒子は A ポリマーを避けるような相で非圧縮な系として扱う）により受ける力は非常に大き互作用パラメータが設定されているので、粒子周囲を Bいことがわかった。また、DIP 粒子の界面幅を設定すポリマーが覆い、凝集構造を形成していることがわかる。

る(4)式におけるパラメータαの影響もあり、例えばαこのような計算も、領域分割した場合に境界に存在す=0.1、R 0 =2、の条件において、DIP 粒子間距離は、4.9

る粒子情報を通信し、MPI+GPU 計算可能とした。

が平均値となった。このような計算を 256（Gore+GPU）
MPI 並列まで実行可能なことを確認した。

図 3 ポリマーブレンド中の DIP 粒子
1（Core+GPGPU）計算、システムサイズ 323(643 メッシュ）、鎖長 N=10、A/B ブレンドに半径 2 の C 粒子64 個を投入し、χAC=1、とした場合の B セグメント分布の断面図。

４．DIP 粒子を粗視化ポテンシャルで接続した計算
図 4 に示すのは、CNT のようなフレキシブルな繊維状のファラーを計算するために、DIP 粒子の粗視化ポテンシャルで結合した計算例におけるフィラーの構造である。結合している DIP 粒子同志をシリンダーとして表示している。粗視化ポテンシャルとしては次のボンドとアングルポテンシャルを導入した。

図 4 DIP 粒子を粗視化ポテンシャルで結合しフレキシブルなファイバーを想定した計算例
64（Core+GPGPU）計算、システムサイズ 2563(5123メッシュ）、鎖長 N=10 の A ポリマー中に半径 2 のDIP を 32 個接続した模擬ファイバーを 2048 存在させた。α=0.1、kb=30、ka=75、θ0=180°とし直鎖状の配置から 20 回緩和した例。

６．FlatMPI と MIP+GPU 計算の比較現在実装中の SCF 法のコードでは、問題によって様々なスキームが作動するので、一概に計算速度を比較することは難しいが、フィラーが存在しない通常の高（様式 20）成果報告書
分子の SCF 計算を MPI+GPU 並列計算と FlatMPI

参考文献

計算で比較すると、NVLINK は利用するが最新の[1] 公益社団法人 新化学技術推進協会 (編集)、“増DirectGPU 等の手法を利用しないプログラムにおいて、補版 高分子材料シミュレーション: OCTA活用事例

システムサイズ 137.63（ミクロ相分離の Gyroid 構造集”、 化学工業日報社 (2017).

512 unit cell）、メッシュサイズ 2563、鎖長 N=20、ブ[2] T. Honda and T. Kawakatsu, “Applications of

ロック比 f=0.35、鎖の数約 13 万本の計算条件にてSUSHI in OCTA system” in Nanostructured

8(Core+GPU)MPI 並列計算での比較結果は、

Soft1 GPU（P100）≒ 72 core Intel Xeon E5-2697A v4 @ 2.60GHz相当であるとの結果となった。この値は最新の CUDAのプログラミング手法を利用してチューニングすればさらに MPI+GPU 計算の速度が向上すると考えられる。
また、大規模 PC クラスターではなくとも数枚の GPUで十分な高速計算を実現できるといえる。

まとめ、今後の課題
高分子の SCF 法に改良を施し、3 種類のフィラー（球、ファイバー、板）が存在する系においても高分子マトリックスの相分離構造の計算ができるようになった。さらに、
DIP 粒子を導入すれば、粒子の位置と相分離構造の同時最適化が可能となった。また、DIP 粒子間に粗視化ポテンシャルを導入し、フレキシブルなフィラーを模擬した 計算も可能とした 。これらのすべての計算はMPI+GPU 計算可能であり、256（Core+GPGP）並列においても並列化効率に問題のない計算を実行できることを確認した。
開発したプログラムは、 http://octa.jp/ より SOBA
Ver.1.や SUSH10.7 として 2019 年 4 月より公開した。
無料でダウンロードが可能である。
今後は、MPI+GPU 計算に粘弾性効果を導入すること、DIP 粒子のパラメータ・サーベイを継続し複合材料を評価できる具体性のあるシミュレーション技術に発展させること、また、CUDA に導入された最新の機能によりより並列計算の速度を向上させることが必要である。

謝辞
この成果は、国立研究開発法人新エネルギー・産業技術総合開発機構（NEDO）の委託業務（P16010）の結果得られたものです。

Matter:Experiments,

TheoryandPerspectives, A. V. Zvelindovsky, ed., SpringerVerlag (2007).
[3]http://octa.jp/ にて公開中（アカウント取得と loginが必要であるが利用は無料） (2019).
[4] H. Morita, M. Toda and T. Honda, Polymer

Journal, 48, 451 (2016).