考虑弹塑性变形阶段的干气密封接触模型

杨小成,丁雪兴,陈金林

(兰州理工大学 石油化工学院，甘肃 兰州 730050)

干气密封源于上世纪60年代，是一种通过气体动压效应来达到密封效果的非接触密封技术，因其具有较高的密封特性，被广泛地应用在一些精密设备中，如航天发动机等[1-3].干气密封摩擦端面之间的接触通常发生在启停阶段，但在实际情况中，由于人为加工误差、外界环境激励等情况下，会导致动静环局部的振动接触，进而引起整个端面摩擦振动，降低干气密封的密封性及可靠性，甚至失效，造成严重泄露等后果[4-5].因此，开展干气密封摩擦界面的接触特性分析，是提高干气密封性能的关键，对密封系统的评估与优化具有重要的理论指导意义.

研究干气密封摩擦界面的接触特性，首先要对其微观接触模型进行探究.目前，国内外学者对于微观接触模型的研究主要分为三种：统计接触模型、分形接触模型及有限元接触模型[6-8].1966年，Greenwood和Williamson等[9]首次提出将结合面上微凸体高度分布假设为随机分布函数，并提出了经典的GW模型，该模型与前人假设微凸体高度分布遵循等高分布的模型相比更为精确，但该模型只考虑了压入硬度较小的接触面，仅适合弹性接触，具有明显的局限性[10]；后来，Chang等[11]通过对微凸体的力学研究，提出了微凸体还存在着不同于GW模型的非弹性变形阶段，并基于塑性流动理论考虑了非弹性变形阶段，建立了经典的CEB模型，但该模型在塑性临界点处却存在着跳跃点，不满足连续光滑的要求；为解决CEB模型的跳跃点问题，Zhao等[12]对其弹性和塑性变形阶段采用多项式拟合，建立了经典的ZMC模型，该模型虽然实现了接触载荷之间的连续光滑过渡，但是通过该模型建立的接触刚度却是不连续的；随着学者不断探究，其接触模型逐渐被完善，如KE模型[13-14]和LIN模型[15]等.但随着这些模型的不断建立，其统计学模型的弊端也逐渐被体现出来：统计学模型的参数在很大程度上取决于所采用设备的分辨率，不同设备的分辨率所得出的统计参数值是不同的，通过设备得出的统计参数不具备唯一性[16-17].为解决这种弊端，Majumdar等[18]学者于1991年建立了具有分形特性的MB模型，该模型通过分形理论来表征接触表面的轮廓曲线，且与GW模型相比，该模型将接触微凸体顶端曲率看做是接触面积的函数，随着接触面积的变化曲率半径不断发生变化，修正了GW模型中曲率半径为常数时造成的接触面积较小时发生弹性变形，接触面积较大时发生塑性变形的结论[19]；后来，很多学者在MB模型的基础上对弹塑性阶段进行了建模[20-22]，但这些模型的弹塑性阶段采用的都是KE模型，该模型将弹塑性划分为两个变形区间，采用幂函数经验公式对有限元求解的理论值进行了拟合.因此，在临界点处会出现接触载荷不连续的情况，不满足变形连续性的要求[23-24].此外，有些模型在满足接触载荷连续性的同时却满足不了接触刚度的连续性，如ZMC模型等.而且目前，对于干气密封摩擦界面接触模型的力学特性相关的研究鲜有报道，文献[6]和[7]都仅考虑弹性阶段和塑性阶段，未对复杂多变的弹塑性阶段进行探究.

针对上述不足，基于分形接触理论及微观接触力学理论，充分考虑微凸体变形的3个阶段，通过余弦函数构建了干气密封全阶段接触模型，并分别与经典接触模型和试验数据作比较，验证了本文中接触模型的正确性.最后经过理论分析与数值求解得到干气密封摩擦界面的力学特性及其影响因素，并探究了干气密封摩擦接触表面的力学特性与其影响因素之间的变化规律.

1 摩擦界面接触模型的建立

