曾凡林，孫 毅，劉一志

（哈爾濱工業(yè)大學航天科學與力學系，150001 哈爾濱）

PVDF單軸拉伸的第一性原理和分子動力學模擬

曾凡林，孫 毅，劉一志

（哈爾濱工業(yè)大學航天科學與力學系，150001 哈爾濱）

為得到聚偏二氟乙烯(PVDF)的大變形拉伸特性，應(yīng)用第一性原理方法模擬了β相PVDF分子鏈的在單軸拉伸下的能量及構(gòu)型變化，得到了β相分子鏈拉伸過程中發(fā)生斷裂時的臨界鍵長、臨界斷裂能量、臨界拉力以及分子鏈彈性系數(shù)的變化規(guī)律.接著用一種新穎的方法構(gòu)造了無定形PVDF的立方元胞模型并用于單軸拉伸的分子動力學模擬.模擬中分別考慮鏈不被拉斷和鏈允許被拉斷兩種情況.第1種情況下，拉伸應(yīng)變率分別為1010，109和108s－1，得到了應(yīng)力應(yīng)變關(guān)系和構(gòu)型變化;第2種情況下，拉伸應(yīng)變率分別為109s－1，5×108s－1和準靜態(tài)，得到了應(yīng)力應(yīng)變關(guān)系和斷裂前后的構(gòu)型.同時分析了一些關(guān)鍵的形變信息和影響模擬結(jié)果的一些關(guān)鍵因素.模擬結(jié)果表明，不僅對于理解PVDF的拉伸變形機制具有較重要意義，而且對于其他無定形材料的拉伸模擬具有參考價值.

PVDF;單軸拉伸;第一性原理;分子動力學;斷裂

作為一種被廣泛應(yīng)用的壓電聚合物，聚偏二氟乙烯(Poly(vinylidene difluoride)，PVDF)還應(yīng)用于對熱穩(wěn)定性、耐化學腐蝕性和耐磨性有特殊要求的領(lǐng)域.因為有與其他聚合物相比更卓越的力學性能以及突出的壓電、熱電性能，PVDF一直以來都引起廣大研究者的注意［1－2］.PVDF 的壓電性能主要源自其分子結(jié)構(gòu)中的β相和γ相.相比于其他的壓電材料例如 PZT，雖然 PVDF的力——電耦合效應(yīng)低于PZT，但是它的可塑性要高出10倍左右［3－4］;此外作為傳感材料更重要的一個特征參數(shù)PVDF的g要比PZT高20倍左右，而且它能被制成厚度＜10 μm的薄膜，振動質(zhì)量極小，具有比PZT大得多的阻尼值，因此具有產(chǎn)生極短波長所需的動態(tài)特性.正因為如此，PVDF正被嘗試應(yīng)用于各種航天智能系統(tǒng)結(jié)構(gòu)中，例如可充氣式展開天線等，有望在各類航天器智能結(jié)構(gòu)中得到應(yīng)用［5－6］.

單軸拉伸是制備β相PVDF的傳統(tǒng)且有效的方法.拉伸過程中，應(yīng)變高達400%，使得PVDF的分子鏈由蜷曲的無定形態(tài)轉(zhuǎn)變?yōu)橹辨湋B(tài).然后在強電場的作用下利用CH2和CF2基團的極性，使得分子鏈中原先呈交錯狀態(tài)的CH2和CF2的相對位置發(fā)生改變，所有的H原子位于一側(cè)，而所有的F原子則轉(zhuǎn)向另一側(cè)，使材料產(chǎn)生宏觀的整體極性，從而具有壓電性能.實驗結(jié)果表明，如果分子鏈之間交聯(lián)度較小的話，較小的拉應(yīng)變即會使鏈和鏈發(fā)生分離從而產(chǎn)生鏈間斷裂.因此在制備無定形的PVDF時，盡可能多地在分子鏈間形成相互交聯(lián)的網(wǎng)狀結(jié)構(gòu)是拉伸法成功的關(guān)鍵.但這樣會產(chǎn)生另外一個問題，即是相互交聯(lián)的分子鏈中不可避免地會存在死纏接.一旦這種情況出現(xiàn)，即使較小的拉伸應(yīng)變也會使得分子鏈斷裂，從而在斷裂處形成空位或銀紋，極化時非常容易被擊穿，使得材料失效.因此研究無定形的PVDF分子鏈間的死纏接以及拉伸過程中的分子鏈斷裂情況，對于PVDF的制備和應(yīng)用是非常有意義的.

