セルロース結晶の第一原理計算とPLA/CNF複合材料の相分離シミュレーション#

Advance/PHASE Logo

植物由来のセルロースナノファイバー（CNF）は、軽量かつ高強度な次世代のグリーンコンポジット材料として期待されており、特にポリ乳酸（PLA）などの生分解性プラスチックとの複合化が盛んに研究されています。複合材料の力学特性を最大限に引き出すためには、柔らかい樹脂マトリックス中に硬いCNFをどのように分散させ、強固なネットワーク（パーコレーション）を形成させるかが鍵となります。本事例では、第一原理計算ソフトウェアAdvance/PHASEを用いて算出したセルロース結晶の弾性率を、メゾスケールの「フェーズフィールド（PF）法」にブリッジするマルチスケール連携手法を紹介します。大きな力学物性のミスマッチが引き起こす「弾性相分離」のメカニズムを定式化し、実際の材料設計に直結する異方的なネットワーク構造の形成プロセスを明らかにしました。

Keywords: 第一原理計算（DFT）, セルロースI $\beta$ , ポリ乳酸（PLA）, フェーズフィールド法, マルチスケール連携, 弾性相分離, Cahn-Hilliard方程式

セルロース結晶の構造の特徴と第一原理計算#

1. セルロースI $\beta$ の特異性と計算の難しさ#

天然セルロースの主成分であるセルロースI $\beta$ 結晶 [1] には、分子鎖方向の強固な共有結合、分子鎖間の方向性を持つ水素結合、そしてシート間に働く広範囲のファンデルワールス（vdW）分散力という、スケールも性質も全く異なる3つの相互作用が混在しています。これは第一原理計算にとって非常に難易度の高い対象系となっています。Liら [2] の報告にもある通り、標準的なGGA（一般化勾配近似）などの汎関数ではシート間引力を表現できず、計算中にセルが破綻しやすいという課題がありました。

2. 異方性を考慮した弾性率のDFT計算#

この難易度の高いセルロースI $\beta$ 系に対し、第一原理計算ソフトウェアAdvance/PHASEを使用してDFT計算を行いました。その際、vdW相互作用に対して、分散力補正（DFT-D2）を適用しました。さらに、等方的な体積変化を強制すると強固な分子鎖（ $c$ 軸）が無理に圧縮され弾性率が過大評価（約40 GPa）される問題を回避するため、「 $c$ 軸を固定し、 $a$ 軸・ $b$ 軸の比率（ $b/a$ ）を各体積で最適化する」という2次元的かつ異方的なアプローチを採用しました。

この手法により得られた状態方程式（EOS）へのフィッティング結果から、最適体積 $V_0 = 639.0$ Å³ と、体積弾性率 $B_0 = 19.1$ GPa が算出されました（図1）。平衡体積 $V_0 \approx 639$ Å³ と格子定数のスケーリング（a軸: x0.9777, b軸: x0.9939）は、Driら [3] の計算結果と高精度に一致しています。体積弾性率（2次元、横方向）の値は、実験的に得られている横方向弾性率（約20 GPa）[4] や、Driらによる第一原理計算の横方向ヤング率（ $a$ 軸：19 GPa）ともよく一致しています。

セルロースIβの状態方程式フィッティング

図1. セルロースI$\beta$の2次元異方性緩和を伴う状態方程式（EOS）フィッティング結果（体積弾性率 19.06 GPa）。

3. 等高線マップから読み解く分子間相互作用の裏付け#

図2は、体積と $b/a$ 比の2変数に対するポテンシャルエネルギー曲面（等高線マップ）です。平坦なエネルギー谷の中に明確なGlobal Minimum（赤い星印）が確認できます。注目すべきは、系を圧縮（体積を減少）させるにつれて、最適となる $b/a$ 比が大きくなる（上へシフトする）点です。

文献 [4] での議論と照らし合わせると、セルロースの $b$ 軸方向は強固な水素結合ネットワークで結ばれている一方、 $a$ 軸方向（シート間）の凝集は主に柔らかな分散力によって担われています。このため、圧縮時には「剛性の低い分散力方向（ $a$ 軸）が優先的に縮み、結果として $b/a$ 比が上昇する」という物理的に妥当な変形メカニズムが、この等高線の傾きに明確に現れています。得られたセルロース結晶の構造を図3に示します。

セルロースIβのエネルギー等高線マップ

図2. 体積とb/a比に対するエネルギー等高線マップ。赤い星印がGlobal Minimumを示す。

セルロースIβの結晶構造

図3. DFT-D2を用いて最適化されたセルロースI$\beta$結晶の構造。左：単位胞（84原子）、右：2x2x2スーパーセル。

