李君，孫毅，曾凡林

(哈爾濱工業(yè)大學(xué)航天科學(xué)與力學(xué)系，黑龍江哈爾濱150001)

多面體低聚倍半硅氧烷(polyhedral oligomeric silsesquioxane，POSS)是一類結(jié)構(gòu)簡式為(RSiO 1.5)n(n≥4且為偶數(shù)，R為有機(jī)取代基團(tuán))的化合物.POSS單體可以與有機(jī)聚合物進(jìn)行化學(xué)連接而生成POSS有機(jī)-無機(jī)納米雜化材料［1］，從而提高聚合物的耐磨性、耐高溫性、硬度、強(qiáng)度等一系列的力學(xué)性能［2-4］.

POSS對基體材料的性能產(chǎn)生影響的機(jī)理缺乏理論.文獻(xiàn)［5-6］通過分子動力學(xué)模擬方法，計算了材料的多種參數(shù)，結(jié)論認(rèn)為POSS單體由于體積巨大，摻雜在材料中阻礙了材料分子的運(yùn)動，從而使材料的性能得到提高.這可以用來解釋增強(qiáng)的原因，但是在聚合物中加入POSS后，有些彈性模量確實得到了增強(qiáng)［7］，但是有些彈性模量反而降低了［7-8］.為了進(jìn)一步分析POSS對材料改性的原因，本文采用分子動力學(xué)方法模擬了聚乙烯(PE)及含有不同POSS比例的雜化聚乙烯(POSS-PE)在不同溫度下的行為，計算得到了4種材料的玻璃態(tài)轉(zhuǎn)變溫度(Tg)及彈性模量.最后通過各部分能量的具體變化來尋求POSS的加入對基體材料力學(xué)性能影響的原因.

1 模擬實驗

1.1 模型

由于PE的結(jié)構(gòu)簡單，模擬起來相對容易，因此選用 PE作為基體材料.而 POSS單體采用T8［CH2CH3］7(CH2)8CH=CH2［9］，其結(jié)構(gòu)如圖1 所示.

圖1 POSS單體結(jié)構(gòu)圖(H原子沒有標(biāo)出)Fig.1 Illustration of the POSS unit(H atoms are not shown)

POSS單體通過一系列化學(xué)反應(yīng)與PE主鏈結(jié)合生成的POSS-PE結(jié)構(gòu)如圖2所示.

圖2 POSS-PE主鏈結(jié)構(gòu)圖(H原子沒有標(biāo)出)Fig.2 Illustration of the POSS-PE chain(H atoms are not shown)

受計算條件所限，在實際模擬中，PE的主鏈由500個乙烯單體聚合而成，而整個PE模型包括6條主鏈，采用周期性邊界條件［10］打包在一個元胞內(nèi).對于POSS-PE，構(gòu)建了3種含有不同POSS比例的模型，即分別在PE主鏈上加入2個，4個和6個POSS單體，從而生成3種POSS-PE主鏈.然后采用周期性邊界條件分別將每種POSS-PE主鏈打包在各自的元胞中，最后生成所需的3種POSS-PE模型.為簡便起見，稱其為PEII、PEIV和PEVI.4個模型的詳細(xì)信息如表1所示.

表1 4種聚合物的參數(shù)Table 1 Details of the four polymers

1.2 勢函數(shù)

勢函數(shù)作為分子模擬的核心，在分子模擬過程中占有非常重要的地位.在本文進(jìn)行的所有模擬中，都采用CVFF力場［11］作為勢函數(shù)，因為其形式簡單，而且被廣泛地應(yīng)用在聚合物的模擬上［12-15］.

CVFF力場的表達(dá)式為