因為傳統(tǒng)的實驗手段很難觀察到PVDF內(nèi)部的分子鏈結(jié)構(gòu)，更難以表征拉伸過程中鏈的解纏接和斷裂過程，因此本文采用計算機模擬的方法.分子鏈的斷裂需要了解單個鏈在拉伸過程的構(gòu)型信息和能量信息從而得到臨界斷裂值，本文采用第一性原理來模擬單鏈的拉伸來進行;鏈間的解纏結(jié)或斷裂過程則通過大規(guī)模的模擬來進行，本文采用分子動力學模擬方法.第一性原理對單鏈的模擬結(jié)果直接應(yīng)用于大規(guī)模的分子動力學模擬中，兩者結(jié)合成功地對PVDF的大變形拉伸過程進行了模擬.該方法對于其他材料的大變形分子動力學模擬同樣具有參考價值.

1 PVDF單鏈拉伸的第一性原理模擬

1.1 模擬模型及方法

本文依據(jù)文獻［7］構(gòu)造了β相PVDF的分子鏈，如圖1所示.模型采用周期性邊界條件以利于實施鏈長方向的軸向加載.為了減少周期性邊界所產(chǎn)生的鏡像分子鏈對當前分子鏈的影響，本文將模擬元胞的垂直于鏈長方向(x和y方向)的尺寸設(shè)置成遠遠大于截斷半徑的1.5 nm.為了使得變形只發(fā)生在鏈長方向，將元胞垂直于鏈長方向(x和y方向)的尺寸固定.為了減少計算量，模型中只取到兩個分子鏈的重復單元，整個模型含有12個原子，其優(yōu)化后的x、y、z方向的尺寸分別為1.5 nm(固定)、1.5 nm(固定)和 0.509 9 nm.

圖1 β相PVDF分子鏈的構(gòu)型及模擬模型

鏈長方向的單軸拉伸過程的第一性原理模擬采用Accelrys公司的 Materials Studio軟件中的Castep模塊來進行.通過改變元胞的z向尺寸來實現(xiàn)加載過程，每步加載應(yīng)變?yōu)?%，加載步數(shù)不限，直至分子鏈發(fā)生斷裂.加載過程中x和y方向的尺寸保持不變，具體的模型及加載信息參考文獻［8］.本文在計算中采用了Vanderbilt超軟贗勢和 LDA 的 Ceperley-Alder形式的函數(shù)［9］，布里淵區(qū)采用2×2×5的 k點取樣(其收斂性設(shè)為1.0 e－8eV/原子)來進行數(shù)值積分.元胞的能量計算中收斂標準設(shè)置為1.0 eV/原子.整個模擬過程參考文獻［10］方法.

考慮到模型只包含單個分子鏈并且分子鏈之間的相互影響忽略不計，此外分子鏈只在鏈長方向受到約束而x和y方向形變自由，因此拉伸過程中能量的變化可看作全部由鏈長方向的變形所貢獻.這樣分子鏈所受到的拉力F可表示為F=－ (?ΔE/?εzz)/d0.其中:d0為元胞的 z向原始尺寸;εzz為當前的拉伸應(yīng)變;ΔE為當前拉伸加載后相對于初始結(jié)構(gòu)的能量變化.同樣，C-C鍵的彈性系數(shù)k可表示為:k= － (?ΔF/?εzz)/d0.其中F 為計算得到的分子鏈所受到的拉力.

1.2 模擬結(jié)果

PVDF單分子鏈在鏈長方向的單軸拉伸下，能量變化如圖2所示.曲線顯示系統(tǒng)能量隨著應(yīng)變的增加而不斷升高，當應(yīng)變εzz=0.3時，能量發(fā)生突降，預示著分子鏈可能已經(jīng)被拉斷，這一點可被隨后的構(gòu)型觀測結(jié)果證實.考慮到本文的模型選取較小，能量的變化主要來自于C-C之間的共價鍵作用，由此得到β相PVDF分子鏈上C-C鍵的臨界斷裂能約為1.14 eV.