フェーズフィールド法による相分離シミュレーション#

1. 解析モデルと理論背景（熱力学パラメータの導入）#

ポリ乳酸/セルロースナノファイバー（PLA/CNF）ブレンドのメゾスケール構造形成を追跡するため、系の全自由エネルギー $F$ を以下の式で定義しました。

$F = \int \left[ f_\text{chem}(c) + \frac{\kappa}{2}|\nabla c|^2 + f_\text{el}(c) \right] dV$

ここで、 $c$ はCNFの体積分率（濃度場、0〜1）です。化学的自由エネルギー $f_\text{chem}$ は系の相分離を駆動するベースであり、Flory-Huggins理論に基づく $\chi$ （カイ）パラメータによって決定されます。本解析では迅速な材料設計スキームを実証するため、Hansen溶解度パラメータ（HSP）の文献値 [5] を用いて $\chi$ パラメータを理論的に導出しました。

$\chi \approx \frac{V_\text{ref}}{RT} \left[ (\delta_{d,1} - \delta_{d,2})^2 + (\delta_{p,1} - \delta_{p,2})^2 + (\delta_{h,1} - \delta_{h,2})^2 \right]$

PLAのHSP値（ $\delta_d \approx 18.6$ , $\delta_p \approx 9.9$ , $\delta_h \approx 6.0$ ）とCNFのHSP値（ $\delta_d \approx 16.1$ , $\delta_p \approx 14.9$ , $\delta_h \approx 23.3$ ）の差、特に水素結合項（ $\delta_h$ ）の大きなミスマッチから、系が強い相分離傾向を持つことが裏付けられます。数値計算上は、安定性を確保するため二重井戸型ポテンシャル $f_\text{chem} = W c^2 (1-c)^2$ を適用し、駆動力 $W$ に導出した $\chi$ を反映させています。

2. 弾性ひずみエネルギーの導入（DFT計算との連携）#

次に、本事例のハイライトである弾性ひずみエネルギー $f_\text{el}(c)$ を導入します。ここで第一原理計算（DFT）の役割は、得られた弾性率の絶対値（19.1 GPa）を直接フェーズフィールド法に代入することではなく、樹脂（約3 GPa）とCNF（19.1 GPa）の弾性率比（約6倍）を正確に評価することにあります。この比をベースに、フェーズフィールド法におけるエネルギーのスケール（化学的駆動力 $W$ など）に合わせて規格化した値として、弾性効果係数0.8を設定しました。

このような大きな弾性非均質性（Elastic Inhomogeneity）が存在する場合、凝集ドメインは系全体のひずみエネルギーを最小化しようとして特定の結晶学的方位へ引き延ばされます。

計算コストを大幅に抑えるため、Khachaturyanの微視的弾性理論 [6] に基づく「異方的な弾性エネルギーペナルティ」をフーリエ空間（波数ベクトル $k$ ）で直接導入する手法を採用しました。弾性ポテンシャル $\mu_\text{el}$ は波数空間において以下のように定式化されます。これは特定の異方性（直交するネットワーク形成など）を再現するための「2次元の現象論的・簡易的な弾性ペナルティ項」です。

$\mu_\text{el}(\mathbf{k}) \propto B_\text{elastic} \frac{k_x^2 k_y^2}{|\mathbf{k}|^4} \cdot \tilde{c}(\mathbf{k})$

3. Cahn-Hilliard方程式と計算アルゴリズム#

濃度の時間発展は、保存系を記述するCahn-Hilliard方程式に従います。フェーズフィールド法における標準的な定式化 [7] に倣い、4階の空間微分に伴う実空間での時間刻み制限を回避するため、高速フーリエ変換（FFT）を用いた擬スペクトル法（Pseudo-spectral method）と、時間積分に関する準陰解法（Semi-implicit method）を組み合わせた自作のPythonソルバーを構築し、高速かつ安定したシミュレーションを実現しました。

$\frac{\partial c}{\partial t} = \nabla \cdot \left( M \nabla \frac{\delta F}{\delta c} \right) \quad \xrightarrow{\text{Fourier}} \quad \frac{\partial \tilde{c}}{\partial t} = -M |\mathbf{k}|^2 \left( \tilde{\mu}_\text{chem} + \kappa |\mathbf{k}|^2 \tilde{c} + \tilde{\mu}_\text{el} \right)$

表1. フェーズフィールド解析に使用した主なパラメータ設定

