氮气在溴化丁基橡胶中渗透行为的分子模拟

席彤彤，张书槐，李幸围，熊 英，郭少云

(高分子材料工程国家重点实验室 四川大学高分子研究所 四川省橡塑材料复合成型技术工程实验室，四川 成都 610065)

溴化丁基橡胶因其优异的气体阻隔性、与其他橡胶良好的相容性和粘接性而广泛应用于各种轮胎的内衬层，从而起到防止气体泄漏，保证轮胎内压的功能[1～3]。目前工业上测量轮胎气密性的方法是直接在轮胎内充入一定压力的气体，测量内压随时间的变化[4]。不同使用环境下轮胎的充气压力不同，如家用汽车轮胎的内压一般约为0.2～0.25 MPa，而飞机轮胎的内压高达1.3～1.5 MPa。上述方法显然难以在实验室条件下实现。而实验室条件下测量聚合物阻隔性能的仪器一般执行GB/T 1038-2000，ASTM D1434，ISO 2556，ISO 15105-2，YBB 00082003 等标准，其测量压力为1个大气压(0.1 MPa)，无法达到或接近航空轮胎内压的高压力来测量材料的气体阻隔性。而高压下测量聚合物阻隔性的仪器需要特殊定制，定制周期长、成本高。与之相比，分子模拟技术能够直观模拟气体在聚合物中的溶解和扩散过程，大大缩短了研究周期。

分子模拟法从原子分子尺度研究体系的微观结构和相互作用，进而得到材料的宏观性能，在药物研发、新材料设计、高分子等领域具有非常广泛的应用[5]。自20世纪90年代以来，研究小分子在聚合物中的渗透行为，如O2，N2，H2，CO2，水，乙醇等小分子在聚硅氧烷类、聚酰胺类、聚酰亚胺类高分子中的扩散行为已逐渐展开并活跃起来[6～10]。因此，为了研究不同压力下气体分子在BIIR中的渗透系数并探索渗透机理，本文采用目前运用最广的美国Accelrys 公司研发的Materials Studio 软件分别以分子力学对初始结构进行优化、以分子动力学和分子蒙特卡洛法对N2在BIIR中的吸附和扩散过程进行模拟，详细研究了N2在BIIR中的渗透过程以及气体压力对其渗透行为的影响，并结合试验数据，验证了分子模拟法在材料气体阻隔性能测试方面的可行性。

1 模拟过程

1.1 计算参数设置

本文使用Materials Studio 软件通过分子力学(MM)、分子动力学(MD)和巨正则蒙特卡洛(GMMC)等方法研究在40 ℃，0.1 MPa 和40 ℃，1.5 MPa 下氮气在溴化丁基橡胶中的溶解和扩散行为，主要采用的软件模块包括：（1）Visualizer模块，用于模型构建；（2）Amorphous Cell模块，用于BIIR无定形晶胞的构建；（3）Sorption 模块，计算氮气在BIIR 中的吸附和溶解；（4）Forcite模块，用于模型的能量计算、几何优化、退火模拟和分子动力学模拟。

在分子模拟过程中时，无论是采用分子力学MM、MD或GMMC方法进行模拟，均需要对计算所采用的力场及电荷计算方式、范德华力计算方式、静电力计算方式等进行设置，本文对相关计算参数的设置如Tab.1所示。

Tab.1 Settings of parameters related to energy calculation

1.2 分子模型构建

本文采用的溴化丁基橡胶为德国朗盛公司生产的BIIR 2030，其分子结构中溴的质量分数为1.8%±0.2%，不饱和度约1.9%。溴化丁基橡胶是丁基橡胶溴化之后的产物，质子核磁共振结构表明，经溴化反应后，约60%和20%的异戊二烯成为仲烯丙基溴化异戊二烯(SABI)和伯烯丙基溴化异戊二烯(PABI)，另外约有20%的异戊二烯未被溴化[11]，如Fig.1所示。