圖2 單軸拉伸載荷下PVDF單分子鏈的能量變化曲線

當εzz=0.29(第29步)和εzz=0.3(第30步)時，分子鏈構(gòu)型如圖3(a)、(b)所示.顯然，拉伸到第30步時(εzz=0.3)分子結(jié)構(gòu)已經(jīng)失穩(wěn)，相鄰兩個C-C鍵顯示巨大差異(鍵長分別為0.292 6、0.131 3 nm)，鍵長較大的那個 C-C 鍵(0.292 6 nm)已經(jīng)被拉斷.而在此前的拉伸步(εzz=0.29)，兩個相鄰C-C鍵的鍵長接近(分別為0.182 5 、0.183 2 nm)，表明此時結(jié)構(gòu)穩(wěn)定，分子鏈并未被拉斷.構(gòu)型分析的結(jié)果與前面能量變化的結(jié)果一致，當應(yīng)變εzz達到0.3時，β相PVDF單分子鏈將會被拉斷，其臨界斷裂鍵長可認為是0.183 2 nm.

圖3 PVDF分子鏈在應(yīng)變?yōu)?.29和0.30時的構(gòu)型

拉伸過程中PVDF單分子鏈的拉力變化曲線和C-C鍵的彈性系數(shù)k的變化分別如圖4(a)、(b)所示.拉力變化曲線顯示，當能量發(fā)生突變時(εzz=0.3，對應(yīng)分子鏈被拉斷)拉力也發(fā)生了突變.而在εzz=0.28時拉力達最大值，為6.75 nN，顯示PVDF單鏈在斷裂之前有一個非常微小的屈服過程，屈服拉力約為6.75 nN.圖4(b)的曲線顯示只有當應(yīng)變非常小時(εzz＜0.03，對應(yīng)著彈性理論里面的小變形)，彈性系數(shù)k才可近似認為是常數(shù)，而當應(yīng)變增大時，k迅速降低，在斷裂之前趨于零.普通的勢函數(shù)描述C-C之間的彈性能時并未考慮這一因素，而是處理成常數(shù)，由此可見普通的勢函數(shù)(分子力場)只適用于彈性范圍內(nèi)的小變形模擬，對于大變形的模擬則必須考慮這一因素.

圖4 拉伸過程中PVDF單分子鏈的拉力變化和C-C鍵的彈性系數(shù)變化曲線

1.3 討論

由于第一性原理模擬的計算量大，本文的模擬所采用的模型較小，模擬所選取的精度也不是非常高.不過考慮到研究β相分子鏈拉伸過程中的斷裂性能并用于分子動力學的模擬中，這樣的計算精度是可以接受的.對單鏈的模擬結(jié)果表明，β相的PVDF分子鏈被拉伸時最終會發(fā)生鏈的斷裂，而在這個過程中，系統(tǒng)能量和C-C共價鍵彈性系數(shù)k的變化都是非線性的.因為高分子材料的分子動力學模擬所用的分子力場如COMPASS、PCFF、MM3、UFF等都只能模擬小的彈性變形，無法描述大變形甚至斷裂過程，本文所得到的這些拉伸結(jié)果及影響參數(shù)都是在利用分子動力學模擬來研究塊體材料大變形拉伸性質(zhì)時需要考慮并注意的問題.

2 無定形PVDF的元胞模型的構(gòu)建

