热处理对木质纤维素Iβ结构及其力学性能影响的分子动力学模拟

王 巍， 王云婷， 李新宁

(东北林业大学 工程技术学院， 黑龙江 哈尔滨 150040)

随着科技发展及人们对节能环保材料的需求，木材因较低的成本、可再生性和较好的机械性能备受关注，而木材的热处理过程对其质量起着至关重要的影响。目前已经有许多中外学者对木材本身的热处理过程进行研究，发现木材热处理能够改变木材的化学成分，并提高其尺寸稳定性[1]。木材的动态弹性模量随着热处理温度的升高和处理时间的延长而逐渐降低[2]，热处理对木材弹性模量的影响随着温度的升高逐渐加大，对剪切模量则没有很大影响[3]；高温热处理也会导致木材力学强度减小和质量损失[4]。因此，从微观层面了解木材的热处理过程可以更好地理解其宏观性能的变化。

纤维素是木材的主要成分之一，大量存在于起着支撑作用的细胞壁中。纤维素微纤维的高轴向模量和强度对植物细胞壁的机械性能有着很大的贡献，木材中的纤维也是从微纤维中获得的拉伸强度和模量。到目前为止，研究者发现纤维素具有5种结晶形式。通过13C CP/MAS固态核磁共振实验发现，天然纤维素的结晶体为纤维素I，并且具有Iα和Iβ2种结晶相[5]。Iα结晶相主要存在于细菌和海藻中，Iβ结晶相主要存在于高等植物和动物被膜中。木材中的纤维素为纤维素Iβ结晶相，因此本研究选取纤维素Iβ作为研究对象。

分子动力学(MD)模拟是一种从微观层面研究物质力学性能的方法。齐晓飞等[6]采用MD模拟方法研究了硝化纤维素/硝化甘油共存体系力学性能随温度变化的关系。陈芳等[7]采用MD模拟方法分析了高聚物力学性能和结合能的变化规律。MD模拟可以成熟地建立宏/微观联系，但是在木材热处理方面鲜有通过该法对木材内部组织结构及力学性能的研究。本研究通过使用Material Studio 7.0软件构建纤维素Iβ模型，分析了不同温度下的结晶纤维素结构和力学参数，以期为理解木材热处理过程提供理论依据和性能预测。

1 模型构建和动力学模拟

1.1 纤维素Iβ模型构建

纤维素的基本单位是脱水葡萄糖，其化学结构是由8 000到10 000个D-吡喃葡萄糖环以β-1,4-糖苷键联结而成的线形高分子化合物[8]，重复单位是纤维二糖。纤维素Iβ链是平行排列的，通过链间OH…O氢键形成稳定的片状结构，氢键的存在使得纤维素分子链的结构更加牢固[8]。

根据Nishiyama通过X射线和中子衍射数据，纤维素Iβ的模型[9]参数为：a=0.778 4 nm，b=0.820 1 nm，c=1.038 nm，α=β=90°，γ=96.5°。使用Build Crystal模块对Iβ晶体纤维素采用周期性边界条件进行建模，并构建了3×3×3的超晶胞结构，空间群结构为单斜晶P21。纤维素Iβ的晶体结构和晶格参数如图1所示。

a. x-y平面x-y plane; b. y-z平面y-z plane; c. x-z平面x-z plane图1 纤维素Iβ模型及晶格参数Fig.1 Cellulose Iβ model and lattice parameters

1.2 动力学模拟

1.2.1模拟设置 使用Accelrys Materials Studio 7.0中的Forcite Plus模块进行动力学模拟，该模块是一个经典分子力学工具，可以进行几何优化、淬火、退火、动力学计算、机械性能计算等。Forcite Plus模块中的Calculation工具可以对得到的结构进行力学性能计算，得到弹性模量、剪切模量、体积模量、刚度矩阵等力学参数。Analysis工具可对结构进行分析，能得到原子径向分布(RDF)、X射线衍射图、均方位移等图像。通过使用Forcite Plus可以计算并预测出不同温度下的晶胞参数、密度、力学性能。该模块支持COMPASS、cvff、pcff、universal、Dreiding等力场，COMPASS力场能准确地模拟出分子和凝聚态的结构和热力学性质[10]，因此本研究选择COMPASS力场作为模拟力场。