根据上述BIIR的结构，通过Visualizer模块构建4 种结构单元并几何优化，使其达到低能量稳定状态。氮气分子的构造过程类似。之后通过Build Polymers 选项将结构单元首尾连接无规共聚，建立聚合度分别为50，100，200，300 的BIIR 分子链并进行几何优化。BIIR 分子链构建完成后，在Amorphous Cell模块的Construction下分别建立4个密度为0.93 g/cm3的无定形晶胞。由于后续计算中，范德华力的截断半径取值为1.25 nm，为了避免BIIR模型中的原子与自己的周期性镜像发生作用而导致计算结果错误，在构建BIIR 无定形模型时，需对不同聚合度BIIR 链的个数进行调节，聚合度为50，100，200，300 的无定形晶胞中放置链的个数分别设为3，2，2，1，使得最终构建BIIR 无定形模型的长宽高均大于截断半径的2 倍，最终得到不同聚合度的BIIR无定形晶胞，如Fig.2所示。

Fig.2 BIIR amorphous cells with polymerization degree of(a)50,(b)100,(c)200 and(d)300

聚合物的密度对气体阻隔性有较大影响。为了选取密度最接近BIIR 真实密度的模型作为后续研究的模型，需要对上述4 种不同聚合度的无定形晶胞进行几何优化、退火及分子动力学模拟。几何优化选择Smart 算法，能量收敛判据设置为4.184×10−3kJ/mol。在几何优化后需要对体系进行退火模拟，目的在于使分子结构不断调整构型，释放内应力，逐渐接近能量最低的真实结构。进行退火模拟时。初末温度分别为300 K 和600 K，循环5 次，升、降温度间隔为50 K，每个温度下运行10 ps的分子动力学平衡。最后将退火后的模型进行NPT系综下的动力学模拟，温度为298 K，压力为0.1 MPa，时间步长设为1 fs。在COMPASS II 力场下进行计算，范德华力采用Atom based，静电力采用Ewald，截断半径为1.25 nm。恒温器采用Andersen，恒压器采用Berendsen，总模拟时长为500 ps。最终得到的4 种模型的密度如Tab.2所示，其中聚合度为200的BIIR模型的密度为(0.922±0.006) g/cm3，与实际BIIR 的密度0.93 g/cm3最为接近，二者之间的误差范围为0.8%，这个误差是可以接受的，因此本文将选取聚合度为200 的BIIR作为后续计算的模型。

Tab.2 Density of BIIR amorphous cell

1.3 BIIR体系对N2的等温吸附过程模拟

前文在进行密度计算时，设置温度和压力条件是常温常压(298 K，0.1 MPa)，为了模拟BIIR在40 ℃，0.1 MPa和40 ℃，1.5 MPa时的各种性质，需通过分子动力学模拟使模型达到对应的条件下的平衡状态。分子动力学模拟采用NPT系综，温度和压力分别设置为313 K(40 ℃)，0.1 MPa和313 K(40 ℃)，1.5 MPa，时间步长设为1 fs，恒温器采用Andersen，恒压器采用Berendsen，总模拟时长为500 ps。经过分子动力学模拟使BIIR无定形晶胞达到40 ℃，0.1 MPa和40 ℃，1.5 MPa的平衡状态后，为获得N2分子在BIIR中的溶解度系数，需要对N2分子在BIIR中的等温吸附进行模拟，其具体过程如下：在Sorption 模块中设置计算任务为“Adsorption isotherm”，吸附物(Sorbate)设置为经过几何优化的N2分子，吸附温度设置为313 K(40 ℃)，压力范围为0.01～5 MPa，蒙特卡洛抽样方法设置为“Metropolis”，平衡阶段和数据产出阶段的计算步数均设置为2×106步，计算得到N2分子在BIIR中吸附等温曲线。

1.4 N2在BIIR体系中的等温等压渗透过程模拟

等温等压实验条件下气体分子在聚合物中的渗透过程包括溶解和扩散两个方面。通过Materials Studio对N2分子在2种BIIR模型中的溶解过程进行模拟，其具体过程如下，在Sorption 模块中设置计算任务为“Fixed pressure”，吸附物(Sorbate)设置为经过几何优化的N2分子，吸附温度和压力分别设置为313 K(40 ℃)，0.1 MPa 和313 K(40 ℃)，1.5 MPa，蒙特卡洛抽样方法设置为“Metropolis”，平衡阶段和数据产出阶段均计算2×106步。