在以往研究中給出了高分子材料的單分子鏈模型和包含多個無定形分子鏈的元胞模型的構(gòu)建方法［11－12］.為了使構(gòu)建的模型更加與實際的分子鏈構(gòu)型相符，本文采用了一種較為新穎的方法來進行.該方法主要采用Materials Studio軟件中的Perl script腳本語言完成.其基本步驟及實現(xiàn)過程為:1)將80個聚合度為200的PVDF分子鏈按照1.60 g·cm－3的密度裝入長方體元胞中建立具有周期性邊界條件的初始長方體模型(其長寬高a、b、c分別為 a=9.335 nm、b=9.335 5 nm、c=12 nm).在這一過程中，分別將每個分子鏈端部的C原子設(shè)置成“活性”原子，意味著在隨后的過程中，這些較短的鏈之間有機會相互連接而形成長鏈.2)在NVT系統(tǒng)下對第1步建立的初始模型進行分子動力學模擬.在模擬過程中，每個“活性”C原子不斷地在其截斷半徑內(nèi)(1 nm)尋找其他分子鏈上的“活性”原子，一旦發(fā)現(xiàn)有滿足條件的原子，兩個“活性”原子發(fā)生相互連接，同時刪掉“活性”標記，變?yōu)椤胺€(wěn)定”原子，這一過程貫穿整個動力學模擬過程，直到該“活性”原子完成連接或自始至終也無法找到可連接的原子.為了避免同一鏈上的兩個“活性”原子相互連接而形成環(huán)鏈，此處又對每個原子加了一個額外標記用以記錄分子鏈信息，那些只有來自不同的分子鏈上的“活性”原子才能發(fā)生反應(yīng)，并且一旦發(fā)生反應(yīng)后，新形成的鏈上的其他“活性”原子重新進行標記用以表明它們位于同一分子鏈上.3)檢查生成的含有不同長度的分子鏈的模型的完整性，對不飽和C原子進行加H處理并對模型進行徹底的NVT模擬進行能量優(yōu)化和弛豫后，得到的模型作為單軸拉伸的分子動力學模擬模型.本文利用這種方法最終建立的模型包含的分子鏈數(shù)目為30(一共是96 060個原子)，最大的聚合度為2 200(意味著該分子鏈由11個初始分子鏈連接而成)，最小的200(該分子鏈未和其他分子鏈相互連接)，更多的分子鏈聚合度為3～8(整體符合正態(tài)分布).這樣的模型與實際情況更加吻合，因為在實際材料中，PVDF的鏈長是不盡相同的，從聚合度為數(shù)萬到數(shù)千的都可能存在，此處建立的模型在一定程度上可以描述分子鏈間的這種差異.圖5顯示了初始模型和經(jīng)過上述方法處理過的模型，可明顯見到分子鏈數(shù)目減少且變得更長.

圖5 “活性”原子發(fā)生反應(yīng)前和發(fā)生反應(yīng)后的無定形PVDF的元胞模型

3 無定形PVDF單軸拉伸的分子動力學模擬

利用無定形PVDF的元胞模型構(gòu)建的模擬模型，采用LAMMPS軟件和PCFF分子力場對無定形PVDF的單軸拉伸進行了分子動力學模擬，模擬的溫度為300 K，最大拉伸應(yīng)變設(shè)為1，常應(yīng)變率模擬.模擬中模型沿著z軸(c方向)被拉伸，兩端各占模型長度10%的部分被固定，3個方向都采取周期性邊界條件.模擬分兩種情況:第1種在拉伸過程不考慮分子鏈斷裂(傳統(tǒng)模擬)，模擬采用的應(yīng)變率分別為1010、109和108s－1;第2種情況考慮了鏈的斷裂情況，其應(yīng)變率分別為109s－1，5×108s－1和準靜態(tài)(模擬應(yīng)變率為零).后一種情況下，依據(jù)PVDF單鏈拉伸的第一性原理模擬結(jié)果，設(shè)置C-C鍵的臨界斷裂鍵長為0.184 nm.實際的模擬中，依據(jù)模擬結(jié)果，實際得到的最大應(yīng)變從0.22～1.00不等.因為較低的應(yīng)變率下模擬需要耗費非常長的時間，本文只進行了最低應(yīng)變率為108s－1的動態(tài)模擬.

