交联对硅橡胶热导率影响的分子动力学模拟*

徐文雪梁新刚† 徐向华祝渊

1)(清华大学航天航空学院,热科学与动力工程教育部重点实验室,北京100084)

2)(南方科技大学深港微电子学院,深圳518055)

(2020年5月17日收到;2020年6月15日收到修改稿)

1 引 言

电子器件功率密度的增加使得如何有效散热成为一个关键问题[1].同时,电子元器件若不能做到有效散热,必然会导致温度升高,降低可靠性.热失效是电子设备最主要的失效模式,据相关数据统计[2],超过半数的电子设备失效是由不良的冷却导致的.因此电子元器件的降温散热十分重要.

电子器件的散热过程中,利用界面材料降低界面热阻是一个重要环节.硅橡胶具有绝缘、耐热等优势,是一种常用的热界面材料基体.硅橡胶基体通常和高热导率材料粒子结合形成复合材料,其拥有许多产品形态[3].影响硅橡胶复合材料导热性能的因素主要有硅橡胶基体、导热填料以及加工工艺等.影响硅橡胶基体导热性能的主要因素有取向度、交联度、分子内部结构的致密性和温度等.导热填料及加工工艺对硅橡胶热导率的影响有大量的实验研究[4,5],但缺少硅橡胶基体结构对硅橡胶导热影响研究.硅橡胶是由高黏滞性的塑性态硅橡胶生胶经过硫化反应得到的,硫化的目的是通过加成、缩合或者氧化等化学反应过程形成交联结构.因此本文研究交联对硅橡胶基体热导率的影响.

有大量关于无定形材料内部交联形成空间网状结构来提升聚合物热导率的研究.Kanari和Ozawa[6]报道了不同交联密度下环氧树脂的热导率,研究表明热导率随着交联密度的增大线性增大.Kikugawa等[7]用分子动力学方法研究了交联对聚乙烯和聚苯乙烯热导率的影响,结果表明热导率会随着交联密度增大而增大,但是两种材料提升的效果不同;此外他们还发现在低分子量形成的聚乙烯中,交联能对热导率有更大的提升.Li等[8]报道了交联官能团密度对有机硅基材料热导率的影响,结果表明交联密度对热导率有提升的效果,文章认为高交联度的复合材料能提供更多传播途径且自由体积更小,所以热导率会更高.目前尚未见硅橡胶交联密度对热导率影响的报道,因此开展硅橡胶交联密度对热导率的研究是有必要的.

2 计算方法

本文采用基于傅里叶导热定律的非平衡分子动力学方法计算热导率—双温度方法[9](dualthermostat method):该方法需要建立一个温度分布,统计维持稳定的温度梯度需要的热流.系统沿热流方向被分为若干层,设置两端作为热端区域,中间层作为冷端区域,分别控制冷/热区域的温度,并统计控温过程中的能量输入与输出,能量的统计是通过粒子的动能变化来计算的,热流方向为图1中箭头的方向.热导率的计算根据以下公式进行[9]:

其中DEtotal表示时间间隔Dt内总能量变化,DEH表示热端区域的动能变化量,DEC表示冷端区域的动能变化量,DEerr表示计算误差导致的能量变化;(2)式中,l是热导率,下标i取值为H或C,分别表示热端和冷端,DEi分别按区域求解即可得到热端/冷端区域的动能变化;A是横截面积,dT(z)/dz表示沿热流z方向的温度梯度.计算中硅橡胶采用COMPASS力场[10],控温方法采用Nosé-Hoover方法[11],非键截断半径与库仑截断半径均为12.5Å,时间步长为0.05 fs[12].

图1双温度方法示意图[13]Fig.1.Schematic diagram of dual-thermostat method[13].

3 交联硅橡胶的建模

甲基乙烯基硅橡胶和交联剂DBPMH的反应[14]是一种常见的交联反应,两者的分子结构如图2所示,图中R表示端基,可以是甲基或者乙烯基,DBPMH含有两个酯键,通过断开酯键与硅橡胶高分子链形成交联键,从而得到交联硅橡胶,在图3中用ROOR表示.反应过程如图3所示,在交联剂的作用下,硅橡胶产生两种交联反应:一种是硅橡胶两个甲基中的碳原子形成交联键,如图3(a)所示;另一种是双键中的碳原子与甲基中的碳原子形成交联键,如图3(b)所示.

图2甲基乙烯基硅橡胶和交联剂DBPMH的分子式Fig.2.Molecular formula of methyl vinyl silicone rubber and crosslinker DBPMH.