式中:Kb、α是與鍵類型相關(guān)的參數(shù)，b0是鍵的平衡長度;Kθ是與鍵角類型相關(guān)的參數(shù)，θ0是平衡鍵角;Kφ、s、n是與二面角類型相關(guān)的參數(shù)，是二面角的大小;Aij和Bij為作用系數(shù)，不同原子之間的作用系數(shù)通常由來計算，rij是 i原子和j原子之間的距離;ε為介電常數(shù)，qi和qj分別是第i和j個原子的剩余電荷.剩余電荷采用鍵增量δij表示原子j對原子i的剩余電荷.在鍵增量δij中，i表示電荷受體，j表示電荷給體，例如δCH=-0.053，表示碳?xì)滏I中氫原子對碳原子的剩余電荷為-0.053，電荷鍵增量 δij是通過從頭計算靜電勢能得到的.原子i的剩余電荷是所有與其成鍵的電荷鍵增量 δij的總和，即同類型原子間鍵增量為0.模擬過程中用到的參數(shù)如表2～6所示.

表2 鍵長能參數(shù)Table 2 Parameters for bond energy

表3 鍵角能參數(shù)Table 3 Parameters for angle energy

表4 二面角能參數(shù)Table 4 Parameters for dihedral energy

表5 范德華能中的作用系數(shù)Table 5 Parameters for van der Waals(VDW)energy

表6 鍵增量Table 6 Bond increments

1.3 模擬過程

最初，4種聚合物都以大小為0.3 g·cm-3的低密度建立.然后采用模擬退火法，讓模型從500 K溫度開始，每隔50K通過分子動力學(xué)方法弛豫1ns.通過這種方法，使體系的結(jié)構(gòu)逐漸得到優(yōu)化，從而保證模型的結(jié)構(gòu)更加合理.最終4種材料在50 K下的密度分別為0.89、0.92、0.93 和 0.95 g·cm-3.得到材料的初始模型后，繼續(xù)以50 K為間隔，從50 K開始一直模擬到400 K.在每個溫度下，時間步長取0.5 fs，截斷半徑取1.2 nm.每次模擬弛豫1 ×107步，結(jié)果的后一半用于體積，彈性模量和能量的計算.

所有的模擬均在NPT系綜［10］下進(jìn)行，壓強(qiáng)為一個大氣壓.所用的算法［16］結(jié)合了 Nosé-Poincaré控溫算法［17］和 metric-tensor控壓算法［18］.而時間積分則采用generalized leapfrog算法［19］.以上算法都是保辛的，這樣得到的計算結(jié)果才是合理的.另外，計算過程中庫倫能采用Wolf方法［20］計算.Wolf方法比Ewald方法［10］更加簡單，同時又不失準(zhǔn)確性，可以大大節(jié)省計算時間.

2 模擬結(jié)果和討論

2.1 玻璃態(tài)轉(zhuǎn)變溫度

玻璃態(tài)轉(zhuǎn)變溫度Tg是聚合物的一個重要性能指標(biāo).以Tg為界，聚合物呈現(xiàn)不同的物理性質(zhì)，在Tg以下，聚合物處于玻璃態(tài)，在Tg以上，聚合物表現(xiàn)出高彈性質(zhì).Tg有多種測定方法，在分子動力學(xué)模擬中，可以通過體積溫度曲線斜率的轉(zhuǎn)折點(diǎn)來確定［21］.

4種聚合物的體積溫度曲線如圖3所示.圖中，每個溫度下的體積為單位質(zhì)量體積的時間平均值，即平均比容.直線通過最小二乘法擬合得到.從圖中每2條直線的交點(diǎn)，計算得到4種材料的Tg分別為305.5、305.8、307.7 和302.9 K.對于玻璃態(tài)轉(zhuǎn)變溫度，如此微小的差別完全可以忽略.也就是說，POSS的加入對于PE的Tg基本沒有影響.

圖3 體積－溫度曲線Fig.3 Specific volume versus temperature for the polymers

與文獻(xiàn)［22-24］的結(jié)果相比，本文得到的玻璃態(tài)轉(zhuǎn)變溫度顯然偏大.一般說來，聚合物的Tg與結(jié)晶度密切相關(guān)，結(jié)晶度越高，聚合物的Tg越大.而支鏈的存在會對主鏈的結(jié)晶產(chǎn)生影響，支鏈越多結(jié)晶度越差［25］.本文建立的模型中沒有引入支鏈，因此會大大提高材料的結(jié)晶度，從而導(dǎo)致計算得到的Tg偏大.