圖6顯示了不同應(yīng)變率下不同情況的模擬中所得到的應(yīng)力－應(yīng)變曲線.當不考慮C-C鍵的斷裂時，應(yīng)變?yōu)?1010，109，108s－1的應(yīng)力 － 應(yīng)變曲線非常相似，應(yīng)力能夠隨著應(yīng)變的增加而達到非常高的值，而這與實際結(jié)果是不相吻合的.事實上在如此高的應(yīng)變率下，PVDF很容易被拉斷.只有在準靜態(tài)拉伸下才可能達到非常高的應(yīng)變.而當考慮C-C鍵的斷裂因素時，應(yīng)變率為109s－1的應(yīng)力－應(yīng)變曲線顯示當應(yīng)變達到較低的0.21時，模型即被拉斷，應(yīng)力急劇下降，這表明模型實際是在沖擊作用下發(fā)生的斷裂，109s－1的應(yīng)變率實際對應(yīng)著很高的沖擊速度;當應(yīng)變率降為5×108s－1時，實際結(jié)果相差不大，發(fā)生斷裂時的應(yīng)變僅僅有少許的增加.這些結(jié)果表明動態(tài)模擬PVDF的拉伸過程，因為計算速度的限制，所選擇的應(yīng)變率太高而難以模擬實際拉伸過程.為此本文采取了準靜態(tài)模擬，并考慮鍵的斷裂因素.從準靜態(tài)所對應(yīng)的應(yīng)力－應(yīng)變曲線上，可以發(fā)現(xiàn)當應(yīng)變較小時，其曲線形狀與其他情況下得到的結(jié)果是一致的，但是當應(yīng)變在0.15附近時，應(yīng)力有一個小的突降，但整個曲線走勢變化不大，表明在這個應(yīng)變下模型內(nèi)部可能發(fā)生了某個分子鏈的斷裂，但是直到應(yīng)變達到0.85時，從應(yīng)力的變化上看，整個模型都未發(fā)生斷裂.這樣的結(jié)果與實際情況是相吻合的.由此可見，考慮鍵斷裂因素的準靜態(tài)模擬是模擬大變形的有效途徑.

進一步地，當不考慮鍵的斷裂時應(yīng)變率為109s－1下和考慮鍵斷裂時準靜態(tài)下的部分模型構(gòu)型如圖7所示.可明顯看出，當應(yīng)變率為109s－1，ε=0.34時，模型中出現(xiàn)明顯的微孔洞，而當ε=0.85時，模型中微孔洞發(fā)展成了比較大的空洞，而且銀紋也已連接成大的微孔洞，模型本質(zhì)上已經(jīng)破壞，但因為鏈不斷裂，所以整個模型被拉成了兩部分，中間由幾個分子鏈保持連接.而當考慮鍵斷裂因素，在準靜態(tài)加載條件下，即使當ε=0.85，模型中也僅見少量微孔洞，而未見整體斷裂跡象，但是通過對每個分子鏈的構(gòu)型進行分析，確實發(fā)現(xiàn)了有3個分子鏈已被拉斷，這也證明了上述應(yīng)力－應(yīng)變分析的結(jié)果.由此可見在大變形的拉伸過程中，PVDF分子鏈確實會發(fā)生死纏接現(xiàn)象而被拉斷，但整體體現(xiàn)得并不明顯，即使對每個分子鏈都采取了較低的斷裂準則(鍵長為0.184 nm)，整體模型在較大應(yīng)變下仍未發(fā)現(xiàn)明顯的損傷破壞.

圖6 無定形PVDF在考慮鍵斷裂和不考慮鍵斷裂時在不同應(yīng)變率下的應(yīng)力－應(yīng)變響應(yīng)

圖7 PVDF的構(gòu)型變化

4 結(jié)論

1)β相PVDF分子鏈在鏈長方向被持續(xù)拉伸會發(fā)生斷裂，斷裂時有一個微小的屈服過程，分子鏈的彈性系數(shù)只有當應(yīng)變很小時才可以近似看作是常數(shù).

2)分子鏈被拉斷時存在臨界拉力、臨界鍵長及臨界斷裂能，這些參數(shù)可作為大變形模擬時分子鏈的斷裂準則.

3)構(gòu)造大規(guī)模無定形材料模擬模型的方法對于構(gòu)造與實際構(gòu)型相吻合的無定形材料的模擬模型具有很好的實用性.

4)應(yīng)變率是在拉伸分子動力學模擬時必須考慮的重要因素，并且只有當考慮鍵的斷裂因素時才能夠較真實地模擬材料的大變形行為，PVDF在大變形的拉伸下，存在少量分子鏈發(fā)生死纏接而被拉斷的現(xiàn)象，但整體體現(xiàn)得并不明顯.

［1］YEE W A，KOTAKI M，LIU Ye，et al.Morphology，polymorphism behavior and molecular orientation of electrospun poly(vinylidene fluoride)fibers［J］.Polymer，2007，48(2):512－521.

