非晶態木質納米纖維素機械特性的建模1)

張秀梅 曹 軍 仇遜超 Mark A Tschopp Mark Horstemeyer Sheldon Shi

(東北林業大學，哈爾濱，150040) (美國密西西比州立大學) (美國北德克薩斯州大學)

納米材料的研究是現在的熱門研究方向，在現代科學中，人們從大量材料與結構的破壞事件中認識到材料的變形與破壞現象源自材料的微觀尺度及原子尺度［1］。將納米技術引進到木材科學中，極大地擴寬了木材科學的研究領域，使研究深度從細胞水平上升到分子水平。纖維素是世界上可再生、最為豐富的天然生物高分子聚合物，植物的主要結構成分是纖維素，其在木材中的所占比例為40%～50%［2］。木質納米纖維素是一種由1，4-β-D-葡萄糖甙鍵鏈接組成的無支鏈葡萄糖聚合體，其天然結構是由規則的晶體結構和無規則的非晶態結構構成，晶體結構決定了木材的彈性，非晶態結構決定了木材的靈活性與柔軟性，其機械特性由晶體結構與非晶態結構所占的比例來刻畫，具有很高的強度以及彈性模量等［3］。自19世紀人類發現木質纖維素中的晶體結構以來，多種不同種類木材的晶體結構已被確定。然而，由于木質纖維素中的非晶態結構的內部原子的排列呈現雜亂無章的分布狀態，使得對其機械特征進行相關的研究，存在一定的困難。

筆者利用分子動力學仿真的方法，對木質納米纖維素中的非晶態結構進行分子建模與拉伸變形仿真研究，通過所建模型仿真研究數據，求得納觀級結構中應力—應變曲線，將其同實驗研究方法中的數據相比較，從而預測出材料宏觀尺度各特性以及本構關系。對非晶態木質納米纖維素模型進行研究，對其結構及拉伸變形過程進行細致的分析，填補了國內外對木質納米纖維素的分子模擬研究的空白，具有重要的意義。

1 建立模型

纖維素是自然界中分布最廣、含量最多的一種多醣，由8 000～10 000個D-葡萄糖以β-1，4-糖苷鍵組成的無支鏈葡萄糖聚合體，相對分子質量約50 000～2 500 000，相當于300～15 000個葡萄糖基，分子式可寫作(C6H10O5)n。1年生或多年生植物，尤其是各種木材都含有大量的纖維素，主要的生物學功能是構成植物的支持組織。自然界中，植物體內約有50%的碳于纖維素的形式存在，一般木材中，纖維素占40%～50%，還有10%～30%的半纖維素和20%～30%的木質素［4-5］。

纖維素的材料力學性能主要依賴于其晶體結構和非晶態結構(無定型結構)的比值，纖維素晶體結構有較強的機械強度，可以作為復合材料的增強基物質;而非晶纖維素展現了重要的黏彈性。反而言之，各種超微結構的形成同時取決于單一分子鏈的構象特點及其在固態狀態下的交互作用。因此，對于木質納米纖維素的晶體和非晶態結構特點的研究有助于從分子水平理解木材的宏觀力學特性，把握木材力學性質的本質性起源。

到目前為止，尚沒有成熟的理論模型來描述納米尺度下木質纖維素的機械特性，對于木質纖維素納米結構的研究，主要是采用實驗的方法對納米纖維進行提取(酸解法、機械方法、生物菌培育)，再利用實驗手段進行分析［6-7］。由于提取到的納米纖維素的幾何尺寸微小，因此對實驗測定的設備和測試技術要求較高，僅僅通過實驗方法并不能完全了解納米材料的微觀結構、性能及內在機制。隨著計算機技術的進步，計算機模擬已廣泛應用于材料科學的研究中，分子動力學模擬可以彌補實驗的不足，獲得有關的原子運動細節，有效地澄清實驗現象，揭示內在機制。

采用美國Accelrys公司開發的Materials Studio軟件建立非晶態木質納米纖維素結構的分子模型，應用周期性邊界條件，并提取初始狀態的原子坐標、晶格長度、速度分布、環境狀態等數據，圖1為木質纖維素單體的分子結構。

圖1 木質纖維素單體的分子結構

圖2為微/納木纖絲非晶態結構的分子模型，其中圖2b為經過能量最小化后的分子模型。

木質納米纖維素結構是由一些長短不等的纖維素分子鏈聚集成的，在大分子鏈排列最致密的地方，分子鏈規則平行排列，定向良好，反映出一些晶體的特征，所以被稱為纖維素的結晶區，分子鏈與分子鏈間的結合力隨著分子鏈間距離的縮小而增大。當纖維素分子鏈排列的致密程度減小、分子鏈間形成較大的間隙時，分子鏈與分子鏈彼此之間的結合力下降，纖維素分子鏈間排列的平行度下降，此類纖維素大分子鏈排列特征被稱為纖維素非晶態結構(有時也稱作無定形區)。結晶區與非結晶區之間無明顯的絕對界限，而是在纖維素分子鏈長度方向上呈連續的排列結構。

