王 巍， 王云婷， 李新寧

(東北林業(yè)大學(xué) 工程技術(shù)學(xué)院， 黑龍江 哈爾濱 150040)

隨著科技發(fā)展及人們對(duì)節(jié)能環(huán)保材料的需求，木材因較低的成本、可再生性和較好的機(jī)械性能備受關(guān)注，而木材的熱處理過(guò)程對(duì)其質(zhì)量起著至關(guān)重要的影響。目前已經(jīng)有許多中外學(xué)者對(duì)木材本身的熱處理過(guò)程進(jìn)行研究，發(fā)現(xiàn)木材熱處理能夠改變木材的化學(xué)成分，并提高其尺寸穩(wěn)定性[1]。木材的動(dòng)態(tài)彈性模量隨著熱處理溫度的升高和處理時(shí)間的延長(zhǎng)而逐漸降低[2]，熱處理對(duì)木材彈性模量的影響隨著溫度的升高逐漸加大，對(duì)剪切模量則沒(méi)有很大影響[3]；高溫?zé)崽幚硪矔?huì)導(dǎo)致木材力學(xué)強(qiáng)度減小和質(zhì)量損失[4]。因此，從微觀層面了解木材的熱處理過(guò)程可以更好地理解其宏觀性能的變化。

纖維素是木材的主要成分之一，大量存在于起著支撐作用的細(xì)胞壁中。纖維素微纖維的高軸向模量和強(qiáng)度對(duì)植物細(xì)胞壁的機(jī)械性能有著很大的貢獻(xiàn)，木材中的纖維也是從微纖維中獲得的拉伸強(qiáng)度和模量。到目前為止，研究者發(fā)現(xiàn)纖維素具有5種結(jié)晶形式。通過(guò)13C CP/MAS固態(tài)核磁共振實(shí)驗(yàn)發(fā)現(xiàn)，天然纖維素的結(jié)晶體為纖維素I，并且具有Iα和Iβ2種結(jié)晶相[5]。Iα結(jié)晶相主要存在于細(xì)菌和海藻中，Iβ結(jié)晶相主要存在于高等植物和動(dòng)物被膜中。木材中的纖維素為纖維素Iβ結(jié)晶相，因此本研究選取纖維素Iβ作為研究對(duì)象。

分子動(dòng)力學(xué)(MD)模擬是一種從微觀層面研究物質(zhì)力學(xué)性能的方法。齊曉飛等[6]采用MD模擬方法研究了硝化纖維素/硝化甘油共存體系力學(xué)性能隨溫度變化的關(guān)系。陳芳等[7]采用MD模擬方法分析了高聚物力學(xué)性能和結(jié)合能的變化規(guī)律。MD模擬可以成熟地建立宏/微觀聯(lián)系，但是在木材熱處理方面鮮有通過(guò)該法對(duì)木材內(nèi)部組織結(jié)構(gòu)及力學(xué)性能的研究。本研究通過(guò)使用Material Studio 7.0軟件構(gòu)建纖維素Iβ模型，分析了不同溫度下的結(jié)晶纖維素結(jié)構(gòu)和力學(xué)參數(shù)，以期為理解木材熱處理過(guò)程提供理論依據(jù)和性能預(yù)測(cè)。

1 模型構(gòu)建和動(dòng)力學(xué)模擬

1.1 纖維素Iβ模型構(gòu)建

纖維素的基本單位是脫水葡萄糖，其化學(xué)結(jié)構(gòu)是由8 000到10 000個(gè)D-吡喃葡萄糖環(huán)以β-1,4-糖苷鍵聯(lián)結(jié)而成的線形高分子化合物[8]，重復(fù)單位是纖維二糖。纖維素Iβ鏈?zhǔn)瞧叫信帕械?，通過(guò)鏈間OH…O氫鍵形成穩(wěn)定的片狀結(jié)構(gòu)，氫鍵的存在使得纖維素分子鏈的結(jié)構(gòu)更加牢固[8]。

根據(jù)Nishiyama通過(guò)X射線和中子衍射數(shù)據(jù)，纖維素Iβ的模型[9]參數(shù)為：a=0.778 4 nm，b=0.820 1 nm，c=1.038 nm，α=β=90°，γ=96.5°。使用Build Crystal模塊對(duì)Iβ晶體纖維素采用周期性邊界條件進(jìn)行建模，并構(gòu)建了3×3×3的超晶胞結(jié)構(gòu)，空間群結(jié)構(gòu)為單斜晶P21。纖維素Iβ的晶體結(jié)構(gòu)和晶格參數(shù)如圖1所示。