模型优化部分，本研究选用Smart方法对模型进行5 000步的结构优化，使原子自由运动趋于平衡，达到能量最小化。在500 K温度下应用正则(NVT，即原子数N、体积V、温度T保持不变)系综进行时间步长为1 fs、总时间长200 ps的模拟。最小化过程中，空间群结构由P21变为P1。温度采用Nose控制方法，压强采用Berendsen控制方法，电荷相互作用加和方法采用Ewald方法，范德华力加和方法采用基于原子计算方法，截断半径设置为1.25 nm，每5 000步输出一个结构。经过该优化步骤，模型已经实现初步弛豫，消除了局部应力，能量达到最小化。

分子模拟部分，本研究选用等温等压(NPT，即原子数N、压强P、温度T保持不变)系综来完成各温度下的模拟。设置初始速度为基于蒙塔卡罗算法的随机速度，进行时间步长为1 fs[11]，时间总长为 1 ns 的模拟。根据木材热处理的温度范围，设置模拟温度为350～550 K，温度间隔为20 K。

1.2.2平衡判定 为确定模拟体系达到稳定且具有统计热力学意义，分子动力学模拟中通常认为温度与能量波动值在5%～10%时即可认为当前体系趋衡[12]。因此在模拟体系中进行了1 ns 的模拟，得到了温度和能量随时间的变化图(图2和图3)。图中能够看出模拟过程中温度和能量波动都很小，表明该体系经过分子动力学模拟已经平衡，建立的模型可用于接下来的结构分析和计算。

能量收敛参数(δE)能够用于衡量模拟的精度，当δE≤ 0.001～0.003时，模拟结果可靠，计算公式如下：

式中：N―分子模拟过程中设定的总步数；E0―模拟体系的初始能量；Ei―模拟体系经i步后的能量。经计算得出本研究模拟体系的δE= 0.001 9，处于可靠的范围当中，证明了该体系能量收敛可靠。

图2 温度-时间变化图Fig.2 Temperature-time variation graph

图3 能量-时间变化图Fig.3 Energy-time variation graph

2 结果与分析

2.1 晶胞体积及密度变化

a. 350 K; b. 550 K图4 不同温度下纤维素Iβ结构Fig.4 Cellulose Iβ structure at different temperatures

为了探讨温度对内部分子结构的影响，本研究模拟了350 K(热处理起始温度)和550 K(热处理最高温度)下的纤维素Iβ结构(见图4)。由图4可知，350 K时纤维素分子链维持相对规整的状态，经550 K高温弛豫后分子链发生明显扭曲，可以看出高温破坏了维持分子链稳定的作用力。

由图5可知，纤维素Iβ的内部结构受温度影响，纤维素Iβ的体积随着温度的升高而逐渐增大。随着温度的升高，链间氢键出现了断裂，减弱了对纤维素分子链排布的约束作用。同时，使片层稳定叠积的弱范德华力减小也使得纤维素链间距离增加。

密度是用于评估模型的一个参数，NPT系综允许晶胞参数发生变化，因此在不同的模拟条件下，晶胞的参数也会发生变化。350 K下模拟的纤维素Iβ密度为1.617 g/cm3，550 K下模拟的密度为1.581 g/cm3。文献中试验测定的密度在1.543～1.643 g/cm3范围[13-15]，和本研究模拟测定的基本相符，说明构建的模型是可行的。伴随着密度的减小的同时，晶胞体积由11.99 nm3增加至12.26 nm3。

2.2 氢键变化