另外，由于CVFF勢缺少非簡諧振動項與交叉能量項，對聚合物的動力學(xué)模擬結(jié)果并不十分精確.盡管如此，勢函數(shù)對結(jié)果的趨勢并沒有影響，依然可以將其應(yīng)用于后面的計算.

2.2 彈性模量

材料的彈性常數(shù)可以通過應(yīng)力應(yīng)變漲落方法［26］來計算:式中:εij和σij為每一步的應(yīng)變和應(yīng)力，＜* ＞表示系綜平均，在計算中可以用時間平均來代替.通過式(2)可以得到材料在各溫度下的剛度矩陣.由剛度矩陣可以很容易地求出材料的彈性模量［27］.4種聚合物彈性模量隨溫度的變化曲線如圖4所示.

從圖4中可以看出，加入少量POSS后，POSSPE的彈性模量會降低，隨著POSS含量的增加，POSS-PE的彈性模量達(dá)到最大，然后繼續(xù)加入POSS，彈性模量反而下降.在Tg之前，彈性模量的變化都保持這種模式，在Tg之后，聚合物的彈性模量變得很小，POSS對PE彈性模量的影響也變得很不明顯.在300 K時，4種材料的彈性模量分別是0.69、0.67、0.77 和 0.75 GPa.PE 在常溫下的彈性模量約為0.9 GPa左右，與本文得到的結(jié)果比較接近，在可接受的范圍之內(nèi).

圖4 彈性模量隨溫度變化曲線Fig.4 Isotropic elastic constants for the polymers with respect to the temperatures

2.3 能量分析

在模擬的過程中，每一步的能量都可以得到，能量的每一部分也可以清楚地獲得.材料的結(jié)構(gòu)改變十分明顯地表現(xiàn)在能量的變化上，因此，可以通過從能量各部分的變化情況來考察材料結(jié)構(gòu)的變化，進(jìn)而與材料的性質(zhì)產(chǎn)生聯(lián)系.

4種聚合物的各部分能量隨溫度變化的曲線如圖5～9所示，其中，圖7、8中的直線是最小二乘擬合的結(jié)果.從圖中可以看出，鍵長能和鍵角能隨溫度線性變化，在Tg前后，曲線斜率沒有任何改變，說明鍵長和鍵角的變化并不影響材料的Tg.而庫倫能幾乎不隨溫度變化，說明材料的庫倫能不受材料結(jié)構(gòu)變化的影響.另外，可以看出各材料庫倫能之間相差為常數(shù)，約125 eV.這應(yīng)該是加入0.4 mol%的POSS產(chǎn)生的庫倫能改變.在以后的討論中，就不考慮庫倫能的變化.

受Tg影響有明顯變化的是二面角能和范德華能.兩部分能量在Tg前后，曲線斜率均出現(xiàn)了轉(zhuǎn)折點(diǎn)，與體積溫度曲線的情形類似.二面角能的變化代表的是主鏈結(jié)構(gòu)的扭轉(zhuǎn)，而范德華能的變化則表示主鏈之間發(fā)生相對運(yùn)動.二者能量的突變表明，主鏈的扭轉(zhuǎn)和主鏈之間的相對運(yùn)動是聚合物產(chǎn)生玻璃態(tài)轉(zhuǎn)變主要原因.這一結(jié)論與文獻(xiàn)［25］中的理論相符合.

圖5 鍵長能隨溫度變化曲線Fig.5 Bond energies for the polymerswith respect to the temperatures

圖6 鍵角能隨溫度變化曲線Fig.6 Angle energies for the polymerswith respect to the temperatures

圖7 二面角能隨溫度變化曲線Fig.7 Dihedral energies for the polymerswith respect to the temperatures

圖8 范德華能隨溫度變化曲線Fig.8 VDwenergies for the polymers with respect to the temperatures

圖9 庫倫能隨溫度變化曲線圖Fig.9 Coulomb energies for the polymers with respect to the temperatures

