机械密封为端面密封，其依赖两个固体表面间的接触和相对运转，以防止流体泄漏。常规设计中通常按照接触端面的名义面积来计算其端面比压。然而，机械加工表面在本质上都是粗糙的，动、静环端面的实际接触仅发生在部分微凸体上，导致实际接触面积仅占名义接触面积很小的一部分^[1]。干摩擦是机械密封最恶劣的润滑状态，其接触状态直接决定了密封的性能与寿命^[2]。张永振^[3]指出在干摩擦状态下，摩擦热量与速度、接触压力存在幂函数的关系，外载荷和速度对端面的接触、摩擦磨损行为的影响尤为突出。所以，研究外载荷和速度对机械密封干摩擦时端面接触特性的影响对了解其密封机理和端面磨损具有重要意义。

目前，对相对运动摩擦端面的动态过程的研究仍然缺乏直接观察和检测评估的手段，因此，为了准确描述摩擦端面的接触特性，需要建立合理的接触模型，前人对此进行了大量研究。Majumdar等^[4-5]首先基于2维W–M分形函数建立了弹塑性接触分形模型，探究了总载荷与真实接触面积的关系式；Wang等^[6-7]在M–B模型的基础上进行改进，得到了M–B分形接触修正模型；孙见君^[8]、魏龙^[9]、丁雪兴^[10]等参照M–B模型建立了机械密封摩擦端面的分形模型，通过2维理论计算对密封端面泄漏、真实接触面积、温度分布进行了研究。基于分形理论建立的接触模型在各个领域都得到了广泛运用^[11-16]。可见，利用分形理论可以实现对粗糙面的唯一表征，不受尺度的影响。

干摩擦的宏观表现行为是众多接触摩擦学行为耦合作用的结果^[3]。粗糙表面上的接触压力、热应力、温度场是相互影响的，因此，考虑微凸体间的相互作用及热–力耦合对分析摩擦端面的接触特性尤为重要。孙星星等^[17]考虑密封动静环的热力变形，数值求解了密封环的温度分布，但忽略了粗糙面微凸体的影响。房桂芳等^[18]依据分形理论提出了机械密封端面摩擦热的耦合计算方法。Li等^[19]提出了考虑磨损和热变形的接触机械密封分形泄漏模型。黄健萌^[20-21]、丁雪兴^[22]等建立了粗糙面/平面的1维滑移接触模型，通过有限元模拟求解了滑动摩擦过程中微凸体变形及热–力耦合对粗糙面应力、温度变化的影响。总的来说，对于机械密封摩擦端面接触的研究主要在理论计算方面，模拟计算也仅采用了直线滑动摩擦的方式；对于相互接触的表面而言，直线运动与回转运动的摩擦学特性是不同的^[3]，因此采用3维回转接触模型进行仿真将更符合实际。

针对碳石墨与碳化钨配对的机械密封，首先，作者基于W–M分形函数建立了模拟其摩擦副端面的3维粗糙实体与理想光滑刚体的转动摩擦接触模型；然后，在考虑热–力耦合及微凸体间相互作用的基础上，通过ABAQUS有限元分析软件求解了其瞬时干摩擦时密封端面间的接触特性，并讨论不同外载荷、滑动速度对其特性的影响，并针对API682中串联式机械密封的第2级接触式干摩擦机械密封（抑制机械密封）进行了分析，以期为探索机械密封在干摩擦时的发热、摩擦行为和密封设计提供借鉴。

1 计算模型的建立 1.1 粗糙表面及转动接触模型的建立

表征粗糙接触界面的形貌是摩擦学的热点，Yan等^[23]的研究给出了能够描述粗糙表面微观特征的3维W–M函数：

$\begin{aligned}[b] &{\textit{Z}}\left( {x,y} \right) = L{\left( {\frac{G}{L}} \right)^{\left( {D - 2} \right)}}{\left( {\frac{{\ln\; \gamma }}{M}} \right)^{1/2}} \times \sum\limits_{m = 1}^M {\sum\limits_{n = 0}^{{n_{\max }}} {{\gamma ^{\left( {D - 3} \right)n}}} } \times \bigg\{ \cos\;{\varphi _{m,n}} -\bigg. \\ &\;\;\;\;\cos\bigg[\frac{{2\text{π} {\gamma ^n}{{\left( {{x^2} + {y^2}} \right)}^{1/2}}}}{L} \times\bigg.\cos \left( {{{\tan }^{ - 1}}\left( {\frac{y}{x}} \right) - \frac{{\text{π} m}}{M}} \right) + {\varphi _{m,n}}\bigg]\bigg\} \end{aligned}$

