第 5 章 分子動力学法を用いた引張解析
5.2 解析手法
本研究ではモデルを作成する際に高分子鎖の化学構造の特徴を保持する範囲 で,数個の原子を1 つのユニットとして扱う United Atom モデルを用いて分子 動力学シミュレーションを行なった.これにより (CH), (CH3)を1つの粒子とし て見なすことができ,計算コストを削減できる.まず分子鎖内ボンド結合 (bond, angle, torsion)および非ボンド相互作用(van der Waals )に対するポテンシャルを式 (5.1)で評価した.[5-2]
𝑈
𝑝𝑜𝑡𝑒𝑛𝑡𝑖𝑎𝑙= 𝑈
𝑏𝑜𝑛𝑑+ 𝑈
𝑎𝑛𝑔𝑙𝑒+ 𝑈
𝑡𝑜𝑟𝑠𝑖𝑜𝑛+ 𝑈
𝑛𝑜𝑛−𝑏𝑜𝑛𝑑𝑖𝑛𝑔・・・ (5.1)
また力場パラメーターにはDREIDING力場を用いた.DREIDINGは,各原子の 混成軌道の種類に基づき,中心原子のみでパラメーター値を決定できる力場で あるため,適用範囲が非常に広く,新しい原子にも拡張が容易な汎用ポテンシ ャルセットである.式(5.1)の右辺の各項はそれぞれ分子鎖内の粒子間の結合長 r, 結合角θ, 2面角φ, 非結合原子距離r’に対するポテンシャルを表し,𝑈𝑏𝑜𝑛𝑑は 結合伸縮ポテンシャル,𝑈𝑎𝑛𝑔𝑙𝑒は結合変角ポテンシャル,𝑈𝑡𝑜𝑟𝑠𝑖𝑜𝑛は結合二面角
126
ポテンシャル,𝑈𝑛𝑜𝑛−𝑏𝑜𝑛𝑑𝑖𝑛𝑔は非結合相互作用を示している.それぞれ式 (5.2) ~ (5.4)で表される.[5-2]
𝑈
𝑏𝑜𝑛𝑑(𝑟) =
12
𝑘
𝑟(𝑟 − 𝑟
0)
2・・・ (5.2) 𝑈
𝑎𝑛𝑔𝑙𝑒(𝜃) =
12
𝑘
𝜃(𝜃 − 𝜃
0)
2・・・ (5.2) 𝑈
𝑡𝑜𝑟𝑠𝑖𝑜𝑛(𝜑) = 𝑘 ∑
𝑁−1𝑛=0𝐴
𝑛𝑐𝑜𝑠
𝑛𝜑 ・・・(5.3)
ここでkr, kθ, kはバネ定数, r0は結合平衡長,θ0は平衡角,N-1は多項式の次数,
Anは原子種の対に対する定数,σは長さのパラメーターでファンデルワールス半 径,ε はエネルギーのパラメーターを表し,rc は cut off (カットオフ距離) を表しており,それぞれの値はTable 5.2.1に示す.ここで,cut offは原 子間の相互作用を無視することのできる原子間距離であり,これを設定するこ とにより計算コストの削減が可能になる.また分子動力学シミュレーションに おいて系の温度Tは式 (5.5)のように表される.
𝑇 =
(3𝑁−𝑁1𝑐)𝑘𝐵
∑
𝑝𝑖2𝑚𝑖
𝑁𝑖=1
・・・ (5.5)
Nは原子数,Ncは拘束される自由度数,kBはボルツマン定数,pi, miは各々,原 子iの運動量と質量を表している.またNcは3次元周期境界条件の時にNc=3と なる.さらに本研究では周期境界条件を用いているため圧力テンソル p はビリ アル定理より式 (5.6)のように表される.
127
𝐩 =
𝟏𝑉(∑
𝑚𝑝𝒊𝟐𝒊
𝑁𝑖=0
+ ∑
𝑁−1𝑖=0∑
𝑁𝑗>𝑖𝐫
𝒊𝒋𝐟
𝒊𝒋) ・・・ (5.6)
式 (5.6)においてVはセルの体積,𝐫𝒊𝒋は原子i, j間ベクトル𝐫𝒊-𝐫𝒋, 𝐟𝒊𝒋は原子jが 原子iに及ぼす力を表している.
5.2.1 モデル作成
モデル作成および解析には株式会社 JSOL 製の材料物性解析ソフトウェア
J-OCTA 2.0を用いた.本研究ではモノマーとしてFigure 5.2.1のようなL乳酸を
作成した.DREIDING 力場に設定した後,構造最適化計算を行った.作成した モノマーモデルの繰り返し数(重合度)や立体構造を設定し,任意の重合度のポリ マーモデルを作成した.本研究ではPLLAの数平均分子量Mn=8000, 16000, 24000 のポリマーモデルの作成を行った.分子の立体構造は180°のモノマー間のねじ れ角を持つアイソタクチック構造とした.作成された一本鎖モデルは Figure
5.2.2 に示す.一本鎖モデル作成後はモノマーモデル作成時と同様の計算を行っ
た.その後,作成したポリマーモデルをユニットセル内に50本ランダムな初期 位置に配置した.その後,外圧を大気圧(約 0.1MPa)に設定し,アンサンブルは 粒子数,圧力,温度一定のNPTを選択し,単位時間ステップΔt=10-15 [s] とし,
5万ステップで緩和計算を行った.圧縮後,セル寸法は約10 nm×10 nm×10 nm となり,セル内密度は約1.29 [g/cm3]となった.その後,粒子数,体積,温度一 定のNVTアンサンブルで100万ステップの緩和計算を行い,セル内を安定状態 にした.作成したモデルをFigure 5.2.3に示す.各モデルはそれぞれ3つ作成し た.
5.2.2 引張解析
5.2.1 で作成したモデルを用いて引張解析を行った.引張解析はセルの z 方向
の全長Lzに対して両端約 10%の原子を固定して行った.固定をする際に,片側
10%には速度を与えずに固定し,反対側の10%には z 軸正方向に速度 Vzを与え
て固定することにより,Figure 5.2.4のように引張試験を模擬した.境界条件に は周期境界条件を適用した.引張速度は1×1010 [1/s]とした.モデルのセル内の 応力とひずみの変動を解析結果から抽出し,応力ひずみ線図を得ることが出来
128
る.