由于在所有溫度下，POSS-PE的能量趨勢是相同的，所以本文以50 K為例，分析POSS單體的加入對基體主鏈能量的影響.主鏈的各部分能量與POSS比例的關(guān)系如圖10所示.從圖中可以看出，加入POSS后，POSS-PE主鏈的鍵角能，二面角能和范德華能與PE相比變化很小，基本可以忽略.對于鍵長能，PEIV和PEVI的鍵長能變化不大，而PEII的鍵長能有顯著提高，說明POSS的加入使基體材料主鏈的鍵長被拉長，從而導(dǎo)致結(jié)構(gòu)變得不穩(wěn)定，使基體的彈性模量減小.PEII相對于PE彈性模量的減小就是由此導(dǎo)致的.

圖10 50 K下主鏈能量隨POSS比例的變化曲線Fig.10 Energies of themain chainswith respect to the content of the POSS units at 50 K

圖11 平均每個POSS單體與主鏈之間的范德華能隨POSS比例的變化曲線Fig.11 VDwEnergies between the POSS unit and the main chains on averagewith respect to the content of the POSS units

為了更深入地分析POSS的加入對PE的影響，計算得到了POSS單體與基體材料主鏈之間的范德華能的變化，如圖11所示.從圖中可以看出，300 K以前，平均每個POSS單體與基體之間的范德華能隨POSS比例的增高而增大，在400 K時，能量的變化沒有之前顯著.范德華能的增加說明POSS與主鏈之間的距離增加.而POSS與主鏈距離的增加說明兩者的結(jié)合減弱.PEIV和PEVI彈性模量在300K之前隨POSS比例增加而減小與這個趨勢相符合.在300 K之后，POSS對PE的彈性模量影響不大，也與400 K時能量的變化吻合.因此，材料彈性模量增強(qiáng)與POSS和主鏈之間的相互作用密切相關(guān).由此說來，PEII的增強(qiáng)效果應(yīng)該最強(qiáng)，而PEVI最弱.但是如前面討論過的，由于PEII中POSS與主鏈的強(qiáng)烈吸引，導(dǎo)致主鏈的鍵長增加，反而劣化了基體的結(jié)構(gòu)而導(dǎo)致彈性模量減小，最終的結(jié)果是PEIV對材料彈性模量的增強(qiáng)效果最大.這樣POSS-PE的彈性模量隨著POSS比例的增加而變化就解釋清楚了.

3 結(jié)論

通過分子動力學(xué)模擬得到的結(jié)果，可以得出以下結(jié)論:

1)POSS的加入對于PE的Tg基本沒有影響.聚合物在Tg前后，產(chǎn)生了較大的主鏈扭轉(zhuǎn)和鏈間相對運(yùn)動.這一結(jié)論與已有的理論相符合.

2)隨著POSS比例的增加，彈性模量先出現(xiàn)減小，然后增加，而后再減小.在Tg之前，彈性模量的改變均保持這種模式.而在Tg之后，PE的模量變得很小，POSS的加入對PE彈性模量的改變可以忽略.說明POSS有效地改變了PE在Tg之前的彈性模量.

3)POSS與主鏈之間的范德華作用是導(dǎo)致POSS-PE彈性模量增加的原因，而PEⅡ相對于PE彈性模量的減小，是因為POSS的加入導(dǎo)致主鏈鍵長的增加，而使基體材料結(jié)構(gòu)變差，彈性模量減小造成的.

［1］SANCHEZ C，RIBOT F.Design of hybrid organic-inorganic materials synthesized via Sol-Gel chemistry［J］.NewJournal of Chemistry，1994，18:1007-1014.

［2］MATHER PT，JEONH G，ROMO-URIBEA，etal.Mechanical relaxation andmicrostructure of poly(norbornyl-POSS)copolymers［J］.Macromolecules，1999，32:1194-1203.

［3］FU B X，HSIAO B S，PAGOLA S，et al.Structural development during deformation of polyurethane containing polyhedral oligomeric silsesquioxanes(POSS)molecules［J］.Polymer，2001，42:599-611.

［4］LIU H Z，ZHENG S X.Polyurethane networks nanoreinforced by polyhedral oligomeric silsesquioxan［J］.Macromolecular Rapid Communications，2005，26(3):196-200.

［5］BHARADWAJR K，BERRY R J，F(xiàn)ARMER B L.Molecular dynamics simulation study of norbornene-POSS polymers［J］.Polymer，2000，41:7209-7221.