a. x-y平面x-y plane; b. y-z平面y-z plane; c. x-z平面x-z plane圖1 纖維素Iβ模型及晶格參數(shù)Fig.1 Cellulose Iβ model and lattice parameters

1.2 動(dòng)力學(xué)模擬

1.2.1模擬設(shè)置 使用Accelrys Materials Studio 7.0中的Forcite Plus模塊進(jìn)行動(dòng)力學(xué)模擬，該模塊是一個(gè)經(jīng)典分子力學(xué)工具，可以進(jìn)行幾何優(yōu)化、淬火、退火、動(dòng)力學(xué)計(jì)算、機(jī)械性能計(jì)算等。Forcite Plus模塊中的Calculation工具可以對(duì)得到的結(jié)構(gòu)進(jìn)行力學(xué)性能計(jì)算，得到彈性模量、剪切模量、體積模量、剛度矩陣等力學(xué)參數(shù)。Analysis工具可對(duì)結(jié)構(gòu)進(jìn)行分析，能得到原子徑向分布(RDF)、X射線衍射圖、均方位移等圖像。通過(guò)使用Forcite Plus可以計(jì)算并預(yù)測(cè)出不同溫度下的晶胞參數(shù)、密度、力學(xué)性能。該模塊支持COMPASS、cvff、pcff、universal、Dreiding等力場(chǎng)，COMPASS力場(chǎng)能準(zhǔn)確地模擬出分子和凝聚態(tài)的結(jié)構(gòu)和熱力學(xué)性質(zhì)[10]，因此本研究選擇COMPASS力場(chǎng)作為模擬力場(chǎng)。

模型優(yōu)化部分，本研究選用Smart方法對(duì)模型進(jìn)行5 000步的結(jié)構(gòu)優(yōu)化，使原子自由運(yùn)動(dòng)趨于平衡，達(dá)到能量最小化。在500 K溫度下應(yīng)用正則(NVT，即原子數(shù)N、體積V、溫度T保持不變)系綜進(jìn)行時(shí)間步長(zhǎng)為1 fs、總時(shí)間長(zhǎng)200 ps的模擬。最小化過(guò)程中，空間群結(jié)構(gòu)由P21變?yōu)镻1。溫度采用Nose控制方法，壓強(qiáng)采用Berendsen控制方法，電荷相互作用加和方法采用Ewald方法，范德華力加和方法采用基于原子計(jì)算方法，截?cái)喟霃皆O(shè)置為1.25 nm，每5 000步輸出一個(gè)結(jié)構(gòu)。經(jīng)過(guò)該優(yōu)化步驟，模型已經(jīng)實(shí)現(xiàn)初步弛豫，消除了局部應(yīng)力，能量達(dá)到最小化。

分子模擬部分，本研究選用等溫等壓(NPT，即原子數(shù)N、壓強(qiáng)P、溫度T保持不變)系綜來(lái)完成各溫度下的模擬。設(shè)置初始速度為基于蒙塔卡羅算法的隨機(jī)速度，進(jìn)行時(shí)間步長(zhǎng)為1 fs[11]，時(shí)間總長(zhǎng)為 1 ns 的模擬。根據(jù)木材熱處理的溫度范圍，設(shè)置模擬溫度為350～550 K，溫度間隔為20 K。

1.2.2平衡判定 為確定模擬體系達(dá)到穩(wěn)定且具有統(tǒng)計(jì)熱力學(xué)意義，分子動(dòng)力學(xué)模擬中通常認(rèn)為溫度與能量波動(dòng)值在5%～10%時(shí)即可認(rèn)為當(dāng)前體系趨衡[12]。因此在模擬體系中進(jìn)行了1 ns 的模擬，得到了溫度和能量隨時(shí)間的變化圖(圖2和圖3)。圖中能夠看出模擬過(guò)程中溫度和能量波動(dòng)都很小，表明該體系經(jīng)過(guò)分子動(dòng)力學(xué)模擬已經(jīng)平衡，建立的模型可用于接下來(lái)的結(jié)構(gòu)分析和計(jì)算。

能量收斂參數(shù)(δE)能夠用于衡量模擬的精度，當(dāng)δE≤ 0.001～0.003時(shí)，模擬結(jié)果可靠，計(jì)算公式如下：