氢键是一种较强的非键作用，形成氢键的原子之间既有范德华力又有库仑力，大小介于化学键作用和非键作用之间，目前为止尚未统一界定氢键形成的标准。通常将和氢原子成键的原子叫做供体，将和氢原子相连的原子称为受体。采用判定氢键存在的几何标准中，氢原子的受体的距离不大于0.28 nm[16]，供体、氢原子和受体的角度不小于110°。氢键的状态和分布对纤维素的结构有着重大的意义，因此分析氢键的数量变化对理解纤维素Iβ是必不可少的。采用Material Studio软件的计算氢键模块进行氢键的判断，通过Perl脚本进行氢键数量的统计。分子模型中的氢键分布如图6所示，不同温度下纤维素Iβ的氢键数量如表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)

分子链内氢键维持分子链稳定性，而分子间氢键维持片层稳定性。由表1可知，纤维素Iβ的氢键总数量随着温度增大而逐渐减小，350 K到550 K之间氢键总数减少了24%，这表示温度对氢键影响不大，氢键具有一定的热稳定性。随着温度的升高，分子链内氢键占氢键总数的百分比逐渐减小，由350 K时的68%降低至550 K时的40%，其中390～450 K范围内链内氢键大幅减少，链间氢键大幅增加。这是由于温度的升高导致链内氢键破裂而形成了新的链间氢键。链内氢键与链间氢键的比值由350 K时2.1∶1变成550 K时1∶1.5。高温重构了氢键网络使纤维素分子链更易发生扭曲，而链间稳定性得到增强，印证了2.1节纤维素Iβ的结构和体积变化。

表1 不同温度下纤维素Iβ氢键数目

2.3 力学性能

MD模拟能够计算不同温度下纤维素Iβ的力学性能。杨氏模量(E)能够说明材料抵抗形变的能力，该值越大，则材料刚性越大，越不容易产生变形;剪切模量(G)是材料在剪切情况下，应力与应变的比值;体积模量(K)指在均匀压缩条件下，材料所受应力和应变的比值。通过Materials Studio软件模拟能够直接得到上述力学参数，如表2所示。

表2 各温度下的纤维素Iβ力学参数

判断变量之间是否相关可使用皮尔逊相关系数(R)进行计算，计算公式如下：

式中：n—数据个数；x—时间；y—力学参数。

相关系数定义规定：若0.7≤|R|≤1，则两变量为高度线性相关；若0.4<|R|<0.7为显著性相关；若|R|≤0.4为低线性相关[17]。

由表2可以看出，杨氏模量随着温度的升高有减小趋势，但是变化不大，变化率约为13%。在390～450 K区间内，杨氏模量呈直线下降，经计算得出R=0.814 7为高度线性相关。因此温度对杨氏模量的影响显著，这种高度相关性是2.2节所述分子链内氢键的减少导致的。剪切模量、体积模量的R值分别为-0.065 1和0.274 0，与温度呈低线性相关。总体来说温度的变化对纤维素Iβ影响比较小。

3 结 论

3.1通过使用Material Studio分子动力学模拟软件对木质纤维素Iβ进行建模，并在等温等压(NPT)系综不同温度下进行分子动力学模拟。结果表明：高温会使纤维素Iβ分子链活动加剧，晶格体积增大，由350 K的11.99 nm3增加至550 K的12.26 nm3。模型密度变化范围(1.581～1.617 g/cm3)和实验得到的密度范围(1.543～1.643 g/cm3)相符，验证了模型的可行性。

3.2计算并统计了不同温度下的氢键数目，结果发现：氢键随着温度的升高而逐渐减少，但是整体变化不大(350～550 K氢键总数减少24%)，体现了纤维素Iβ具有一定的热稳定性。分子链内氢键部分断裂而形成了链间氢键，印证了纤维素Iβ的结构和体积变化。

3.3在350～550 K木材热处理的温度范围内计算了纤维素Iβ的力学参数，氢键网络的重新排布改变了纤维素的力学性能，温度的升高使杨氏模量逐渐减小，变化率约为13%。相比于杨氏模量，剪切模量和体积模量受温度的影响较小，基本处于在一个范围内波动的状态，没有明显的变化趋势。