(1)

式中：Z(x，y)为3维粗糙表面形貌的高度；D和G分别为分形维数和分形粗糙度；L为取样长度；γ为表面轮廓频率密度的特征参数，通常取1.5；M为曲面褶皱的重叠数，值为10； ${\varphi _{m,n}} $ 为随机相位；n_max为最大频率，其值为：

${n_{\max }} = {{\rm{int}}} \left[ {{{{{\lg}}\left( {\frac{L}{{{L_{\rm{s}}}}}} \right)} / {\lg\; \gamma }}} \right]$

(2)

式中，L_s为截止长度（近似取材料的原子间距）。

对于同一粗糙表面，Z(x，y)可由与尺度无关的分形参数D与G唯一确定，因此式（1）具有在不同长度尺度上描述粗糙表面的能力^[23]，并且随着D增大，G减小，表面粗糙度降低，趋于光滑^[24-25]。本文中，D取2.35，G取2.45×10^–9 m，此时表面粗糙度Ra的大小为0.8 μm。

由于接触问题属于典型的边界非线性问题，且粗糙表面的特征尺度在微米级，考虑到计算成本，研究中将微元体的模型尺寸也控制在微米级，其中，静环模型：外径D_Ai=260 μm，内径D_Ao=380 μm，厚度h_A=60 μm；动环模型：外径D_Bi=250 μm，内径D_Bo=390 μm，厚度h_B=60 μm。通过增大转速的数量级可以实现模型微元体的真实线速度V，根据文献[26]粗糙表面轮廓曲线的局部与整体的自相似性可知，微元体的接触行为可以反映实际端面之间的接触行为。

机械密封中，动、静环材料常采用硬软配对，摩擦副为硬质合金（WC）对碳石墨（M106K），两者的弹性模量相差1～2个数量级^[22]，见表1。因此，研究中将静环（实体A）与动环（实体B）的接触简化为分形粗糙表面与刚性理想光滑平面的接触^[9]；同时，由于密封环模型轴对称，为了减少计算时间，静环采用1/8周期模型，如图1所示。其中：实体A中，A₂为顶面，A₁为底面（分形粗糙面），A₃、A₄为侧面，A₅为内圆柱面，A₆为外圆柱面；实体B中，B₁为顶面（接触面），B₂为底面，B₃为内圆柱面，B₄为外圆柱面；RP为参考点；径向X、周向Y、轴向Z（实体A、B厚度方向）为坐标系。

1.2 位移边界条件

计算时设置两个时间段：在第1个时间段内，将均布载荷由0线性增加到P，加载到实体A的A₂面上，同时对其4个侧面施加沿X、Y方向的位移约束。对于实体B，约束其Z方向的位移，同时将实体B用参考点RP（0，0，0）进行刚性约束（设置实体B为刚性体，同时可以通过RP的运动控制实体B的运动），通过约束RP的6个自由度使B保持静止，让两实体建立稳定的接触关系。在第2个时间段内，保持A₂面上载荷P不变，同时给定实体B的参考点RP绕Z轴的转动速度ω，使两实体间发生转动摩擦。

1.3 模型假设

由于计算条件的限制，本文在计算中做出如下假设：

1）滑动摩擦过程遵循库仑定律，在计算时间范围内认为摩擦系数保持不变，不受外载荷、滑动速度变化的影响。

2）构成实体A、B的材料各向同性，且在瞬态计算中，由于温升不大，假设材料的热物性参数不会改变。

3）忽略材料磨损的影响，认为摩擦消耗的功全部转化为热能，并被摩擦副吸收，则局部热流密度可表示为：

$q\left( {r,\varphi ,t} \right) = f{P_{\rm{b}}}\left( {r,\varphi ,t} \right)V\left( {r,\varphi ,t} \right)$

(3)

式中，q为热流密度，f为摩擦系数，P_b为粗糙表面微元体的接触压力，V为相对滑动速度，t为时间，r为半径。

4）考虑摩擦热在两实体间的分配，认为接触区域为理想热传导，即接触点处两实体的瞬时温度相等，这样摩擦热流密度就能根据材料的热物性参数及散热环境在两实体间自由分配^[20]，同时忽略热辐射的影响。

1.4 热边界条件

实体A、B的热传导方程可表示为：