［6］BIZET S，GALY J，GERARD J F.Molecular dynamics simulation of organic-inorganic copolymers based on methacryl-POSS and methyl methacrylate［J］.Polymer，2006，47:8219-8227.

［7］FINA A，TABUANI D，CAMINO G.Polypropylene-polysilsesquioxane blends［J］.European Polymer Journal，2010，46:14-23.

［8］LEU C，CHANG Y，WEIK.Synthesis and dielectric properties of polyimide-tethered polyhedral oligomeric silsesquioxane(POSS)nanocomposites via POSS-diamine［J］.Macromolecules，2003，36:9122-9127.

［9］TSUCHIDA A，BOLLN C，SERNETZ F G，et al.Ethene and propene copolymers containing silsesquioxane side groups［J］.Macromolecules，1997，30:2818-2824.

［10］FRENKEL D，SMIT B.分子模擬-從算法到應(yīng)用［M］.北京:化學(xué)工業(yè)出版社，2002:27-29.

［11］MAPLE J R，DINUR U，HAGLER A T.Derivation of force fields formolecularmechanics and dynamics fromab initio energy surfaces［C］//Proc Nati Acad Sci.［s.l.］，1988，85:5350-5354.

［12］LAUTENSCHL?GER P，BRICKMANN J，VAN RUITEN J.Conformation and rotational barriers of aromatic polyesters［J］.Macromolecules，1991，24:1284-1292.

［13］BHOWMIK R，KATTIK S，KATTID.Molecular dynamics simulation of hydroxyapatite-epolyacrylic acid interfaces［J］.Polymer，2007，48:664-674.

［14］GREATHOUSE JA，ALLENDORFmD.Force field validation formolecular dynamics simulations of IRMOF-1 and other isoreticular zinc carboxylate coordination polymers［J］.JPhys ChemC，2008，112:5795-5802.

［15］HEINZ H，KOERNER H，ANDERSON K L，et al.Force field for mica-type silicates and dynamics of octadecylammoniumchains grafted to montmorillonite［J］.ChemMater，2005，17:5658-5669.

［16］HERNANDEZ E.Metric-tensor flexible-cell algorithmfor isothermal-isobaric molecular dynamics simulations［J］.J ChemPhys，2001，115(22):10282-10290.

［17］BOND SD，LEIMKUHLER B J，LAIRD B B.The nosepoincaremethod for constant temperaturemolecular dynamics［J］.JComput Phys，1999，151:114-134.

［18］SOUZA I，MARTINS J L.Metric tensor as the dynamical variable for variable-cell-shape molecular dynamics［J］.Phys Rev B，1997，55(14):8733-8742.

［19］SUN G.Construction of high order symplectic Runge-Kutta methods［J］.JComput Math，1993，11(3):250-260.

［20］WOLF D，KEBLINSKIP，PHILLPOT S R，et al.Exact method for the simulation of coulombic systems by spherically truncated，pairwise r-1 summation［J］.J ChemPhys，1999，110(17):8254-8282.

［21］HAN JC，GEE R H，BOYD R H.Glass transition temperatures of polymers frommolecular dynamics simulations［J］.Macromolecules，1994，27:7781-7784.

［22］GAUR U，WUNDERLICH B.The glass transition temperature of polyethylene［J］.Macromolecules，1980，13:445-446.

［23］NG SC，HOSEA T JC，GOH SH.Glass transition of polyethylene as studied by brillouin spectroscopy［J］.Polymer Bulletin，1987，18:155-158.

［24］EBERLE R，ANDERSSH，WEISHAUPT K，et al.Anisotropic effects of the glass transition in oriented polyethylene［J］.Europhys Lett，1998，43:201-206.

［25］何平笙.新編高聚物的結(jié)構(gòu)與性能［M］.北京:科學(xué)出版社，2009:159-182.

［26］CUIZ，SUN Y，LI J，et al.Combination method for the calculation of elastic constants［J］.Phys Rev B，2007，75(21):214101.1-214101.6.

［27］HADDADIK，BOUHEMADOU A，LOUAIL L，et al.Ab initio study of structural，elastic，and high-pressure properties of rubidiumhalides RbX(X=F，Cl，Br and I)［J］.Phase Transit，2009，82(3):266-279.