式中：N―分子模擬過(guò)程中設(shè)定的總步數(shù)；E0―模擬體系的初始能量；Ei―模擬體系經(jīng)i步后的能量。經(jīng)計(jì)算得出本研究模擬體系的δE= 0.001 9，處于可靠的范圍當(dāng)中，證明了該體系能量收斂可靠。

圖2 溫度-時(shí)間變化圖Fig.2 Temperature-time variation graph

圖3 能量-時(shí)間變化圖Fig.3 Energy-time variation graph

2 結(jié)果與分析

2.1 晶胞體積及密度變化

a. 350 K; b. 550 K圖4 不同溫度下纖維素Iβ結(jié)構(gòu)Fig.4 Cellulose Iβ structure at different temperatures

為了探討溫度對(duì)內(nèi)部分子結(jié)構(gòu)的影響，本研究模擬了350 K(熱處理起始溫度)和550 K(熱處理最高溫度)下的纖維素Iβ結(jié)構(gòu)(見(jiàn)圖4)。由圖4可知，350 K時(shí)纖維素分子鏈維持相對(duì)規(guī)整的狀態(tài)，經(jīng)550 K高溫弛豫后分子鏈發(fā)生明顯扭曲，可以看出高溫破壞了維持分子鏈穩(wěn)定的作用力。

由圖5可知，纖維素Iβ的內(nèi)部結(jié)構(gòu)受溫度影響，纖維素Iβ的體積隨著溫度的升高而逐漸增大。隨著溫度的升高，鏈間氫鍵出現(xiàn)了斷裂，減弱了對(duì)纖維素分子鏈排布的約束作用。同時(shí)，使片層穩(wěn)定疊積的弱范德華力減小也使得纖維素鏈間距離增加。

密度是用于評(píng)估模型的一個(gè)參數(shù)，NPT系綜允許晶胞參數(shù)發(fā)生變化，因此在不同的模擬條件下，晶胞的參數(shù)也會(huì)發(fā)生變化。350 K下模擬的纖維素Iβ密度為1.617 g/cm3，550 K下模擬的密度為1.581 g/cm3。文獻(xiàn)中試驗(yàn)測(cè)定的密度在1.543～1.643 g/cm3范圍[13-15]，和本研究模擬測(cè)定的基本相符，說(shuō)明構(gòu)建的模型是可行的。伴隨著密度的減小的同時(shí)，晶胞體積由11.99 nm3增加至12.26 nm3。

2.2 氫鍵變化

氫鍵是一種較強(qiáng)的非鍵作用，形成氫鍵的原子之間既有范德華力又有庫(kù)侖力，大小介于化學(xué)鍵作用和非鍵作用之間，目前為止尚未統(tǒng)一界定氫鍵形成的標(biāo)準(zhǔn)。通常將和氫原子成鍵的原子叫做供體，將和氫原子相連的原子稱為受體。采用判定氫鍵存在的幾何標(biāo)準(zhǔn)中，氫原子的受體的距離不大于0.28 nm[16]，供體、氫原子和受體的角度不小于110°。氫鍵的狀態(tài)和分布對(duì)纖維素的結(jié)構(gòu)有著重大的意義，因此分析氫鍵的數(shù)量變化對(duì)理解纖維素Iβ是必不可少的。采用Material Studio軟件的計(jì)算氫鍵模塊進(jìn)行氫鍵的判斷，通過(guò)Perl腳本進(jìn)行氫鍵數(shù)量的統(tǒng)計(jì)。分子模型中的氫鍵分布如圖6所示，不同溫度下纖維素Iβ的氫鍵數(shù)量如表1所示。

圖5 不同溫度下纖維素Iβ的密度和體積變化Fig.5 Density and volume change of cellulose Iβ at different temperatures

圖6 分子模型中的氫鍵分布(虛線表示氫鍵)Fig.6 Hydrogen bond distribution in the molecular model(dotted lines indicate hydrogen bonding)