干气密封结构示意图如图1(a)所示，结构图中包含有两对摩擦副，每组摩擦副由动环和静环组成，如图1(b)所示.动环和静环为取得更好的动压效果，常采用硬硬组合或硬软组合.本文中选取常见的硬软组合，即碳化硅和石墨组合[6-7].由于接触变形发生在较软材质的接触面上，因此，本文中只考虑静环上的接触变形，将静环表面等效为柔性粗糙平面，而将硬质动环等效为刚性光滑平面，如图2所示.

Fig.1 Schematic diagram of dry gas seal structure and seal ring图1 干气密封结构及密封环示意图

Fig.2 Schematic diagram of the contact between smooth rigid surface and rough surface图2 光滑刚性表面与粗糙表面接触示意图

1.1 单微凸体接触模型

本文中接触模型的假设如下：假设表面形貌为各向同性；微凸体与平面接触为半球体与平面接触；微凸体之间不发生相互作用力，且微凸体基底位置不发生任何变化；粗糙表面接触过程中不考虑表面热变形等一些极端情况.如图3所示，以微凸体中心为原点，建立x-z坐标轴，图中R为微凸体变形前最大曲率半径，变形前的接触高度为ω，l为微凸体基底尺寸，微凸体受刚性光滑平面挤压作用将会发生屈服变形，变形量为δ，l′为微凸体在受到挤压变形时的最大接触直径.

Fig.3 Schematic diagram of asperity deformation图3 微凸体变形示意图

由W-M函数可知，粗糙表面可以近似等效为由一系列不同尺度的波长相互叠加而成[18]，则粗糙表面轮廓曲线可以表示为[25]

式中，G为特征尺度，D为分形维数.

微凸体最大曲率半径为

微凸体变形前的高度为

微凸体在受到载荷后，会发生接触变形，产生的变形量为

当微凸体受到的载荷小于微凸体的屈服极限时，发生弹性变形，该阶段的力学特性可由Hertz接触力学理论分析求解；当微凸体继续受到载荷作用时，变形量和接触面积也会增加，当接触面积大于弹性临界接触面积时，微凸体进入弹塑性变形阶段，其内部部分区域发生塑性变形，但其周围仍然存在着大量的弹性区域[26].弹性区域随着微凸体所受载荷的不断增大，其弹性比例不断减少，塑性区域不断增加；最后，弹性区域全部转换为塑性区域，微凸体开始进入完全塑性变形阶段.

1.1.1 完全弹性变形阶段

微凸体受刚性光滑平面挤压作用将会发生屈服变形进入弹性变形阶段.根据Hertz弹性接触变形理论可得，该阶段的接触特性(Ae、Fe和Ke)与变形量δ之间的关系如下所示[9]：

式中，E为复合弹性模量，由微凸体两接触表面确定，具体表达式为E=((1-ν12)/E1+(1-ν22)/E2)，其中ν1和ν2为摩擦界面的泊松比，E1和E2为摩擦界面的弹性模量，下标e表示微凸体在弹性变形阶段的接触特性(接触面积、接触载荷及接触刚度).

联立式(5)、式(6)、式(7)及式(2)可得弹性阶段的接触特性(Fe和Ke)与接触面积a之间的关系：

弹性变形的临界值δec表达式为[5]

式中，K为最大接触因子，是摩擦界面泊松比的函数，K=0.545+0.41ν；H为摩擦界面的接触硬度；θ为材料属性，与接触材料自身有关，θ=H/E.

联立式(5)、式(6)及式(8)可得微凸体弹性临界接触面积Aec、临界接触载荷Fec及临界接触刚度Kec：

式中，下标ec表示微凸体在弹性变形阶段的临界接触特性(接触面积、接触载荷及接触刚度).

1.1.2 完全塑性变形阶段

当变形量大于等于弹塑性变形临界阶段时，即δ≥δpc(δpc=110δec)，微凸体变形超过其屈服极限，开始进入塑性变形阶段，该阶段的接触载荷Fp和接触刚度Kp可由Abbott等[27]的塑性接触理论推导得出：

式中，下标p表示微凸体在塑性变形阶段的接触特性(接触面积、接触载荷及接触刚度).