进一步对吸附了N2的BIIR 模型进行分子动力学模拟，研究N2在BIIR中的扩散行为。为了使系统充分平衡，保证所获数据的可靠性，本文采用的分子动力学模拟过程包括以下3 个阶段：(1)预平衡：对“Fixed pressure”模式下获得的吸附有N2的BIIR 模型的模型进行500 ps 的NPT 分子动力学模拟，温度和压力分别设置为313 K (40 ℃)，0.1 MPa 和313 K(40 ℃)，1.5 MPa，时间步长设为1 fs，恒温器采用控温更为平稳的Nose，恒压器采用Berendsen；(2)深度平衡：对上一步获得的最终构型，运行总时长为2000 ps的NPT分子动力学计算，其余模拟条件同预平衡；(3)数据采集：经过预平衡和深度平衡2步平衡后，系统达到了充分平衡，此时再运行总时长为2000 ps 的NPT分子动力学计算(模拟条件同深度平衡)，累计输出101 帧(包括初始的第1 帧)的轨迹文件，作为结果文件，用于后续分析。

2 结果与讨论

2.1 自由体积

自由体积的存在为聚合物分子链的构象调整和小分子在聚合物中的扩散提供了空间，是影响气体分子在聚合物中扩散的关键因素。Materials Studio软件中的Atom Volume&Surface 工具能够通过“硬球探针法”获得模型中的自由体积含量。BIIR 分子链中所有的原子被假定为具有范德华半径的硬球，然后将一个探针放在由硬球组成的聚合物分子表面滑动，探针所能达到的表面构成Connolly 表面，由Connolly 表面包围的体积即为聚合物的自由体积。通过上述方法，选择不同半径的探针计算得到的BIIR 模型中的自由体积分数如Fig.3 所示。探针半径越大，探针所能达到的体积就越少，因而得到的BIIR 的自由体积分数越小。对比40 ℃，0.1 MPa 和40 ℃，1.5 MPa下的BIIR自由体积分数可以发现，随着气体压力的增加，BIIR 中的自由体积分数降低。当以氮气的范德华半径0.15 nm 作为探针半径时，40 ℃，0.1 MPa和40 ℃，1.5 MPa下BIIR的自由体积分数分别为7.0%和6.7%，说明在这2种条件下N2分子在BIIR 中的活动空间仅占BIIR 总体积的7.0%和6.7%，此时BIIR 模型中的自由体积分布如Fig.4 所示。进一步对比2种条件下BIIR模型中的自由体积和占有体积可以发现，当压力提高时，BIIR中的占有体积变化不大，BIIR 晶胞的总体积和自由体积分别下降了4.7%和4.5%。这说明随着外界气体压力的提高，BIIR通过调整构象使分子排列的有序性提高，从而压缩了BIIR中自由体积的含量。

Fig.3 Free volume distribution of BIIR under different pressures

Fig.4 Simulation morphology of free volume in BIIR under(a)40 ℃，0.1 MPa and(b)40 ℃，1.5 MPa

2.2 溶解度系数

通过BIIR体系对N2的等温吸附过程进行模拟，获得不同压力下N2在BIIR 中的等温吸附曲线如Fig.5所示。随着气压的增加，BIIR中氮气的吸附量逐渐增加，且1.5 MPa 下增加幅度更缓慢，说明1.5 MPa 下BIIR 的吸附能力更弱。为了定量描述不同压力下N2在BIIR 中的吸附性质，需计算2 种状态下的N2的溶解度系数。根据亨利定律，在恒定温度下达到吸附平衡后，气体在聚合物中的浓度（c）与气体压力（p）成正比，即满足c=sp，式中：s为气体在聚合物中的溶解度系数。亨利定律只有在气体压力很低时才成立，在一般情况下，氮气在聚合物中的等温吸附过程满足Langmuir吸附等温方程，

此时p/与p呈线性关系，（为被吸附的气体在标准状态下的体积；b为吸附平衡常数；m为达到吸附饱和时被吸附气体的体积）。因此为获得N2在BIIR中的溶解度系数，可以首先以对作图并进行线性拟合(如Fig.6)，得到对应的Langmuir 吸附等温方程和参数的 m,b值，根据得到的Langmuir 吸附等温方程，当p趋近于0时，极限，

此时满足亨利定律，因此通过体积和浓度的换算即可得到N2在BIIR中的溶解度系数，

式中，V为BIIR无定形晶胞的体积。