［2］WANG Mian，SHI Jiahua，PRAMODA K P，et al.Microstructure，crystallization and dynamic mechanical behaviour of poly(vinylidene fluoride)composites containing poly(methyl methacrylate)-grafted multiwalled carbon nanotubes［J］.Nanotechnology，2007，18(23)，235701.

［3］黃維恒，聞建勛.鐵電壓電高分子材料［M］//聞建勛，王賢姍.高技術(shù)有機高分子材料進展.北京:化學工業(yè)出版社，1994:526－553.

［4］陳貽瑞，王健.基礎(chǔ)材料與新材料［M］.天津:天津大學出版社，1994.

［5］TUZZOLINO A J，ECONOMOU T E，CLARK B C，et al.Dust measurements in the coma of comet 81P/Wild 2 by the dust flux monitor instrument［J］.Science，2004，304(5678):1776－1780.

［6］McDONNELL J A M，BURCHELL M J，GREEN S F，et al.The stardust dust flux monitor［J］.Advances in Space Research，2000，25(2):335－338.

［7］KARASAWA N，GODDARD W A.Force-fields，structures，and properties of poly(vinylidene fluoride)crystals［J］.Macromolecules，1992，25(26):7268 －7281.

［8］曾凡林，孫毅.PVDF單鏈拉伸的第一性原理模擬［C］//第七屆海峽兩岸工程力學研討會論文摘要集.哈爾濱:中國力學學會，2011:111－116.

［9］CEPERLEY D M，ALDER B J.Ground-state of the electron-gas by a stochastic method［J］.Physical Review Letters，1980，45(7):566 －569.

［10］FUJINAMI Akinori，OGATA Shigenobu，SHIBUTANI Yoji.Ab initio study of the tensile behavior of single polyimidemolecularchain ［J］. Polymer， 2004，45(26):9023－9028.

［11］曾凡林，孫毅，周玉.有機材料的分子模擬模型(I):單分子鏈的構(gòu)建［J］.哈爾濱工業(yè)大學學報，2009，41(10):95－99.

［12］曾凡林，孫毅，周玉.有機材料的分子模擬模型(II):模擬元胞的構(gòu)建［J］.哈爾濱工業(yè)大學學報，2010，42(2):281－287.

Ab initio and molecular dynamics study of uniaxial tensile behaviors of PVDF

ZENG Fan-lin，SUN Yi，LIU Yi-zhi

(Dept.of Astronautic Science and Mechanics，Harbin Institute of Technology，150001 Harbin，China)

To obtain the behaviors of Poly(vinylidene difluoride)(PVDF)under large tensile deformations，this paper investigates the tensile behavior of a single molecular chain of rod-like PVDF with β phase structure at zero temperature using a Density Functional Theory(DFT)calculation.We obtain the critical bond length，critical fracture energy，critical tensile force and the elastic constant when the β phase PVDF chain is stretched to break.Then a novel method is proposed to build the cell model of the amorphous polymer，which is used for the molecular dynamics(MD)simulations of uniaxial tensile of PVDF.During the simulations，two cases that the chains will never break and are permitted to break are considered.In the first case，the simulations with three different strain rates 1010s－1，109s－1，and 108s－1were carried out and the strain rate dependence of the stress-strain relationship，the structure deformation were analyzed as well.In the second case，the simulations with a strain rate of 109s－1，5 ×108s－1and a quasistatic strain loading were performed and the stressstrain relationship，the cell structures before and after the breaking were researched.Some key factors to influence the simulation results and some important deformation information were analyzed.This is helpful to elucidate the inherent deformation mechanisms during the tensile procedure.

PVDF;uniaxial tensile;DFT;molecular dynamics;fracture

O631

A

0367－6234(2012)09－0046－05

2011－08－26.

國家自然科學基金資助項目(11102053);哈爾濱工業(yè)大學科研創(chuàng)新基金資助項目(HIT.NSRIF.2010070);哈爾濱市科技創(chuàng)新人才研究專項資金資助項目(2012RFQXG001).

曾凡林(1978—)，男，副教授，碩士生導師;

孫 毅(1961—)，男，教授，博士生導師.

曾凡林，zengfanlin@hit.edu.cn.

(編輯 張 紅)