圖2 構建的木質納米纖維素非晶態結構的分子模型

對納米木纖維非晶態結構進行建模時，首先建立一個單元的無規則晶胞，為了后期驗證不同鏈長對于建模仿真的機械特性的影響，所建模型中定義了3種不同長度的纖維素分子鏈，分別是由20、50、100個葡萄糖單體連接而成。定義纖維素分子鏈排列的致密程度為1.5 g/cm3，圖2a所示為100個葡萄糖單體構成的纖維素分子鏈，然后沿著x、y、z方向進行周期性拓展，從而構建出3×3×3的木質納米纖維素非晶態結構的模型。由于生成的無規則單胞，分子可能不是等價地分布在單胞中，從而造成了真空區，為了矯正此現象，就要進行能量最小化來優化此模型。圖2b為周期性擴展后并進行能量最小化后得到的分子模型，也是本研究用于拉伸變形仿真模擬的模型。

2 分子動力學

2.1 運動方程

分子動力學是一套分子模擬方法，依靠牛頓力學來模擬分子或原子的運動，通過求解各分子或原子的運動方程，得到其不同時刻的位置、動能和能量等，從而統計出材料的宏觀熱力學量和其他宏觀性質。

分子動力學的研究對象是由一定量的分子或原子組成的系綜。假設在包含有N個原子的系綜中，定義 r1，r2，…，rN為原子的位置矢量，U(r1，r2，…，rN)為描述N個原子間的作用力勢能函數，Fi=-▽riU(r1，r2，…，rN)為第 i個原子的受力，則用于描述第i個原子運動的牛頓方程為：

式中：mi為第i個原子的質量。

在實際計算中，根據加載情況給出相應的邊界條件與起始條件，系統在此時的狀態是強加的。它必須通過很多松弛步，逐漸達到真實的平衡位置，即滿足運動微分方程及初始與邊界條件，而達到總勢能最小的位置。為了得到較精確的結果，顯然要求從前一步到下一步的遞推關系給出的算法能進行高效率的計算，并且計算格式在數值上是穩定的。常用的有限差分算法有Verlet算法、Leap-frog算法、Velocity Verlet算法和 Beeman’s算法［8］。

2.2 初始條件和邊界條件

在進行分子動力學仿真前，需要對原子的空間位置ri、初始速度vi進行初始化。本研究利用Materials studio軟件構建木纖維非晶態納米結構的分子模型。然后通過Matlab軟件將該模型轉換為Lammps可讀取的文件，從而獲得了各個原子的空間位置與初始速度，實現了原子的初始化。

在分子動力學仿真中，通過選取一定的邊界條件，可以用有限的模擬區間得到具有宏觀尺寸的物理結構與材料的機械特性。本研究選取周期性邊界條件，引入周期性邊界條件的作用是：在原子的運動過程中，若有一個或幾個粒子跑出模型，則必有一個或幾個粒子從相反的界面回到模型中，從而保證該模擬系統的粒子數恒定;采取最近鏡像方法計算原子間的作用力，使得所模擬系統邊界處的原子的受力比較全面，從而消除了邊界效應。

2.3 反應性力場

本研究中，利用反應性力場(ReaxFF)進行拉伸變形仿真，模擬原子間的相互作用。反應性力場是加州理工 Adrivan Duin、William A.GoddardⅢ等人開發的力場，它是結合鮑林鍵長鍵能關系得到的半經驗勢，呈現出簡單、快速、放熱反應勢能的反應性力場，對描述鍵之間的斷裂與再結合很有利。大多數高能材料反應區域比較寬，反應時間比較長，由于過分消耗計算資源，模擬原子間的穩定爆破與凝聚幾乎是不可能的，為此產生了反應性力場經驗鍵級勢能，并且被延伸到扭轉效應，使它可用于依賴外部環境的原子間的相互作用，適合模擬木質納米纖維素非晶態結構的分子模型在拉伸變形過程中鍵之間的各種變化［9-11］。反應性力場將整個系統的能量劃分為幾個不同的能量部分，如下式所示：

式中：Ebond為共價相互作用能;Eunder和 Eover為對等/過等能;Eval為鍵角相互作用能;Epen為補償能;Etors為二面角扭轉能;Econj為鍵結合能;Evdwaals為范德華作用能;Ecoulomb為靜電作用能。可知，反應性力場由原子內相互作用和原子間相互作用兩大部分構成，即力場的勢能包括成鍵和非鍵相互作用，其中范德華作用能和靜電作用能屬于非鍵相互作用勢能。