根据拟合曲线方程计算可以得到在40 ℃，1.5 MPa 下，N2在BIIR 中的溶解度系数为6.62×10−8cm3(STP)/(cm3(BIIR)·Pa)，低于40 ℃，0.1 MPa下的溶解度系数3.05×10−7cm3(STP)/(cm3(BIIR)·Pa)。结合之前关于自由体积的计算结果可以推断，气体分子的吸附主要发生在BIIR 的自由体积内，与常压相比，高压下BIIR分子链的排列紧密程度和有序性增加，BIIR 中自由体积的含量减少，供气体分子吸附的位点数量降低，因此高压下气体分子在BIIR中的溶解能力更低，溶解度系数也更小。

2.3 BIIR体系对N2的吸附

Fig.6 Langmuir equation linear fitting graph

通过等温等压下BIIR 体系对N2的溶解吸附过程进行模拟，分析返回的吸附构型中能量最低的结构中N2的吸附情况和自由体积分布，结果如Fig.7所示。其中深色面所围成的空间代表BIIR 中的自由体积。从图中可以看出，在40 ℃，1.5 MPa 和40℃，0.1 MPa 下N2分子均主要分布在BIIR 的自由体积中。进一步分析2 种条件下N2在BIIR 中的吸附量，结果如Tab.3 所示，40 ℃，1.5 MPa 下N2在BIIR 中的平均吸附量大于40 ℃，0.1 MPa时N2在BIIR中的平均吸附量。在晶胞体积相近的情况下，40 ℃，1.5 MPa 时,1 个晶胞中最多吸附5 个N2分子，40 ℃，0.1 MPa 为2 个，Fig.8 为2 种条件下的最大吸附量时吸附的结构图。这也与实际情况下的吸附规律相同，当气体压力增加时，吸附解析平衡向吸附的方向移动，因而吸附的气体分子越多，吸附量越大。

Fig.7 Distribution of N2 and free volume in BIIR under different conditions of(a)40 ℃，0.1 MPa and(b)40 ℃，1.5 MPa

Fig.8 Simulation morphology of the maximum adsorption model of N2 in BIIR under different conditions of(a)40 ℃，0.1MPa and(b)40 ℃，1.5 MPa

Fig.9 Mean square displacement curves of N2 diffusion in BIIR under different pressures

2.4 扩散系数

对上述吸附了N2的BIIR 模型进行分子动力学模拟研究N2的扩散行为，得到均方位移曲线(Mean squared displacement)与时间的关系曲线如Fig.9 所示。纵坐标MSD表示t表示时刻(r(t))位置与最初位置(r(0))之间的位移的平方，即

式中：＜＞指系综平均。

根据Einstein方程，当气体分子进行无规行走的布朗运动时，气体的扩散系数（D）与均方位移的关系为

均方位移MSD与时间成正比。因此通过图中对均方位移曲线的直线部分进行线性拟合，所得直线(图中虚线)斜率值的1/6 可近似为N2在BIIR 中的扩散系数。根据上述的拟合计算可以得到，在40 ℃，1.5 MPa 下，N2在BIIR 中的扩散系数为1.77×10−6cm2/s，略低于40 ℃，0.1 MPa 下的扩散系数(4.57×10−6cm2/s)。结合之前关于自由体积的计算结果可以推断，N2分子在BIIR中主要通过自由体积进行扩散，与常压相比，高压下BIIR 分子链的排列紧密程度和有序性增加，BIIR中自由体积的含量减少，气体分子的扩散变得更困难，因此高压下，N2分子在BIIR中的扩散系数低。

2.5 渗透系数

Tab.3 Adsorption capacity of N2 in BIIR under different pressures