联立式(14)和式(15)可得微凸体弹性临界接触面积Apc、接触载荷Fpc及接触刚度Kpc：

式中，下标pc表示微凸体在塑性变形阶段的临界接触特性(接触面积、接触载荷及接触刚度).

1.1.3 弹塑性变形阶段

通过Kogut等[13]对微凸体变形阶段的研究可知，当变形量处于δec≤δ≤δpc时，微凸体接触面积则为aec≤a≤apc，微凸体进入弹塑性变形阶段.该变形阶段是弹性变形与塑性变形的混合阶段，无法直接精准的进行定量描述.因此，就需要对该阶段进行力学特性分析，通过插值函数来拟合出满足该阶段的力学特性函数.

考虑到弹塑性区域的两个特殊点即弹性变形临界点和塑性变形临界点，因此，接触载荷与真实接触面积a关系式可用余弦函数方程表示为

式中，下标ep表示微凸体在弹塑性变形阶段的接触特性(接触面积、接触载荷及接触刚度).

由式(12)可得接触载荷与接触面积函数表达式在弹性和塑性临界点处的斜率为

为保证弹塑性变形阶段满足光滑连续性，根据光滑连续条件，函数式应满足以下关系式

联立式(19)、(20)和(21)即可解出A、B和C的表达式：

将A、B和C的求解值带入到(19)式中可得：

同理，通过以上步骤，可以求解得到(24)式：

接触刚度Kep可由定义求解：

1.2 粗糙表面接触模型

1.2.1 面积分布密度函数

对于面积分布密度函数，国内外学者将其分为两类，一类是基于统计学原理的高斯分布，另一类是基于分形原理的岛屿面积分布函数.岛屿面积分布函数最早由Mandelbrot在研究海岸线与岛屿面积关系时提出，后经Wang等[28]学者的不断改进，逐渐趋于完善.本文中采用的为改进后的岛屿面积分布函数，如下所示：

式中，a为真实接触面积，al为单微凸体的最大接触面积；φ为区域扩展系数；区域扩展系数φ与分形维数D满足超越方程，其关系如下所示[29]：

1.2.2 粗糙表面整体建模

根据前文所建立的单微凸体接触模型与面积分布函数可得粗糙表面整体接触模型.当微凸体最大接触面积大于塑性接触面积临界值时(即a≥apc时)，粗糙表面整体接触特性如下所示：

式中，下标er、epr和pr分别表示粗糙表面在弹性阶段、弹塑性阶段和塑性阶段的整体接触特性(接触面积、接触载荷及接触刚度).

2 模型验证

本文中先采用理论验证，通过与经典模型对比分析，来证明本文中的模型在弹性、弹塑性及塑性变形全阶段的有效性；最后，再将本文中模型与Jiang等[30]的试验数据及未考虑弹塑性阶段的JZZ模型作对比.其中，表面接触参数选取文献[30]中的实测值进行模型验证，其它材料属性参数选取文献[6]和文献[7]中石墨环的属性值.具体模型参数列于表1中.

表1 接触模型参数Table 1 Contact model parameter

2.1 理论验证

本文中通过与GW模型、KE模型以及ZMC模型三大经典模型作比较，验证本文中所建立模型的正确性.如图4所示，本文中模型与GW模型、KE模型以及ZMC模型三大经典模型变化趋势一样，无量纲接触面积与无量纲接触刚度呈非线性增加.图中GW模型只考虑了弹性变形阶段，未对其他变形阶段进行探究，因此该模型的变化趋势与其他模型有较大的偏差；KE模型虽然在弹性阶段和弹塑性阶段实现了刚度的连续性，但是在塑性临界点处(a/ac=220)却存在着跳跃点，不满足变形全阶段连续性的要求；ZMC模型虽然在KE模型塑性临界点处(a/ac=220)实现了连续性，但是该模型在ZMC塑性临界点处(a/ac=108)却是不连续、不光滑的.而本文中模型对接触刚度的预测曲线无论在KE塑性临界点处(a/ac=220)还是ZMC塑性临界点(a/ac=108)处都实现了连续性和光滑性，符合接触变形的要求.