2.4 系綜

分子動力學仿真需要在一定的系綜下進行，所謂系綜是指在一定的宏觀條件下，大量結構與性質完全相同的、處于各自獨立、各種運動狀態的系統的集合。系綜用于研究系統的微觀狀態與宏觀熱力學性質的對應規律。分子動力學仿真的不同系綜有：仿真系統中原子數N、體積V和能量E都保持不變的微正則系綜(NVE);仿真系統中原子數N、體積V和溫度T都保持不變的正則系綜(NVT);仿真系統中原子數N、溫度T和壓力P都保持不變的等溫等壓系綜(NPT);以及仿真系統中原子數N、壓力P和焓值H都保持不變的等焓等壓系綜(NHP)。

本研究對能量最小化后的納米木纖維非晶態結構的分子模型，在微正則系綜(NVE)條件下，實現熱平衡仿真。對熱平衡仿真后的模型，在等溫系綜條件下，進行拉伸速率恒定的拉伸變形仿真。

3 仿真模擬

3.1 仿真步長

分子動力學仿真中，最重要的工作為如何選取合適的積分步長，在節省時間的同時也保證計算的精確性。太長的步長會造成分子間的激烈碰撞，體系數據溢出;太短的步長會降低模擬過程搜索相空間的能力。本研究選取兩種積分步長，即Δt1=0.5 fs，Δt2=1 fs，觀察不同積分步長下，系統溫度、總能量以及動能、勢能的波動情況，從而確定積分步長的選取。兩種積分步長模型溫度的波動情況如圖3所示，其中圖3a選取的積分步長為Δt1=0.5 fs，圖3b選取的積分步長為Δt2=1 fs。

圖3 兩種積分步長系統溫度的波動情況

熱平衡的仿真在基于反應性力場的微正則系綜(NVE)條件下進行，由圖3可知，在積分步長選取為Δt1=0.5fs的條件下，模型未經過震蕩達到了熱力學平衡態;而在積分步長選取為Δt2=1fs時，模型經過多次震蕩才趨于平衡。模型的能量(總能量、動能及勢能)的波動情況如下圖4所示，其中圖4a選取的積分步長為Δt1=0.5fs，圖4b選取的積分步長為 Δt2=1fs。

圖4 兩種積分步長系統總能量、動能及勢能的波動情況

能量變化的仿真在基于反應性力場的微正則系綜(NVE)條件下進行，由圖4可知，在積分步長選取為Δt1=0.5fs的條件下，模型的能量(總能量、動能及勢能)在沒有經過明顯震蕩的情況下較快速的達到平衡態;而在積分步長選取為Δt2=1fs的條件下，系統的能量(總能量、動能及勢能)有明顯的震蕩，并且在原子間距離過近時，易發生原子回彈現象。因此，在今后采用基于反應性力場的分子動力學仿真研究中，應選取0.5fs的積分步長。

3.2 物理與機械特性

通過計算模型的密度、彈性模量、剪切模量以及泊松比，來反映構建的木質納米纖維素非晶態結構的分子預測模型的物理與機械特性，將計算得到的數據與在文獻中采用實驗研究方法獲得的結果進行比較，從而判斷該預測模型的精準性，比較結果如表1所示。可以觀察到，對構建的木質納米纖維素非晶態結構分子模型，采用基于反應性力場的分子動力學進行原子間相互作用的仿真，從而計算得到木纖維非晶態結構納觀尺度的物理與機械特性，與文獻［12］中給出的相關特性值十分相似。因此，利用構建的理想模型，采用基于反應性力場的分子動力學模擬原子間的作用，對納米木纖維非晶態結構進行進一步的相關研究是可行的，并且研究結果是完全具有可信度的。計算理想模型得到的相關特性值與文獻中的相關特性值存在細微差異的原因是：構建的理想模型是均質化的，而采用實驗研究方法獲得的微/納木纖絲，其纖維素分子鏈、晶體與非晶態結構的分布均是非均質的;并且采用實驗研究方法獲得的微/納木纖絲可能會含有一定的濕度及雜質，從而造成了細微的差異。

表1 微/納木纖絲非晶態結構的物理與機械特性

3.3 拉伸仿真

在對納米木纖維非晶態結構的分子模型進行基于反應性力場分子動力學拉伸變形仿真的過程中，只對x軸方向施加單軸應力，而保證由y軸、z軸構成的二維平面的應力為0。拉伸仿真是在模型經過能量最小化及熱平衡仿真后，在等溫系綜下進行的，分子動力學計算使用開源代碼程序LAMMPS。拉伸仿真前后的原子結構如下圖5所示。

