陶瓷轴承Si3N4-GCr15摩擦副的润滑[1]及冷却问题一直是制造领域关注的热点。石墨烯由规则六边形经sp2杂化的碳原子构成，其特殊结构使其具有极其优良的性质，具有高的面积体积比，良好的导电性导热性，单层石墨烯导热率可达5 300 W/(m·K)[2],已超越碳纳米管、石墨，金刚石等碳同素异形体的极限，更远超银和铜等金属材料，在强化传热领域优势显著。尝试将石墨烯作为添加剂制备石墨烯润滑油，用于Si3N4-GCr15摩擦副润滑同时，充分利用石墨烯导热特性增强润滑油在摩擦区域的导热特性具有重要意义。石墨烯在强化传热冷却方面的研究取得了较大进展。从最初的将石墨烯简单地与聚合物共混来提高导热性，到制备石墨烯纳米流体强化对流传热，通过表面修饰或材料复合强化相变传热。Ota J等[3]应用机械化学方法将石墨烯作为介质添加到烃类油中，石墨烯在油中的分散提高了基础油的导热系数，并且表现出了良好的散热性。
对石墨烯具体形态的研究主要通过分子动力学模拟，王建等[4]采用非平衡分子动力学模拟(NEMD)方法研究了石墨烯/聚乙烯纳米复合材料的界面导热性能,主要考察了石墨烯层数、尺寸对界面热阻的影响。实验研究结果表明:当石墨烯层数为单层时,界面热导率为46.79 W/(m·K),随着石墨烯层数的增加,会使热导率下降;但石墨烯层数超过四层后,界面热导趋于恒定接近39.00 W/(m·K)。惠治鑫等[5]通过模拟发现，硅原子的加入改变了石墨烯声子的模式、平均自由程和移动速度，使得单层硅功能化石墨烯模型的热导率随着硅原子数目的增加而急剧地减小。在300～1 000 K 温度变化范围内，单层硅功能化石墨烯的热导率随温度升高呈下降趋势，具有明显的温度效应。Hu等[6]利用分子动力学模拟了对称石墨烯纳米带与不对称石墨烯纳米带的热导率。发现相同尺寸的对称石墨烯纳米带的热导率大于不对称石墨烯纳米带的热导率，原因是在不对称石墨烯纳米带中存在热整流。Shin Pon Ju等[7]研究了石墨烯和折 叠石墨烯的导热率。单层石墨烯的导热是2 212 W/(m·K)，折叠后的导热率约是71.4 W/(m·K)，这说明通过折叠会大幅度降低导热率。Chenhan Liu等[8]研究了多层的夹层界面的热阻，随着界面压力的增加，热阻降低，当压力达到10 GPa时，热阻降低了4倍，压力施加到内边界时，热阻稍微上升，两个方法都施加时热阻的变化和外边界一样。
目前石墨烯模拟导热方面研究较多，但将石墨烯用在某工作场合进行导热的研究较少。石墨烯在强化传热方面的作用原理复杂，状态和机制多变，尚没有完善的传热理论模型，进一步明确纳米尺度强化传热的机理，建立完善的理论模型，对指导石墨烯强化传热具有重要意义。
本文通过建立分子动力学模型，从微观尺度分析层数，堆叠方式、环境温度及环境压力对石墨烯导热特性的影响。同时，基于实验，从宏观尺度分析石墨烯对Si3N4-GCr15摩擦副温度场的影响。
1 石墨烯导热特性分子动力学模拟
1.1 非平衡分子动力学模拟
分子动力学是为求解体系中分子或原子间作用势和系统外加约束共同作用下，随时间推进的微观过程，最后通过统计方法得到系统宏观平衡参数或输运参数。为了计算具有温度梯度的非平衡系统，本文采用非平衡态分子动力学方法对图1所示动力学模型进行热导率分子动力学模拟[9]。模拟结果可被用来检验理论的正确性，又可以与实验值相比较，以检验和改进模型，并提供更多的数据。许多在实验中无法实现的情况，都可以通过模拟来观察。热导率模拟原理是一种速度交换算法，可以使得加入热流时不改变系统的动量。这种方法速度的变化公式：
vi,new=vG+α(vi,old-vG)
(1)
温度和速度的关系满足玻尔兹曼统计规律：
(2)
其中，KB代表玻尔兹曼常数。控温系数：
(3)
控温部分原子的质心速度：
(4)
所以：
(5)
动量守恒，这时使用 “+” 时施加于热端。使用 “—” 时施加于冷端，相对变化能量定义为：
(6)
通过这个算法，相当于在热端加入从冷端提取的能量。所以热流的表达式为：
(7)
然后利用傅里叶导热定律来计算热导率。
(8)
1.2 模拟过程
模拟过程本文使用MS (materials studio)来构建仿真模型。首先根据石墨稀的几何结构，利用MS软件构建石墨稀的单胞模型,，然后构建成石墨稀片层。接下来构建GCr15晶胞，利用层建模将石墨烯与GCr15构建成Graphene-GCr15晶胞，如图1所示，构建出不同石墨烯层数晶胞，不同石墨烯堆叠方式，不同石墨稀形态的晶胞。
图1 单层Graphene-Gr15模型示意图
本文应用MS软件forcite模块的Universal力场来描述粒子间相互作用，选择该力场就是选择了一种原子参数定义方式和一种势函数表达式，包括非键作用的范德华力，静电作用以及化学键的作用等。在MS软件中Quality为模拟的截断半径，设为1.25 nm，范德华力计算用截断方法，静电力的计算用Ewald方法。
在进行分子动力学模拟时需调整模型内部分子的位置坐标优化体系结构[10]。首先采用Verlet速度下降法对模型粗略的优化，然后应用共轭梯度法进行精细的结构优化，从而寻找到最佳动力学模拟构型。为了获得稳定准确的动力学结果，先将模型在NPT系综下进行 100 000 步的动力学模拟，利用控制温度的方法使得体系达到一个能量稳定的状态，温度控制选Andersne控温方法，积分的时间步长0.1×10-3 ps，积分时间为100 ps，压力控制选Berendsen控压方法，积分的时间步长0.1×10-3 ps，积分时间为100 ps，热导率的计算需要准确的热流密度值，因此保持体系内部总能量不改变情况下需要在NVE系综下进行100 000步动力学模拟，参数设定与在NPT系综步骤相同。
2 石墨烯导热特性分子动力学模拟结果
在非平衡态分子动力学模拟过程中,判断模拟结果可靠性的一个重要原则就是沿着热流方向建立的温度梯度是否呈线性分布。为了说明模拟结果的正确性，图2给出了体系在不同时间段建立的温度分布。从图中可以看出，在 3个不同时间段统计的模拟温度曲线基本一致，说明本文中分子动力学模拟的时间是足够的，模拟结构经过100 ps 后可以达到平衡状态。热浴区域温度跳动较大，这是由于导热声子在热域发生了强烈的散射，导致温度曲线出现了非线性。通过模拟结构的温度曲线获取温度梯度时，本文只考虑了模拟结构中间的线性区域，所以使用 Fourier定律计算热导率是可行的。通过分子动力学得到不同石墨烯形态下的分子动力学模拟结果。
图2 不同模拟时间各区的温度分布曲线
2.1 石墨烯层数对导热率的影响
石墨烯的层数[11]会影响石墨烯的导热率，不同层数石墨烯的Graphene-GCr15模型分子动力学模拟计算结果如图3所示，不同层数石墨烯模拟出导热率结果不同，从图中看出层数越多最终平衡时温度越高。
图3 不同层数Graphene-Gr15模型模拟计算结果
图4表示不同层数情况下计算的导热率，层数越少导热率跳动越大，单层情况下导热率跳动较大，单层石墨烯确实具有较高导热率，随着石墨烯层数增多导热率下降较为明显，图中单层石墨烯情况下导热率比双层石墨烯高近40%，比3层情况下高近60%。
图4 298 K下不同层数的模型模拟计算的导热率曲线
2.2 石墨烯堆叠方式对导热率的影响
3层石墨烯3种不同堆叠方式如图5，石墨烯片层的位置变化导致导致声子耦合变弱，声子色散能力变弱[12]。3种情况的Graphene-GCr15模型导热率计算值如图6所示，从曲线可以看出AAA堆叠方式导热率最高，ABA方式导热率下降了很多，相同模拟时间下导热率下降近40%，ABC方式下导热率下降更多。
图5 3种堆叠方式示意图
图6 3种堆叠方式的导热率计算值曲线
2.3 温度对导热率的影响
图7是单层石墨烯的Graphene-GCr15模型在不同温度下计算的导热率。从图中看出初始时传热过程较为缓慢，导热率较低，因为在初始弛豫时热域冷域温度梯度较大，初始不能保持较高导热率[13]，随着计算时间的延长，在弛豫过程中导热率会逐渐回升并趋于平稳。在温度较低时，导热率变化不明显，在温度为448 K和498 K时，导热率变化较大。结果表明温度对石墨烯导热率存在影响，高温情况下石墨烯会有较高导热率，控温498 K时导热率比低温情况导热率高。图8是不同层数石墨烯Graphene-GCr15模型在不同温度情况下的模拟导热率曲线，图中随着温度升高，单层石墨烯导热率升高最为明显，双层石墨烯情况下导热率渐趋于平稳，3层情况导热率基本不变。层数越多对声子传输的影响越大，进而会严重影响传热过程。
图7 不同温度下计算的导热率曲线
图8 不同层数石墨烯Graphene-GCr15模型在不同温度下的模拟导热率曲线
2.4 压力对导热率的影响
图9是单层石墨烯Graphene-GCr15模型在不同环境压力下的模拟曲线，GCr15与Si3N4陶瓷摩擦的工况为0.2～0.4 MPa[14]，图中曲线说明在单层Graphene-GCr15模型在不同环境压力下表现出不同的导热情况，高的环境压力会表现出高的导热率，在0.4 MPa时的导热率比低的环境压力下导热率要高，在达到一定模拟时间下导热率渐趋于平稳，不同压力下导热率情况不同，在压力为0.4 MPa时导热率最高。,在气压较高情况下，冷却效果会更好，但在压力过高时，导热率不会一直升高。
图10是描述了不同层数石墨烯-GCr15模型随着压力变化的模拟导热率曲线，图中看出3种情况导热率趋势表现不同，随着压力升高导热率随着升高,但层数越多导热率变化不明显甚至下降，层数达到3层时，随着压力升高导热率不在变化。单层石墨烯的高导热率在随着压力提高表现的更加明显，从已知文献[15]得知，由于层数增加在高压力下石墨烯极容易变形及堆积折叠等，会对导热过程产生反作用。
图9 不同压力下计算的导热率曲线
图10 不同层数石墨烯-GCr15模型随着压力变化的模拟导热率曲线
3 石墨烯在Si3N4-GCr15摩擦副中导热特性实验分析
3.1 实验方案
实验采用美国Rtec多功能摩擦磨损试验机(以下简称试验机)实现旋转及往复运动以及Spindle Error Analyzer Version 8.6的TMP190模块测试温升。在图11(a)旋转实验过程中采用FLIR红外摄像仪实时录像，在图11(b)往复实验中由于有加热温箱，利用测温传感器监测温度。测温传感器1测量摩擦副位置1即Z方向的温度，测温传感器2测量摩擦副位置2即X方向温度。
图11 旋转摩擦实验(a)和往复摩擦实验(b)示意图
实验选用GCr15作为试件材料，试件为直径50 mm，高5 mm的扁平型柱体。陶瓷球材质为Si3N4陶瓷材料，其直径为9.5 mm，实验所用中国青岛德通纳米科技有限公司生产石墨烯参数如表1。利用试验机旋转模块，将试件装夹在卡盘上，在试件旋转的同时保持陶瓷球在夹具上不动并与试件接触产生滑动摩擦。实验时利用机械方法将石墨烯压在轴承钢表面。GCr15轴承钢试件及Si3N4陶瓷球，在陶瓷球与试件接触处上方10 mm放置测温探头用以记录实验过程中的温升情况。假设此处温度与摩擦处温度一致。在上方放置红外摄像头，观测摩擦过程中温度变化，考虑到实验环境相同并无不同干扰因素，故实验过程中的热量散失问题及对流换热忽略不计[16]。
表1 石墨烯参数
项目数据成份100%石墨烯面密度/(mg·m-2)0.77形状蜂窝状比表面积/(m2·g-1)887.6状态粉末孔容/(cm3·g-1)0.22
本实验为了获得摩擦副间工作过程温度变化，实验依据牛顿冷却定律，即
Q=A(Tw-Tf)
(8)
式中，Q为对流传热的热流密度(J/(m2·s))；A为换热系数；Tw为试件温度，Tf为流体温度(K)。试件温度与流体温度由实验过程中温升测试仪记录数据，热通量由摩擦力计算得出，即
(9)
式中，Fx为摩擦力(N)；D为摩擦位移(m)；S为试件对流换热表面积(m2)；t为对流换热时间(s)。Fx由Rtec试验机测得，D由圆形磨痕周长与转数求得，实验时间10 min。根据以上理论依据设计实验方案状态参数如表2。在旋转试验台完成有无石墨烯的对比试验，对比摩擦区温度变化；在往复试验台完成温度对石墨烯导热的影响，选取五组温度梯度，由于存在温箱，红外照相不适用，选用传感器测温。
表2 实验方案状态参数
序号载荷/N石墨烯质量分数/%实验温度/ ℃时间/min试验台16002510旋转2600.52510旋转3600.525,55,85,115,14510往复
3.2 实验结果
实验得到在旋转试验台上两种情况下的红外影像，及摩擦区域温度变化的3D图像。图12是在不存在石墨烯情况下的红外摄像，中心区域温度较高，温度较于集中在300～400 mm处，0～300 mm及400～600 mm之间不存在高温，磨痕以外区域温度相对较低，说明热量难以向四周传导，散热情况不好。图13是在摩擦区存在石墨烯的情况，图中可以看出中心区域温度仍较高，相比图12此时四周温度较高，在有石墨烯的区域能表现出较高温度，说明此时导热状态较之前良好。磨痕以外区域具有较高温度形成了分明的温度梯度。
图12 无石墨烯情况的红外影像
图13 有石墨烯情况的红外影像
图14是在往复试验台不同实验温度下的摩擦区温度曲线，在不同温度情况下，温升趋势不同。从图中看出X方向温度变化较Z方向剧烈，同温度下温升较为明显。在25 ℃及55 ℃情况下，温度变化梯度较小，随着温度升高，曲线变化明显，在实验温度145 ℃时，温度梯度近30 ℃，说明在高温情况下摩擦区热量向四周散发较快，导热效果较好。所以控制一定的环境温度范围对发挥石墨烯导热作用具有促进作用，使得摩擦副保持在合理温度区间内，避免温差较大发生严重变形，提高电主轴精度。
图14 温度曲线
图15是每次实验的摩擦因数曲线，可以看出，摩擦因数在不同情况下表现不同，在145 ℃时摩擦因数较低，说明摩擦情况良好，导热情况比低温时好，所以具有一定温度条件对石墨烯导热具有促进作用。
图15 摩擦系数曲线
对于实验可能产生的误差，主要有两点原因：1)实验温度没有严格控制，实验过程中冷却不充分，第二次实验时陶瓷球及试件仍具有一点初始温度，实验初始温度很难控制。2)受限于实验仪器精度，红外摄像仪长时间使用自身机体产生热量会有一定初始温升，导致后面实验具有一定误差。
4 结论
1) 石墨烯在摩擦过程中能够起到强化导热作用，散热效果好。
2) 单层石墨烯情况下导热率最高，层数增多导热率下降，四层或五层石墨烯下降更加明显。3层情况下AAA堆叠对石墨烯导热性能最高。
3) 高温情况下石墨烯有较高导热率，控温498 K时导热率比低温情况导热率高；高的环境压力表现出高的导热率，在0.4 MPa时导热率比低的环境压力下导热率高。
4) 石墨烯薄膜可以改变摩擦区水平方向的温度梯度，加速热量传导，使摩擦区导热性更好，高温下这种性能更明显。
[1] 陈杰,张君松,丁行武,等.面接触摩擦副织构化表面润滑特性优化研究[J].兵器装备工程学报,2018,39(05):192-195.
[2] BALANDIN A A,GHOSH S,BAO W,et al.Superior Thermal Conductivity of Single-Layer Letters Graphene[J].Nano,2008,8(3):902-907.
[3] Ota J,Hait S K,Sastry M I S,et al.Graphene Disper-sion in Hydrocarbon Medium and Its Application in Lubricant Technology[J].RSC Advances,2015,5:53326-53332.
[4] 王建,苗瑞珍,刘伟帅，等.石墨烯/聚乙烯界面导热性能的分子动力学模拟[J].分子科学学报，2016，32(02)：156-161.
[5] 惠治鑫,贺鹏飞,戴瑛，等.硅功能化石墨烯热导率的分子动力学模拟[J].物理学报,2014,63(7):184-190.
[6] HU Jiuning,RUAN Xiulin,CHEN,et al.Thermal conductivity and thermal rectification in graphene nanoribbon-s:A molecular dynamics study[J].Nano Letters，2009.
[7] JU S P,CHEN K Y,LIN M C,et al.Investigation of the thermal behaviors of folded graphene by nonequili-brium moleculardynamics simulation[J].Carbon,2014,77(10):36-43.
[8] LIU C,WEI Z,CHEN W,et al.Pressure effects on the thermal resistance of few-layer graphene[J].Physics Letters A,2015,380(1/2):248-254.
[9] 魏志勇,毕可东,陈云飞.石墨烯纳米带热导率的分子动力学模拟[J].东南大学学报(自然科学版),2010,40(02):306-310.
[10] 张程宾,程启坤,陈永平.分形结构纳米复合材料热导率的分子动力学模拟研究[J].物理学报,2014,63(23):255-260.
[11] 郑伯昱,董慧龙,陈非凡.基于量子修正的石墨烯纳米带热导率分子动力学表征方法[J].物理学报,2014,63(07):278-284.
[12] YANG P,TANG Y,YANG H,et al.Thermal management Performance of bent graphene nanoribbons[J].Rsc Advances,2013,3(38):17349-17354.
[13] ZHAN H,ZHANG Y,BELL J M,et al.Suppressed ThermalConductivity of Bilayer Graphene with Vacancy-Initiated Linkages[J].Journal of Physical Chemistry C,2015,119(4):1748-1752.
[14] 李特,雷春丽.考虑气隙变化的高速电主轴热特性仿真[J].浙江大学学报(工学版),2016,50(5):941-948.40-42.
[15] CHOY C L,GREIG D.The low-temperature thermal conductivity of a semi-crystalline polymer,polyethylene terephthalate[J].Journal of Physics C:Solid State Physics,2001,8(19):3121.
[16] YARIFARD M,DAVOODI J,RAFII-TABAR H.Computation of the thermal resistance in graphene sheets with a rectangular hole[J].Computational Materials Science,2017,126:29-34.