图3甲基乙烯基硅橡胶和DBPMH的反应(a)硅橡胶两个甲基中的碳原子形成交联键;(b)双键中的碳原子与甲基中的碳原子形成交联键Fig.3.Reaction of methyl vinyl silicone rubber with DBPMH.

COMPASS力场[15]是收费力场,公开的COMPASS力场缺少C=C双键的力场参数,且甲基乙烯基硅橡胶含有极少量的C=C双键,通常低于1%;再者,本文主要研究内容是交联形成的结构对硅橡胶的影响.综上所述,采用图3(a)所示交联反应过程,忽略乙烯基形成的交联结构和甲基形成的交联结构的区别,针对该反应本文做出如下简化假设:

1)与硅氧链高分子相比,交联剂的分子相对较小,交联剂活动性远远大于硅氧链高分子,即交联反应不会由于交联剂的缺少而进行不下去;

2)交联分子非常容易形成自由基,建模过程中直接设置活性原子形成;

3)各个自由基反应活性相同,反应是扩散控制的,即距离决定是否反应.

图4交联反应的流程Fig.4.Process of crosslinking reaction.

图5 设定的反应活性原子R1和R2 Fig.5.Set reactive active atoms R1 and R2.

为了让交联过程较为合理,采用交联-松弛交替运行的策略模拟交联过程,该流程要求两个反应活性原子的间隔不能相距太远,交联反应的流程如图4所示.交联第一步需要构建含有反应活性原子的无定形硅橡胶模型.不同链上的反应活性原子分别被命名为R1或R2,为了区别非活性原子,R1和R2被加大显示,如图5所示.在Materials Studio软件的Amorphous模块中封装含有反应活性原子的高分子链得到计算模型,接下来在Forcite Plus模块中采用几何优化和退火对模型进行初步的能量最小化,充分弛豫后得到初始结构进入交联反应.设置计算的基本参数,统计最大可能交联键数目和已经形成的交联键数目,并判断需要交联的键数目是否小于最大可能交联键数目,不满足直接退出,满足则进入下一步.对模型进行初步优化,在NVT,NPT系综下松弛.判断当前搜索半径是否小于设定的最大搜索半径,不满足则直接退出交联反应,满足则进入下一步优化平衡,并执行建立新键命令.建立新键命令是以每一个活性原子为原点,判断是否有未反应的活性原子在搜索半径之内,距离小于搜索半径则建立两原子的近邻接触,判断两个原子是否均为未成键的原子,如果是的话则施加一个逐渐加大的简谐约束力,减少近邻接触原子间的距离,使原子对的距离缩短至成键距离,建立新键;否则跳过该近邻接触判断下一个近邻接触,若未找到满足条件的近邻接触,则跳出该循环,增大搜索半径进入下一个循环.建立新键以后对模型进行几何优化并施加微扰,统计模型的热力学性质,判断当前交联密度enow和目标交联密度etarget的大小:如果enow

图6建立新键的过程示意图Fig.6.Process diagram of creating new bond.

关于交联密度的定义目前不是很统一.实验研究方面由于无法得到内部分子链的具体结构,所以只能通过静态模量、平均分子质量等宏观参数来表示交联密度.基于不同的模型或者简化理论,存在不同的交联密度表达式:Popovic和王韶晖[16]通过静态模量方程计算交联密度g,表达式为g =r/Ma,r是聚合物的密度,Ma是两交联键之间的数均分子质量;胡冬红等[17]用Mooney-Riving方程求解高分子链上两个交联点间的有效分子量Mc,用单位质量高分子链中交联键的数目表示交联密度D,D=1/(2Mc).此外不同的材料、不同掺杂时,颗粒大小、颗粒表面性质等都会影响交联的程度,所以不同材料间的交联密度难以比较.数值模拟中交联密度是用来表征分子交联的程度,表征相对意义,所以本文交联密度的定义是已经反应的交联键数目与最大可能交联键数目的比值,与Kikugawa等[7]的定义本质相同.

4 计算与讨论

4.1 交联键位置对硅橡胶热导率的影响

图7是按照上述交联流程,起始搜索半径设置为3.5Å,终止搜索半径设置为8.0Å,目标交联密度是100%的时候得到的若干个交联结构的示意图,绿色加粗的原子对即为交联键上的原子.每个计算盒子内保证105个反应活性点R1和105个反应活性点R2,所以每个计算模型中添加了图5所示的长链共30条, 模型的尺寸在32Å ×32 Å ×70Å左右.由于建模中会用到NPT和NVT系综进行松弛,所以每个模型的密度和尺寸会略有不同.