圖5 拉伸仿真前后的分子結構

拉伸變形仿真的過程中，選取的溫度為T=100 K，低于玻璃化轉變溫度;選取的拉伸速率為1010s-1，仿真中的原子數為28 404個。圖5a為拉伸仿真前的原子結構，圖5b為拉伸仿真后的原子結構，可以看出，沿拉伸x方向出現拉長現象，而y和z方向則出現長度縮減，正是由于采用ReaxFF力場進行仿真研究，結果符合實際拉伸過程中的變化。

對木質納米纖維素非晶態結構分子模型進行基于反應性力場分子動力學拉伸變形仿真后，繪制其應力—應變曲線，如圖6所示。

圖6 應力—應變曲線

圖6中的黑色曲線為原子受到x軸方向施加的單軸應力的瞬間分布情況的平均應力響應曲線。可以觀察到，納米纖維素非晶態結構分子模型首先經歷彈性變形，接著在0.5 GPa左右出現應力平臺，平臺期的應變長達35%，在平臺期結束后繼續拉伸時，應力重新上升到0.6 GPa左右，并發生彈性屈服。拉伸仿真后繪制的應力—應變曲線與在文獻中采用實驗研究方法獲得的應力—應變曲線具有相同的變化規律。

4 結論

首先對木質納米纖維素非晶態結構進行周期邊界條件建模，然后對該模型進行能量最小化和熱平衡處理，利用基于反應性力場對模型進行拉伸變形仿真，最終繪制出預測模型的應力—應變曲線。通過對構建的納米木纖維非晶態結構的熱平衡及能量的仿真，確定了在今后的采用基于反應性力場的分子動力學仿真研究中，應選取0.5fs的積分步長;計算模型的物理與機械特性，并與采用實驗研究方法獲得的相應特性進行比較，表明了采用本研究方法對木質納米纖維素非晶態結構進行進一步的相關研究是可行的，并且研究結果完全具有可信度;對模型的應力—應變曲線進行分析，表明了采用本研究方法與采用實驗研究方法獲得的應力—應變曲線具有相同的變化規律。本研究屬于交叉和前沿的研究，填補了國內外對木質納米纖維素非晶態結構的分子模擬研究的空白，對提高木材以及木質復合材料從微觀到宏觀的力學性能認識、完善木材納米微觀力學具有重要促進作用。預期研究成果可用于各種工程產品，為終端客戶提供科學試驗數據以及可視化的機械特性效果，推動納米技術在木材科學領域的研究，在新的層次上驗證可持續發展的理論。

［1］黃克智，吳堅，黃永剛.納米力學的興起［J］.機械強度，2005，27(4)：403-407.

［2］邱堅，李堅.納米科技走進木材科學［J］.國際木業，2003(1)：10-11.

［3］Lima M，Borsali R.Rodlike cellulose microcrystals：Structure，properties，and applications［J］.Macromolecular Rapid Communications，2004，25(7)：771-787.

［4］Kamide K.Cellulose and cellulose derivatives：molecular characterization and its applications［M］.The Netherlands：Elsevier，2005.

［5］O'Sullivan A C.Cellulose：the structure slowly unravels［J］.Cellulose，1997，4(3)：173-207.

［6］Yoshiharu N，Junji S，Henri C，et al.Crystal structure and hydrogen bonding system in cellulose Iβfrom synchrotron X-ray and neutron fiber diffraction［J］.Journal of the American Chemical Society，2003，125(47)：14300-14306.

［7］Samira E H，Yoshiharu N，Jean L P，et al.The shape and size distribution of crystalline nanoparticles prepared by acid Hydrolysis of native cellulose［J］.Biomacromolecules，2008，9(1)：57-65.

［8］Frenkel D，Smit B，Ratner M A.Understanding molecular simulation：from algorithms to applications［M］.Orlando，USA：Academic Press，1996.

［9］Hasan M A，Sagar A P，Adri C D，et al.Reactive molecular dynamics：numerical methods and algorithmic techniques［J］.International Journal of Mechanical Sciences，2010，22：1-29.

［10］Kimberly C，Adri CD，William A G.ReaxFFreactive force field for molecular dynamics simulations of hydrocarbon oxidation［J］.The Journal of Physical Chemistry A，2008，112(5)：1040-1053.

［11］Adri CD，Siddharth D，Francoins L，et al.ReaxFF：A reactive force field for hydrocarbons［J］.The Journal of Physical Chemistry A，2001，105(41)：9396-9409.

［12］Chen W，Lickfield GC，Yang CQ.Molecular modeling of cellulose in amorphous state.Part I：model building and plastic deformation study［J］.Polymer，2004，45(3)：1063-1071.