渗透系数（P）、溶解度系数（S）(和扩散系数（D）的关系为P=D·S，Tab.4 列出了不同条件下N2在BIIR 中的溶解度系数、扩散系数、渗透系数。通过对比可以发现，40 ℃，1.5 MPa 下的溶解度系数、扩散系数、渗透系数均比40℃、0.1 MPa 条件下低。这是因为，N2在BIIR中渗透的溶解过程和扩散过程均发生在BIIR的自由体积中，而气体压力的增加导致BIIR的自由体积分数降低，从而导致一方面N2在BIIR中的吸附位点减少，溶解度系数降低；另一方面使N2的扩散空间受到限制，扩散系数降低。在溶解、扩散两方面因素的综合影响下，40 ℃，1.5 MPa下N2在BIIR 中的渗透系数比40 ℃，0.1 MPa 条件下降低约92%。

为验证该结论，本团队定制了1 台高压气体渗透仪，测试了溴化丁基橡胶样品在不同压力下的渗透系数，结果如Tab.4所示，BIIR的渗透系数从40 ℃，0.1 MPa的2.87×10-13cm3(STP)·cm/(cm2·Pa·s)下降到40 ℃，1.5 MPa的2.35×10-14cm3(STP)·cm/(cm2·Pa·s)，下降幅度约为95%，与模拟值相近，证实了通过分子模拟法判断材料气体阻隔性能的可行性。

Tab.4 Solubility coefficient,diffusion coefficient and permeability coefficient of N2 at different pressures in BIIR

2.6 扩散机理

通过对N2在BIIR体系中扩散的轨迹进行模拟，得到不同条件下N2在BIIR 中的位移-时间曲线(Fig.10)和三维运动轨迹图(Fig.11)。由Fig.10 中N2分子的位移随时间的变化关系可以明显看出，N2分子在BIIR 中有2 种运动状态：1)N2分子的位移在一段时间内仅发生小范围波动，波动范围仅0.1～0.5 nm，这种状态可以理解为N2分子在某一较小空间内的振动；2)N2分子的位移在约20 ps的短时间内实现在2 种振动状态间的跳跃，每次跳跃的距离约为1 nm，2 次跳跃的时间间隔约500 ps，这种状态可以理解为气体分子在不同空间之间的跃迁。从Fig.11的N2分子的运动轨迹图也能更直观地看出，N2分子在扩散过程中交替进行着小范围的振动和长范围的跃迁。因此，通过以上分析可以得到，N2在BIIR 中的扩散符合“空穴跳跃扩散理论”：N2分子在BIIR中扩散时，大部分时间被限制于聚合物分子链间的空穴中，气体分子在空穴中的来回振动对扩散的贡献较小，当BIIR的链段运动导致分子链的构型调整使相邻的空穴之间形成瞬时通道时，气体小分子将通过瞬时通道从一个空穴跃迁至另一个空穴，从而完成扩散。

Fig.10 Displacement of N2 in BIIR under(a)40 ℃，0.1 MPa and(b)40 ℃，1.5 MPa as function of simulation time

Fig.11 Diffusion trajectory of N2 in BIIR under(a)40 ℃，0.1 MPa and(b)40 ℃，1.5 MPa

3 结论

通过Materials Studio软件分别对40 ℃，0.1 MPa和40 ℃，1.5 MPa下N2分子在BIIR中的渗透行为进行模拟，研究压力对N2在BIIR 中的溶解度系数、扩散系数、渗透系数的影响及其扩散行为，结果表明：

(1) 40 ℃，0.1 MPa 和40 ℃，1.5 MPa 下BIIR 的自由体积分数分别为7.0%和6.7%，说明气体压力的增高使BIIR分子的链排列更紧密，BIIR的自由体积被压缩。

(2)分析N2分子在BIIR中的吸附图和吸附量发现，N2分子在BIIR中的吸附主要发生在BIIR的自由体积中，由于气体压力的增高使吸附脱吸平衡向吸附的方向移动，因此40 ℃，1.5 MPa 下N2在BIIR 中的吸附量比在40 ℃，0.1 MPa下高。

(3)自由体积从溶解和扩散两方面影响N2分子在BIIR 中的渗透。由于高压下BIIR 的自由体积的降低，N2分子在BIIR中的吸附位点减少，因而40 ℃，1.5 MPa下N2分子在BIIR中的溶解度系数比40 ℃，0.1 MPa下低；高压下自由体积的降低压缩了N2分子在变化中的扩散空间，从而使40 ℃，1.5 MPa下N2分子在BIIR 的扩散系数比40 ℃，0.1 MPa 低。综合溶解和扩散两方面，40 ℃，1.5 MPa下N2分子在BIIR中的渗透系数比40 ℃，0.1 MPa下低，试验结果也证实了这一变化趋势，表明分子模拟法可用于判断材料的气体阻隔性能。

(4)分析N2在BIIR 中的位移发现：N2在BIIR 的扩散过程满足“空穴跳跃扩散理论”，即N2分子在BIIR中扩散时大部分时间在BIIR中的空穴中振动，当分子链结构的调整使相邻的空穴之间形成瞬时通道后，N2分子通过通道从一个空穴向另一个空穴跃迁，从而实现扩散。