Fig.4 Relationship between dimensionless contact stiffness and contact area图4 无量纲接触刚度与接触面积关系图

2.2 试验验证

Jiang等[30]为探究接触刚度与接触压力之间的变化关系，建立了JZZ经典模型，并对其模型的预测结果与试验数据作了比较，发现JZZ模型与试验数据虽然存在一定的误差，但是趋势相同，证明了所建立JZZ模型的正确性.本文中引用Jiang等[30]的试验数据及所建立的JZZ模型与本文中模型作对比.如图5所示，本文中模型与Jiang等[30]所建立的JZZ模型的变化趋势相同，更接近试验数据的趋势.本文中模型在法向接触压力小于1.8 MPa时，预测值大于试验值，当法向接触压力大于1.8 MPa时，预测值小于试验值.而JZZ模型预测值均小于试验值.JZZ模型与本文中模型相比，接触刚度的趋势与试验值偏离较大，产生这种较大偏差的原因是JZZ模型未充分考虑微凸体弹塑性变形阶段，忽略了其弹塑性过渡阶段对微凸体的力学作用.因此，通过Jiang等[30]的试验数据及JZZ模型可知，本文中模型具有较好的有效性及预测性.

Fig.5 Comparison between experimental data and this model图5 试验数据与本文模型对比图

3 理论结果与分析

基于前文推导及验证的接触模型，将实际工况下的接触参数代入本文中模型中，通过计算结果可以分析各参数对摩擦界面的力学特征，并探究接触参数与力学特征之间的影响规律.本文中摩擦界面采用石墨和碳化硅材料，具体参数列于表2中[6-7,30].

表2 接触模型参数Table 2 Contact surface structure parameters and material parameters

3.1 分形维数对接触载荷与刚度的影响

为探究干气密封的真实接触面积Ar和分形维数D对接触载荷Fr与接触刚度Kr的影响，采用表2中的接触参数进行计算，其中，特征尺度G取1×10-12m，分形维数D从1.2起每隔0.15取值，分别为1.2、1.35、1.5、1.65和1.8五种情况，结果如图6和图7所示.从图6可以看出无量纲接触载荷Fr*与分形维数D、无量纲真实接触面积Ar*呈正相关.分形维数D较小时，无量纲接触载荷Fr*较为相近，几乎呈1条水平直线，而分形维数D较大时，无量纲接触载荷Fr*呈数量级的增大，之间的差距也逐渐增大，为突出这种数量级的变化趋势，取无量纲接触载荷Fr*的对数lg(Fr*)时，如图6左上角小图所示，取分形维数D=1.2时，无量纲接触载荷Fr*的最大数量级为1个数量级，而取分形维数D=1.8，无量纲接触载荷Fr*的最大数量级达到3个数量级，即无量纲接触载荷Fr*的数值变化范围在0到103之间.

Fig.6 Relationship between Fr* and Ar*(G =10-12)图6 Fr*与Ar*关系图(G =10-12)

Fig.7 Relationship between Kr* and Ar*(G =10-12)图7 Kr*与Ar*关系图(G =10-12)

无量纲接触刚度Kr*和无量纲接触载荷Fr*相同，都与无量纲真实接触面积Ar*和分形维数D呈正相关，如图7所示，无量纲接触刚度Kr*随无量纲接触面积Ar*的增大开始急剧增加，最后趋于平缓.无量纲接触刚度Kr*在分形维数D较小时，呈线性增长，当分形维数D较大时，无量纲接触载荷Kr*呈非线性增长，且呈数量级的增大.为突出这种数量级的变化趋势，取无量纲接触刚度Kr*的对数lg(Kr*)，如图7中的小图lg(Kr*)-Ar*所示，取分形维数D=1.2时，无量纲接触刚度Kr*很小，当分形维数D=1.2时，无量纲接触刚度Kr*最大数量级达到了3个数量级，且与无量纲接触载荷Fr*相比，无量纲接触刚度Kr*的数量级变化范围较小，均在1个数量级内.