パラメータ	設定値 / 物理的根拠
CNF初期体積分率 ( $c_0$ )	0.25 (25%) / 複合材料として実用的な充填率
熱力学的駆動力 ( $W$ )	HSP文献値より算出した $\chi$ パラメータに基づく
勾配エネルギー係数 ( $\kappa$ )	0.5 / 界面張力に基づく標準化パラメータ
弾性効果係数 ( $B_\text{elastic}$ )	0.0 (従来モデル) または 0.8 (Advance/PHASE連携による非均質弾性)
計算グリッド / ステップ	128 x 128 / 20,000ステップ (粗大化過程まで追跡)

なお、本事例のようなプロトタイプ計算には軽量な自作スクリプトが適していますが、より大規模な3次元計算や、熱拡散・複雑な結晶塑性といったマルチフィジックス連成が要求される高度な解析には、有限要素法ベースのマルチフィジックス汎用シミュレーションフレームワークであるMOOSE（Multiphysics Object-Oriented Simulation Environment）[8]を用いた拡張・置き換えも可能です。

解析結果と考察#

弾性効果が微細構造に与える劇的な変化#

図4は、フェーズフィールド法による相分離シミュレーション（20,000ステップ後）の濃度分布結果です。赤色がCNFリッチな析出相、青色がPLAマトリックス相を示しています。

(A) 弾性効果なし（従来の熱力学モデル： $B_\text{elastic} = 0.0$ ）
Advance/PHASEの力学情報を組み込まず、HSPから得た化学的自由エネルギーと界面張力のみで計算した結果です。界面エネルギーを最小化しようとするため、CNFドメインは綺麗な等方的な円形（球状）に成長し、PLAの海に孤立して浮かぶ「海島構造」が形成されます。オストワルド成長により円は徐々に粗大化しますが、ネットワークは形成されません。

(B) 弾性効果あり（Advance/PHASE連携モデル： $B_\text{elastic} = 0.8$ ）
第一原理計算から得た「非常に硬いCNF」と「柔らかいPLA」の弾性ミスマッチを組み込んだ結果です。球状に凝集しようとする界面張力よりも、系全体のひずみエネルギーを逃がそうとする力学的な駆動力が打ち勝ち、ドメインが縦横方向へ大きく引き延ばされています。25%という非対称なブレンド比率でありながら、ドメイン同士が互いに繋がり合い、繊維状・網目状の「パーコレーション・ネットワーク」を形成する「弾性相分離」の様子が再現されました。

相分離シミュレーションの比較

図4. 相分離シミュレーション結果の比較（CNF 25%, 20,000 steps）。（上）：弾性効果を無視した従来の等方的な海島構造。（下）：弾性率の効果を反映させた結果、ひずみ緩和のためにCNFが異方的に連なり、補強効果の源泉となるパーコレーション・ネットワークが自発的に形成されています。

まとめ#

本事例では、第一原理計算ソフトウェアAdvance/PHASEを用いたDFT計算において分散力と構造異方性を正確に考慮することで、セルロースI $\beta$ の妥当な弾性率（約19 GPa）を導き出しました。さらにその力学情報をメゾスケールのフェーズフィールド法へと連携させ、PLA/CNF複合材料における連続した補強ネットワークの自発的な形成プロセスを計算機上で実証しました。

本解析の詳細や、研究への適用に関するご相談はこちら

お問い合わせ

参考文献#

Y. Nishiyama, P. Langan, and H. Chanzy, "Crystal Structure and Hydrogen-Bonding System in Cellulose I $\beta$ from Synchrotron X-ray and Neutron Fiber Diffraction", J. Am. Chem. Soc. 124, 9074 (2002).
Y. Li, M. Lin, and J. W. Davenport, "Ab Initio Studies of Cellulose I: Crystal Structure, Intermolecular Forces, and Interactions with Water", J. Phys. Chem. C 115, 11533 (2011).
F. L. Dri, L. G. Hector Jr., R. J. Moon, and P. D. Zavattieri, "Anisotropy of the elastic properties of crystalline cellulose I $\beta$ from first principles density functional theory with Van der Waals interactions", Cellulose 20, 2703 (2013).
Y. Nishiyama, "Molecular interactions in nanocellulose assembly", Phil. Trans. R. Soc. A 376, 20170047 (2018).
C. M. Hansen, "Hansen Solubility Parameters: A User's Handbook", CRC Press, 2007.
A. G. Khachaturyan, "Theory of Structural Transformations in Solids", John Wiley & Sons, 1983.
L. Q. Chen, "Phase-Field Models for Microstructure Evolution", Annu. Rev. Mater. Res. 32, 113 (2002).
D. Gaston, C. Newman, G. Hansen, and D. Lebrun-Grandie, "MOOSE: A parallel computational framework for coupled systems of nonlinear equations", Nucl. Eng. Des. 239, 1768 (2009).