图7 不同交联密度下的结构示意图Fig.7.Structure diagram under different crosslinking density.

图8冷热端能量随模拟时间的变化Fig.8.Energy variation at hot and cold ends.

根据上述交联方法得到了不同交联密度的硅橡胶结构.按照双温度模拟方法,计算得到了不同交联密度下硅橡胶的热导率,图8是交联密度为0.95%时能量随时间的变化趋势图,可以看出散热已进入稳定阶段,其他交联密度下结果类似.图9是交联密度对热导率的影响,可以看出,随着交联密度的变大,热导率逐渐升高,交联密度达到80%的时候,热导率提高了40%.由于硅橡胶密度变大会导致单位体积内原子数量变多,所以激发的声子模式也会变多,因此单位体积比热也会变大,而这对热导率是有利的.图10是不同交联密度下结构的质量密度,可以看出,交联密度为80%的时候,质量密度只提高了5%.如果把硅橡胶密度对热导率的影响近似看成是线性的,把密度的倍数乘以交联密度为1%时的热导率得到一个由质量密度增大导致的热导率的大小,将该值与模拟值对比,结果如图11所示.可以看出,由质量密度增大导致的热导率增大远远不及计算得到的热导率.所以可以推测,交联形成的空间网状结构导致传热模式的改变,进而导致热导率有一个较大的提升.这进一步验证了Li等[8]认为提高有机硅基复合材料的交联度能提供更多传播途径且自由体积更小,所以对热导率有提升效果的结论.

图9交联密度对热导率的影响Fig.9.Influence of crosslinking density on thermal conductivity.

图11由密度导致的热导率增大与模拟对比Fig.11.Comparison between the increase of thermal conductivity caused by density and simulation results.

4.2 交联键位置对硅橡胶热导率的影响

上一节中设置的反应活性点的位置是均匀地散布在链长上的,如图12(a)所示.为了研究反应活性点的位置对热导率是否有影响,本节对只在端部设置反应活性点和只在中间设置反应活性点的两种情况进行比较.图12中每幅图下面的红色点代表设置了反应活性点的位置,蓝色点代表未设置反应活性点.每个计算盒子内保证有100个反应活性点R1和100个反应活性点R2,所以每个计算模型中共添加了50条硅氧链.端部设置反应活性点的分子链形(图12(b))在交联反应中形成的交联密度分别有:12%,19%,27%,38%,51%,61%和71%;中间位置设置反应活性点的分子链形(图12(c))在交联反应中形成的交联密度分别有:1%,9%,20%,31%,41%,52%和63%.反应活性点在不同位置时交联率对热导率的影响结果如图13所示,可以看出,交联点的位置不论在端部还是中部,热导率都会随着交联密度的增大而增大,端部和中部设置交联点对热导率的影响差异较小.由此可推得反应活性键的位置对热导率影响较小,选择不同位置上的交联剂对热导率影响较小.

图12反应活性原子不同位置分布(a)均匀;(b)端部;(c)中间Fig.12.Positions distribution of reactive atoms.

图13反应活性点在不同位置时交联密度对热导率的影响Fig.13.Effect of crosslinking density on thermal conductivity at different positions of reactive atoms.

4.3 反应活性点间隔对硅橡胶热导率的影响

前面两节的建模中,每隔2个硅原子设置一个活性反应点.为了验证反应活性点的间隔对热导率是否有影响,本节设计了反应活性点间隔稍大的两种分子链结构,即每隔4个和6个Si原子设置一个反应活性点,如图14所示.间隔为4个Si原子的分子链上有4个活性反应点,间隔为6个Si原子的分子链上有3个活性反应点.为了保证反应活性点的数目基本一致,间隔为4的计算盒子内封装了50条高分子链,保证有100(4×25=100)个反应活性点R1和100个反应活性点R2,模型的尺寸在33Å×33Å×118Å左右;间隔为6的计算盒子内封装了66条高分子链, 保证有99(3 ×33 =99)个反应活性点R1和99个反应活性点R2,模型的尺寸在37Å×37Å×117Å左右.间隔为4时,模拟交联反应形成的交联密度分别有:16%,31%,49%,65%;间隔为6时,模拟交联反应形成的交联密度分别有:3%,20%,39%,55%,72%.

图14不同反应活性点间隔的分子链(a)间隔为4;(b)间隔为6Fig.14.Molecular chains with different reactive atoms intervals.

由于不同间隔下系统的体积不一样,单位体积内的交联点数不一样.间隔为2的时候,活性点的设置最密集,间隔为4和6时,活性点的设置则稀疏一些.为了保证交联密度的可比性,采用表1所列的方式对交联密度进行归一化处理.