分子鏈內(nèi)氫鍵維持分子鏈穩(wěn)定性，而分子間氫鍵維持片層穩(wěn)定性。由表1可知，纖維素Iβ的氫鍵總數(shù)量隨著溫度增大而逐漸減小，350 K到550 K之間氫鍵總數(shù)減少了24%，這表示溫度對(duì)氫鍵影響不大，氫鍵具有一定的熱穩(wěn)定性。隨著溫度的升高，分子鏈內(nèi)氫鍵占?xì)滏I總數(shù)的百分比逐漸減小，由350 K時(shí)的68%降低至550 K時(shí)的40%，其中390～450 K范圍內(nèi)鏈內(nèi)氫鍵大幅減少，鏈間氫鍵大幅增加。這是由于溫度的升高導(dǎo)致鏈內(nèi)氫鍵破裂而形成了新的鏈間氫鍵。鏈內(nèi)氫鍵與鏈間氫鍵的比值由350 K時(shí)2.1∶1變成550 K時(shí)1∶1.5。高溫重構(gòu)了氫鍵網(wǎng)絡(luò)使纖維素分子鏈更易發(fā)生扭曲，而鏈間穩(wěn)定性得到增強(qiáng)，印證了2.1節(jié)纖維素Iβ的結(jié)構(gòu)和體積變化。

表1 不同溫度下纖維素Iβ氫鍵數(shù)目

2.3 力學(xué)性能

MD模擬能夠計(jì)算不同溫度下纖維素Iβ的力學(xué)性能。楊氏模量(E)能夠說(shuō)明材料抵抗形變的能力，該值越大，則材料剛性越大，越不容易產(chǎn)生變形;剪切模量(G)是材料在剪切情況下，應(yīng)力與應(yīng)變的比值;體積模量(K)指在均勻壓縮條件下，材料所受應(yīng)力和應(yīng)變的比值。通過(guò)Materials Studio軟件模擬能夠直接得到上述力學(xué)參數(shù)，如表2所示。

表2 各溫度下的纖維素Iβ力學(xué)參數(shù)

判斷變量之間是否相關(guān)可使用皮爾遜相關(guān)系數(shù)(R)進(jìn)行計(jì)算，計(jì)算公式如下：

式中：n—數(shù)據(jù)個(gè)數(shù)；x—時(shí)間；y—力學(xué)參數(shù)。

相關(guān)系數(shù)定義規(guī)定：若0.7≤|R|≤1，則兩變量為高度線性相關(guān)；若0.4<|R|<0.7為顯著性相關(guān)；若|R|≤0.4為低線性相關(guān)[17]。

由表2可以看出，楊氏模量隨著溫度的升高有減小趨勢(shì)，但是變化不大，變化率約為13%。在390～450 K區(qū)間內(nèi)，楊氏模量呈直線下降，經(jīng)計(jì)算得出R=0.814 7為高度線性相關(guān)。因此溫度對(duì)楊氏模量的影響顯著，這種高度相關(guān)性是2.2節(jié)所述分子鏈內(nèi)氫鍵的減少導(dǎo)致的。剪切模量、體積模量的R值分別為-0.065 1和0.274 0，與溫度呈低線性相關(guān)。總體來(lái)說(shuō)溫度的變化對(duì)纖維素Iβ影響比較小。

3 結(jié) 論

3.1通過(guò)使用Material Studio分子動(dòng)力學(xué)模擬軟件對(duì)木質(zhì)纖維素Iβ進(jìn)行建模，并在等溫等壓(NPT)系綜不同溫度下進(jìn)行分子動(dòng)力學(xué)模擬。結(jié)果表明：高溫會(huì)使纖維素Iβ分子鏈活動(dòng)加劇，晶格體積增大，由350 K的11.99 nm3增加至550 K的12.26 nm3。模型密度變化范圍(1.581～1.617 g/cm3)和實(shí)驗(yàn)得到的密度范圍(1.543～1.643 g/cm3)相符，驗(yàn)證了模型的可行性。

3.2計(jì)算并統(tǒng)計(jì)了不同溫度下的氫鍵數(shù)目，結(jié)果發(fā)現(xiàn)：氫鍵隨著溫度的升高而逐漸減少，但是整體變化不大(350～550 K氫鍵總數(shù)減少24%)，體現(xiàn)了纖維素Iβ具有一定的熱穩(wěn)定性。分子鏈內(nèi)氫鍵部分?jǐn)嗔讯纬闪随滈g氫鍵，印證了纖維素Iβ的結(jié)構(gòu)和體積變化。

3.3在350～550 K木材熱處理的溫度范圍內(nèi)計(jì)算了纖維素Iβ的力學(xué)參數(shù)，氫鍵網(wǎng)絡(luò)的重新排布改變了纖維素的力學(xué)性能，溫度的升高使楊氏模量逐漸減小，變化率約為13%。相比于楊氏模量，剪切模量和體積模量受溫度的影響較小，基本處于在一個(gè)范圍內(nèi)波動(dòng)的狀態(tài)，沒(méi)有明顯的變化趨勢(shì)。