因此，可通过增加干气密封摩擦接触面的真实接触面积Ar和分形维数D来提高密封性能，且分形维数D的改变对干气密封性能的提升要比表面真实接触面积Ar的影响更大，尤其在分形维数D较大时.因此，可保障在一定的真实接触面积的情况下，尽量提高接触面的分形维数，这样能大幅度提升干气密封的密封性能.

3.2 特征尺度对接触载荷与刚度的影响

为探究干气密封的真实接触面积Ar和特征尺度G对接触载荷Fr与接触刚度Kr之间的影响规律，取分形维数D为1.5时，特征尺度G每隔1个数量级取值，分别取1×10-14、1×10-13、1×10-11和1×10-10m进行计算分析，如图8和图9所示.由图8可知，无量纲接触载荷Fr*随着无量纲真实接触面积Ar*的增大而逐渐增大，且增加幅度依次增大.特征尺度G越大，无量纲接触载荷Fr*的值越小，变化范围也越小；当特征尺度G取1×10-10m时，无量纲接触载荷Fr*整体趋势几乎是1条水平的直线，而当特征尺度G取1×10-14m时，无量纲接触载荷Fr*变化幅度较大，且增长趋势与指数函数较为相似.这是因为特征尺度D可以表征接触面的粗糙程度，当特征尺度G越大时，表面粗糙度就越大，表面微凸体高度分布差异也就越大，单位面积发生塑性接触的微凸体就越多，接触面承载能力就越差，因而无量纲接触载荷Fr*随特征尺度G的增大而减小.

Fig.8 Relationship between Fr* and Ar* (D =1.5)图8 Fr*与Ar*关系图(D =1.5)

Fig.9 The relationship between Kr* and Ar* (D=1.5)图9 Kr*与Ar*关系图(D =1.5)

特征尺度G对无量纲接触刚度Kr*的变化规律与特征尺度G对无量纲接触载荷Fr*的变化规律相似，都与特征尺度G呈负相关，与无量纲真实接触面积Ar*呈正相关.但两者变化趋势有所不同，无量纲接触刚度Kr*随无量纲接触面积Ar*的增大呈先急剧增加，后趋于平缓.随特征尺度G的增大，无量纲接触刚度Kr*等数量级减少，且无量纲接触刚度Kr*的变化范围都保持在1个数量级之内，如图9所示.

因此，可通过增加干气密封摩擦接触面的真实接触面积Ar和减小特征尺度G来提高密封性能，且摩擦接触面的特征尺度G在真实接触面积Ar较大时相比Ar较小时影响更大.因此，可保障在一定的真实接触面积的情况下，尽量减少接触面的特征尺度，这样能大幅度提升干气密封的密封性能.

4 结论

a.本文中基于分形理论和微观接触力学理论，考虑传统的弹性变形和塑性变形阶段外，基于余弦函数考虑了弹塑性变形阶段，建立了具有光滑连续特性的干气密封摩擦界面全阶段接触模型，修正了传统接触模型在弹塑性变形阶段的不连续性与不光滑等缺陷，并通过与经典模型和试验数据的双重验证，证明了本文中模型的正确性.

b.无量纲接触载荷Fr*和无量纲接触刚度Kr*都与分形维数D、无量纲真实接触面积Ar*呈正相关，且他们的数量级与分形维数D也呈正相关，分形维数D较小时，无量纲接接触刚度Kr*增长较缓慢，数量级较为接近.但当分形维数D较大时，无量纲接接触刚度Kr*呈非线性增长，数量级之间的差距也逐渐增大，且无量纲接接触载荷Fr*变化范围相对较大.

c.无量纲接触载荷Fr*与无量纲接触刚度Kr*都与特征尺度G呈负相关.当特征尺度G以1个数量级递增时，无量纲接触载荷Fr*与无量纲接触刚度Kr*的变化范围都在1个数量级之内.

d.可通过增加干气密封摩擦接触面的分形维数D和真实接触面积Ar及减小特征尺度G来提高密封性能.一般而言，干气密封摩擦界面的分形维数D应大于1.5，特征尺度G应尽量小于1×10-12m.