$\begin{gathered} \rho c\frac{{\partial T}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {kr\frac{{\partial T}}{{\partial r}}} \right) + \frac{1}{{{r^2}}}\frac{\partial }{{\partial \varphi }}\left( {k\frac{{\partial T}}{{\partial \varphi }}} \right) + \frac{\partial }{{\partial {\textit{z}}}}\left( {k\frac{{\partial T}}{{\partial {\textit{z}}}}} \right) \end{gathered} $

(4)

式中，ρ为实体A、B材料的密度，c为比热容，k为热导率，T为温度，t为时间，r为半径。

初始条件：T_A=T_B=298 K，t=0。

耦合条件（假设4））：T_A=T_B，q_A + q_B = q（实体A、B接触点）。

接触面（A₁、B₁）的热边界条件为：忽略热辐射，非接触区域的传热为对流换热；接触区域为摩擦热流的输入^[21]：

$\left\{\!\!\!\! \begin{array}{l} {k_{{{\rm{A}}}}}\dfrac{{\partial {T_{{{\rm{A}}}}}}}{{\partial {\textit{z}}}} = - (1 - g(m)){h_{{{{\rm{A}}1}}}}\left( {{T_{{{\rm{A}}}}} - {T_{{m}}}} \right) + g(m){q_{{{\rm{A}}}}},\\ {k_{{{\rm{B}}}}}\dfrac{{\partial {T_{{{\rm{B}}}}}}}{{\partial {\textit{z}}}} = - (1 - g(m)){h_{{{{\rm{B}}1}}}}\left( {{T_{{{\rm{B}}}}} - {T_{{m}}}} \right) + g(m){q_{{{\rm{B}}}}} \end{array} \right.$

(5)

式中，接触点处g(m)=1，非接触点处g(m)=0。

由于本文研究的是瞬态问题，在很短的时间内，热量只在两实体接触面的表层及亚表层传递，因此可以认为密封环非接触面（A₂、B₂）是绝热的；同样地，干摩擦条件下可以认为实体A及实体B的侧面A₃～A₆、B₃～B₄是绝热的。

1.5 求解方法

采用ABAQUS/Explicit求解计算模型，对于复杂的接触问题，显示方法求解效率较高且结果可靠。其采用中心差分的方法求解运动方程，用一个时间段的运动学条件计算下一个时间段的运动学条件，通过自动划分足够小的时间增量段来确保模型状态稳定，从而获得高精度的解^[27]。求解前需要设置接触对，将刚性光滑面B₁设置为主面，分形粗糙面A₁设为从面，采用罚函数法作为求解接触问题的算法，两表面间的相对滑动通过有限滑移公式描述。

考虑热–力耦合计算，使用C3D8T六面体热–力耦合单元对整个计算模型进行网格划分，并对实体A的粗糙面（A₁）的表层及亚表层进行网格细化，对光滑实体B使用较粗的网格划分方式。为了证明求解结果与网格大小无关，对模型进行网格无关性验证分析，由于主要研究对象是粗糙表面，所以采用调整粗糙表面网格大小的方式验证接触面积的变化。定义无量纲接触面积为A_r/A_a，其中，A_r为真实接触面积，A_a为粗糙面A₁的名义接触面积。ABAQUS软件中可设置接触面积的监控，计算完成后输出A_r随时间的变化，则无量纲接触面积随网格尺寸的变化如图2（a）所示。可以发现，随着网格尺寸减小，无量纲接触面积趋于稳定，实际粗糙面网格尺寸为1.736 μm，可以较好地保证计算结果的准确性。网格划分完成后，光滑实体单元数目为10400，粗糙实体单元数为20475，如图2（b）、（c）所示。

2 算例分析 2.1 计算参数选择

动环、静环的材料如表1所示，两实体的干摩擦系数选择为f=0.25^[28]，静环材料的屈服强度为σ_y=200 MPa^[9]。在初期计算中，外载荷P过小时，数据变化不明显，故P选取变化范围为0.2～40.0 MPa；滑动速度V为光滑实体B外径处的线速度，变化范围为20～50 m/s。

在模拟计算过程中不考虑材料非线性问题，只考虑几何非线性与接触边界非线性，由于计算过程难以收敛，同时粗糙表面实际处于不断磨损的状态中，表面形貌也不断更新^[24]，因此，只求解转动发生后60 μs以内的瞬态过程。为了使两实体发生平稳接触，将转动前的计算时间设置为40 μs，外载荷以线性增加的方式施加在A₂面上直到载荷等于P为止，之后保持载荷不变使两实体发生相对转动，速度为V；密封环初始温度与环境温度设置为298 K。每次计算改变外载荷P、滑动速度V两个参数中的一个值，计算得到密封环不同工况下的7组摩擦数据。

2.2 结果分析与讨论 2.2.1 真实接触面积

定义无量纲载荷为P/E´，其中：

${E^{'}} = {\left( {\frac{{1 - \nu _1^2}}{{{E_1}}} + \frac{{1 - \nu _2^2}}{{{E_2}}}} \right)^{ - 1}}$

(6)

得到本文分形模型与魏龙分形模型^[9]、M–B分形模型^[5]理论值的对比如图3所示。

从图3中可以发现，随着外载荷增大，真实接触面积近似呈线性增大，外载荷增加3倍时，接触面积增大了2.33倍，这是因为外载荷的增加导致粗糙表面受力增加，接触面主峰微凸体发生弹塑性变形，次峰微凸体逐渐参与接触，增加了接触面积。同时，可以看出：随着载荷的增大，模拟分析的结果与理论值之间具有相同的变化趋势，在载荷较低时，它们间差异较小（串联式机械密封第二级接触式干运转机械密封恰恰是在低载荷下工作）；随着载荷的进一步增大，它们间的误差也逐渐增大（无量纲载荷为0.0015时，M–B模型计算得到的无量纲接触面积是本文模拟计算的1.364倍，是魏龙模型的1.186倍）。产生这一差异的原因在于：魏龙分形模型、M–B分形模型的建立是基于2维粗糙轮廓的，而不是真实的3维粗糙表面，它们未考虑粗糙表面相邻微凸体之间的相互作用。文献[29]的研究表明，魏龙分形模型、M–B分形模型计算得到的真实接触面积与实验值相比较大。本文所采用的3维有限元模拟考虑了粗糙表面相邻微凸体之间的相互作用，所以其结果更接近实际，从而验证了模拟分析的正确性。

图4为P=30 MPa时，不同滑动速度下，滑动0.3 μs及60.0 μs时粗糙面的无量纲接触面积。可以发现，不同滑动时间下，接触面积随速度增加的斜率相同，速度越大，接触面积越大。但在滑动速度增加1.5倍的情况下，接触面积仅增大为原来的1.002倍，因此，粗糙面的真实接触面积的大小主要取决于外载荷，而滑动速度影响较小。

2.2.2 接触压力

摩擦磨损是表面行为，因此，研究粗糙表面的真实接触面积、接触压力、温度的分布极其重要。图5（a）～（d）所示为V=30 m/s时，P=10、20、30、40 MPa下，滑动60 μs时刻分形粗糙面的接触压力分布云图。由图5可以看出：随着外载荷增大，真实接触面积A_r上较高的接触压力节点在增加，粗糙表面上最大接触压力值位于节点13538；当外载荷为10、20、30、40 MPa时，节点13538的接触压力分别为415.88、416.64、423.25和431.59 MPa。对比不同载荷下各个接触区域上接触压力的变化，可以发现外载荷增大3倍后，节点13538上接触压力值增加的幅度并不是很大，这说明相对运动的摩擦副的接触压力峰值达到一定值后变化缓慢，该值对外载荷的变化不敏感，应与摩擦副的材料力学特性有关。

粗糙表面节点13538的接触压力随载荷增加变化缓慢，其原因是接触压力在达到一定值后，材料发生弹塑性变形，次峰微凸体进入接触并承担了额外的载荷，使得节点13538所在微凸体受力得到缓解。表现为粗糙表面的接触压力峰值维持在一个定值范围内，而其他接触节点的接触压力继续增加。换句话说，在本文的研究参数范围内，接触压力峰值是“自限性”的，这种特征保证端面在较大的外载荷下的接触是由点接触到局部面接触、一个节点接触到众多节点接触。更多的节点参与接触，端面受力条件得到改善，这种“自限性”对于密封环端面微凸体受力状况的改善是有利的。

图6（a）～（d）所示为P=30 MPa时，V=20、30、40、50 m/s下，滑动60 μs时刻分形粗糙面的接触压力分布云图。

由图6可以发现，在外载荷相同时，接触压力的分布基本相同，最大接触压力分别为430.34、423.25、428.86、420.26 MPa，均位于粗糙表面节点13538，并随着速度增加具有微小的下降趋势，这是因为滑动期间速度增大时引起的微凸体局部摩擦热流密度变大，微凸体热变形更大使接触节点局部接触面积略微增大，导致节点最大接触压力略有减小。

2.2.3 温度分布

图7（a）～（d）所示为V=30 m/s，P=10、20、30、40 MPa下，滑动60 μs时刻分形粗糙面的温度分布云图。由于粗糙表面微凸体的形状大小、高低不一，滑动接触正是发生在这些不连续的接触微凸体上。摩擦所引起的热量产生、释放也发生这些区域，使得接触节点处温度较周围较高，形成“热点”。

由图7可以看出，P=10 MPa时高温区主要分布在粗糙表面微凸体接触点位置，且端面上存在明显的温度梯度分布，会有热应力产生。随着载荷增大，粗糙表面接触微凸体及接触点越多，在滑动瞬间形成的“热点”也越多，接触中心区域的温度也越高，粗糙表面整体温度也增高。另外，可以发现，不同载荷下最高温度值所在的节点不同，不是单纯地位于接触压力最大的位置。这表明粗糙表面的最高温度还受到局部热流密度的影响；即最高温度受到局部接触压力和局部滑动速度和局部热传递的共同影响。

随着外载荷的增加，虽然粗糙表面最大接触压力改变较小，但载荷变化使得粗糙表面接触点的数量明显增多，滑动过程中粗糙表面上的摩擦热源增多，导致了端面整体温度升高，故材料的耐摩擦磨损性能也受到影响。

图8所示为外载荷为10 MPa、滑动速度为30 m/s，摩擦滑动300 μs条件下的粗糙面各节点的温度变化曲线（节点分布如图7（a）所示）。由图8可以看出：0～40 μs为外载荷加载阶段，粗糙面各节点温度保持不变；当两表面开始滑动，40～42 μs内由于速度突增及热传导的滞后性出现了温度急升；随后42～340 μs，在自身热传导的作用下，各节点温度开始缓慢上升，温度上升斜率也是一致的，且滑动60 μs的温升斜率与滑动300 μs的温升斜率也相同。这也验证了本文采用摩擦滑动60 μs的计算时间进行研究是可行的。

API682中规定串联式机械密封第2级接触式干运转机械密封的介质压力应控制在0.07 MPa以内。为了阐释该密封在有气体冷却的配合下实现25000 h运转的可能性，本文模拟计算了外载荷0.2 MPa、滑动速度30 m/s条件下滑动60 μs时粗糙面温度分布云图，如图9所示。由图9可以看到，摩擦滑动60 μs后，粗糙面的最高温度仅从环境温度25 ℃升高到25.42 ℃，温升仅有0.42 ℃，粗糙面温度变化不明显。图10为滑动速度30 m/s时，不同外载荷下粗糙面最高温度的变化，随着外载荷的增加，最高温度是线性增大的，与式（3）描述的规律相符。可见降低机械密封的外载荷可以有效防止其干摩擦时石墨环温升过高，API682标准中建议介质压力≤0.07 MPa是有必要的。

图11（a）～（d）所示为P=30 MPa，V=20、30、40、50 m/s下，滑动60 μs时刻分形粗糙面的温度分布云图。由图11可以发现，在相同的外载荷下，对应的粗糙表面最高温度值分别为44.84、49.12、52.99和57.52 ℃，随着滑动速度增加温度呈上升趋势，粗糙表面的温度不均匀趋势明显加剧，各区域温度值均有所提高。根据式（3），速度增加，摩擦发热增加，粗糙表面的温度梯度也随之变大，热应力也会加剧。

3 结　论

在P=0.2～40 MPa、V=20～50 m/s，f=0.25和粗糙表面Ra=0.8 μm的工况下，针对碳石墨–硬质合金配对摩擦副在干摩擦状态下的接触特性进行了分析，研究发现：

1）粗糙端面真实接触面积的大小主要取决于外载荷。外载荷增加真实接触面积近似线性增大，但滑动速度影响较小。

2）随着外载荷与滑动速度的增加，石墨粗糙表面微凸体接触压力峰值在415～432 MPa范围内，表现出一种“自限性”。

3）粗糙端面温度分布不均匀，局部温度较高；端面最高温度并非位于接触压力最大的位置，而是受到局部接触压力和局部滑动速度及局部热传递的共同影响。

4）外载荷、速度增大均会使端面整体温度增加。

5）API682标准中要求抑制机械密封的介质压力≤0.07 MPa是有必要的。降低介质压力可以降低外载荷，从而可以有效防止干摩擦时密封环温升过高。