表1归一化交联密度的计算Table 1.Calculation of normalized crosslinking density.

不同反应活性点间隔对热导率的影响如图15所示,间隔为2和4的情况下,热导率随着归一化交联密度的提高而提高.间隔为6的情况下,当归一化交联密度小于20%时,热导率不随归一化交联密度变化;归一化交联密度大于20%时,热导率增加.对比间隔2,4,6时的热导率,发现相同归一化交联密度下间隔为6的热导率最大.交联间隔越大,热导率越高,为了提高交联硅橡胶的热导率,应当适当增加活性点的间距.

图15不同反应活性点间隔对热导率的影响Fig.15.The effect of reaction active atoms intervals on thermal conductivity.

不同的交联间隔下的密度有所不同,密度随归一化交联密度的变化趋势如图16所示,对比热导率随着归一化交联密度的变化趋势可以看出,密度的变化趋势大体与热导率的变化趋势一致,通常是较高的密度对应较高的热导率.把质量密度和热导率放在一张图中,如图17所示,可以看出热导率和质量密度正相关,相关系数为0.838.实际应用中应将硅橡胶材料紧密加载,这是因为紧密加载会使其密度增高进而使其热导率增加,此外紧密加载时接触热阻也会较小.

图16不同归一化交联密度下的质量密度Fig.16.Density under different normalized crosslinking density..

图17 质量密度和热导率的相关性Fig.17.Correlation between mass density and thermal conductivity.

4.4 交联密度对声子态密度的影响

晶格振动的量子称为声子,通常用分子动力学模拟的信息得到声子的相关信息[18],最常见的是由速度自相关函数的傅里叶转换得到声子态密度[19].还有一种方法是基于振动散射理论得到声子的信息[20],它是利用分子动力学计算中的位移进行傅里叶变换得到倒空间的格林函数[21],再进一步得到力常数矩阵和动力学矩阵,通过计算动力学矩阵的特征值得到声子态密度[22],该方法能考虑到声子的非谐波效应.本文声子态密度的计算方法采用后者, 这一部分工作在LAMMPS软件的PHONON命令中完成,结合后处理工具PHANA可以得到声子态密度的信息.

图18间隔为2时不同交联密度下的声子态密度分布(a)Si原子;(b)O原子Fig.18.Distribution of phonon density of states under different crosslinking densities at atoms interval 2:(a)Si;(b)O..

图19间隔为4时不同交联密度下的声子态密度分布(a)Si原子;(b)O原子Fig.19.Distribution of phonon density of states under different crosslinking densities at atoms interval 4(a)Si;(b)O.

图18所示为间隔为2的情况下交联密度分别为1%,10%,75%,80%时Si和O原子振动的声子态密度;图19所示为间隔为4的情况下交联密度为16%和65%时的声子态密度.可以发现,随着交联密度增加低频区声子数变少,在高频区声子数增加.交联密度增加,使得原子间的相互作用力加强,高频声子增加;但是另外一方面,交联密度的增加,使得从低温区向高温区热量传递的路径更短,因而使得热导率随交联密度的增加而增加.

5 结 论

本文根据甲基乙烯基硅橡胶的交联反应原型作出三条简化假设,并根据交联反应的流程,建立了不同交联密度下的模型,计算了热导率.研究发现,热导率随交联密度的提高而提高,交联密度达到80%的时候,热导率提高了40%;通过对比分析发现,由结构的质量密度导致的热导率增大远远不及计算得到的热导率,因此可以认为交联形成的空间结构网络对热导率的促进作用更大.对不同交联键位置的硅橡胶热导率的比较表明,其热导率都随交联密度变大而变大.当设置的总的最大交联数相同时,端部设置反应活性点和中间位置设置反应活性点对热导率基本没有影响,选择不同位置上的交联剂对热导率影响较小.不同交联点间隔对热导率影响的研究表明,间隔较小的情况下,热导率随着交联密度的提高较明显;间隔较大的情况下,热导率随交联密度先不变后增加.交联间隔越大,热导率越高,为了提高交联硅橡胶的热导率,应当适当增加活性点的间距. 分析不同间隔下的质量密度与热导率发现两者正相关,相关系数为0.838.对声子态密度进行了计算,发现交联密度变大,低频声子占比变少,高频声子占比变多,这是因为交联使原子间相互作用增强,所以频率变高.热导率随交联密度变大而变大的原因是交联键的加入缩短了热传导的路径.