盾构机唇形密封流固热耦合仿真研究*

纪佳馨 苏世东 项 冲 郭 飞

（1. 中国石油大学（华东） 机电工程学院 山东青岛 266580；2. 清华大学高端装备界面科学与技术全国重点实验室 北京 100084）

随着国家城市化进程的发展， 盾构掘进技术在城 市地铁、 铁路交通等领域得到广泛应用。 鉴于盾构机工作过程中特殊的隧道环境， 主驱动唇形密封系统[1-3］一旦失效， 会导致外界环境中的土砂倒灌进而破坏主轴承， 严重影响工程进度。 唇形密封件的性能直接决定了整台盾构机的施工效率， 且在施工过程中应避免唇形密封件的更换， 为免影响施工效率， 其故在盾构机主驱动密封设计阶段， 预测唇形密封性能是否满足要求尤为重要。

国产盾构机主驱动密封件多依赖国外进口， 且价格高昂， 为提高国产大型盾构机的国际竞争力， 亟待开展盾构机主驱动密封研究。 谭锋等人[4］利用正交试验法针对盾构机主驱动唇形密封开展结构优化研究，根据密封工作要求确定最优密封结构， 优化后密封件接触宽度降低23%。 张中华等[5］在主驱动密封失效分析的基础上， 总结提出密封安装结构和工艺优化方案， 改善盾构机主驱动密封性能。

目前， 针对盾构机主驱动密封的研究集中于结构和安装优化， 尚无预测密封性能的研究手段。 盾构机主驱动唇封密封介质为润滑脂， 工作时唇口温度可达50～60 ℃， 为更好地预测唇封的密封性能， 本文作者考虑润滑脂流变特性、 唇口温度对流场分析、 密封材料的影响， 建立盾构机主驱动唇形密封流固热耦合仿真模型； 然后利用流速分离法推导润滑脂二维雷诺方程， 采用赫兹接触模型计算粗糙峰接触压力， 结合有限元软件开展热力耦合分析， 实现唇封温度场及摩擦力矩、 泄漏率等关键性能参数的定量预测。 该模型对盾构机主驱动密封设计研究具有重要意义。

1 数值仿真模型

唇封系统主要由密封唇、 旋转轴和密封介质三部分组成， 其相互作用共同决定唇封系统的密封性能。如图1 所示为唇封与轴接触区的微观示意图， 可以看出建立盾构机主驱动唇形密封数值仿真模型， 需要综合分析密封过程涉及的物理要素， 包括密封介质的流体力学分析， 微观粗糙峰与轴之间的接触力学分析和密封唇的固体力学分析。 图中所示x方向为周向，y方向为轴向。

图1 接触区微观示意Fig.1 Microscopic schematic of the contact zone

为方便模型求解， 忽略某些对计算结果影响不大的因素， 提出了以下合理的简化和假设：

（1） 唇封宏观性质不受表面微观变形影响；

（2） 唇封与轴接触特征沿圆周方向周期性分布；

（3） 唇封两侧充满密封介质， 以方便计算泄漏量；

（4） 转速恒定， 唇封与轴长时间磨合， 密封状态稳定；

（5） 膜厚远小于轴径， 故忽略油膜曲率的影响。

1.1 流体力学分析

雷诺方程是解决流体黏性流动的基础理论方程，对于润滑脂这种典型的非牛顿流体[6-7］， 开展流体压场分布的研究仍依靠雷诺方程的建立与求解。 区别于润滑油等牛顿流体， 润滑脂不满足牛顿内摩擦定律，若沿用传统推导方式， 其复杂的本构方程会导致雷诺方程在二维层面难以计算。

文中利用流速分离法[8］推导得到的普适流体润滑方程可以有效解决一般非牛顿流体润滑问题， 普适流体润滑方程如下：

式中：h为流体膜厚；ρ为流体密度；t为时间；f-1为非牛顿流体本构方程的反函数；pf为流体压力；U0、Uh为周向上两边界速度；V0、Vh为轴向上两边界速度；W0、Wh为径向上两边界速度。

目前广泛应用于润滑脂研究的本构模型为Herschel-Bulkley 模型[9］， 然而密封过程中膜厚一般是微米级别的， 给定转速下对应的剪应变率远大于屈服应力， 为计算方便， 可将Herschel-Bulkley 本构模型简化成Ostwald 模型[10-11］， 如下所示：

式中：τf为流体剪应力；η为稠度系数；u为流体周向速度；n为流变指数。

润滑脂的平衡方程为

结合流速分离法可得Ostwald 模型对应的反函数为

上述反函数代入到普适流体润滑方程， 即得到面接触条件下简化的润滑脂二维雷诺方程：

预设量纲一化变量， 防止pf、h等参数量级上的差异导致方程求解误差过大进而不收敛：

式中：pr为参考流体压力；Pf为量纲一流体压力；hr为参考流体膜厚；H为量纲一流体膜厚；lx、ly为唇封和轴接触宽度。

对润滑脂二维雷诺方程进行量纲一化处理：

上述方程中A～D为系数项，Pf为未知项，H为常数项。 若计算网格节点为m×m， 为提高数值求解速度， 可利用矩阵法将上式转换成系数项矩阵（m2×m2） ×Pf列向量（m2×1） ＝常数项列向量（m2×1）的形式， 通过矩阵求逆运算准确快速地得到Pf分布。

在已知pf和h后， 即可通过润滑脂在轴向上的质量流量公式计算泄漏率：

式中：Qy为轴向质量流量。

1.2 接触力学分析

由于物体表面无法绝对光滑， 唇封与轴实际配合时， 会存在微观层面的粗糙峰接触挤压行为， 其产生的接触压力和流体压力共同起到支撑负载的作用。

文中选择经典的赫兹接触模型[12-14］计算粗糙峰接触压力， 该模型为两个弹性接触的表面提供封闭形式的表达式。 其中接触压力计算公式如下：

式中：pc为粗糙峰接触压力；E′ ＝E／（1- ν2） 为等效弹性模量， 其中E为弹性模量，ν为泊松比；d*为粗糙峰的几何重叠深度；R′为等效曲率半径；r*为到粗糙峰接触圆中心点的距离；a*为接触区域半径。

粗糙峰接触产生的摩擦剪应力计算公式如下：

式中：τc为粗糙峰剪应力；f为摩擦因数。

1.3 变形力学分析

实际装配条件下， 唇封还会受预紧力作用发生宏观变形， 产生静态接触压力psc。 对于唇封来说， 正是流体压力、 粗糙峰接触压力和静态接触压力三者的共同影响使得唇口发生变形。

考虑到唇口变形量非常小， 符合小变形理论， 可以认为唇封的变形与施加的载荷呈线性关系， 因此可以采用影响系数法[15］描述变形与载荷的对应关系。其中接触表面任一计算节点i法向变形的计算公式为

式中：δn为法向变形量；In为法向影响系数矩阵；psc为静态接触压力。

可以看出节点i的法向变形量由该点所在轴向列m个节点的法向受力状况共同决定。

同理， 接触表面任一计算节点i切向变形的计算公式为

式中：Is为切向影响系数矩阵。

通过变形力学分析， 计算不平衡应力值下唇封的变形量， 用以修正膜厚， 反过来膜厚又会改变流体压力和粗糙峰接触压力， 引起唇封变形， 故需经过多次迭代， 最终使得应力平衡、 膜厚稳定。

1.4 固体力学分析

文中利用有限元软件ABAQUS 对唇封开展固体力学分析， 建立实际装配受压条件下唇封系统几何模型， 结合材料单轴实验数据， 获取静态接触压力和影响系数矩阵。

仿真具体流程： 首先参考真实尺寸构建唇封二维轴对称模型， 模型包括唇封、 旋转轴、 定位夹具等；其次根据丁腈橡胶单轴试验测试数据定义材料属性，材料评估后选定的应变能函数为Ogden，N＝3， 其中模型系数Mu1 ＝-36.194 852 8， Mu2 ＝12.546 384 6，Mu3 ＝26.979 046 5， Alpha1 ＝2.965 641 8， Alpha2 ＝4.312 813 73， Alpha3＝1.569 973 01，D1＝0.012 090 53，D2＝D3＝0； 然后模拟装配过程， 保证唇封的过盈量准确， 并通过直接加压法施加流体静压， 网格划分后仿真计算， 得到唇封理论的结构变形和应力应变状态， 进而提取唇封与轴之间的静态接触压力。 为获得影响系数矩阵， 对任一接触区节点施加单位力， 获得所有节点的位移值， 最终得到的位移矩阵即为影响系数矩阵。

1.5 热力学与传热学分析

实际工作状态下， 唇封和轴之间会发生摩擦生热现象， 造成唇口局部温升， 同时密封介质和轴之间进行对流换热， 带走部分热量， 如图2 所示。

图2 唇封系统热行为Fig.2 Thermal behavior of lip sealing system

为考虑温度对唇封力场的影响， 文中在有限元软件ABAQUS 中进行了完全热力耦合分析。 首先在初始流固耦合计算结果的基础上考虑摩擦生热作用得到热载荷大小及作用位置； 然后在ABAQUS 中定义部件的材料属性， 包括随温度变化的热传导率、 膨胀系数及应变能函数； 接着设置对流换热系数及环境温度， 施加表面热通量， 进行温度位移耦合分析， 最终得到唇封温度场及应力应变状态。 新得到的静态接触压力重新代入流固耦合过程计算热载荷， 重复上述热力耦合过程， 直至唇封温度及力学状态不发生变化，完成流固热耦合分析。

接触区域的表面热通量计算公式如下：

式中：Q为表面热通量；Fs为流固耦合计算所得摩擦力；v为轴的线速度；S为密封接触区面积。

摩擦热按照一定比例分配到两个物体上， 即热流分配系数， 计算公式为

式中：Kq为热流分配系数；λs、λr分别为固体、流体热导率；ρs、ρr分别为固体、 流体密度；Cs、Cr分别为固体、 流体比热容。

旋转轴与流体发生强制对流换热行为， 认为周向上轴与周围流体相对速度保持一致， 故可以使用等效外掠平板模型[16］对其换热系数进行计算。

式中：Re＝vL／νr为雷诺数， 其中L为特征尺度，νr为流体运动黏度；Pr＝νrρrCr／λr为普朗特数。

1.6 流固热耦合模型

文中通过上述流体力学、 接触力学、 固体力学和热力学与传热学的相互耦合， 建立了盾构机主驱动唇形密封流固热耦合仿真模型。 图3 展示了所建立的数值仿真模型的具体计算流程。

图3 盾构机主驱动唇形密封流固热耦合数值仿真模型计算流程Fig.3 Calculation flow of numerical simulation model of fluidsolid-heat coupling for main drive lip seal of shield machine

主要步骤如下：

（1） 输入唇封基本参数， 利用有限元软件ABAQUS 进行固体力学分析， 获得静态接触压力和影响系数矩阵；

（2） 根据接触区粗糙峰形貌分布， 定义初始膜厚， 结合有限元仿真结果， 分别对流体压力及粗糙峰接触压力进行计算， 根据耦合关系进行迭代求解， 调整油膜厚度；

（3） 计算表面粗糙峰的切向变形， 更新接触区粗糙峰形貌分布， 直至粗糙峰的切向变形误差达到收敛标准；

（4） 根据摩擦力大小及作用范围， 进行密封系统热力耦合分析， 得到摩擦生热条件下静态接触压力分布； 重新进行上述过程， 直至唇口温度误差满足收敛条件， 流固热耦合计算完成， 输出泄漏率等密封性能参数。

2 结果与讨论

根据提出的流固热耦合仿真模型， 对盾构机主驱动唇形密封进行分析， 涉及的基本参数如表1 所示。

表1 仿真模型输入参数Table 1 Input parameters of simulation model

初始有限元仿真结果如图4 和图5 所示， 分别是静态接触压力分布和影响系数矩阵。

图4 静态接触压力分布Fig.4 Static contact pressure distribution

图5 影响系数矩阵Fig.5 Influence coefficient matrix

数值求解受网格密度、 收敛准则和表面复杂性影响， 根据网格无关化条件， 计算网格设置为97×97；为提高计算结果的准确度， 设置收敛精度为5×10-4，同时利用确定性模型定义初始膜厚， 即膜厚由常数项以及唇和轴表面微变形产生的波动项组成， 膜厚公式如下：

式中：A1为唇口表面粗糙峰半幅值；λ11为唇口表面形貌x方向波长；λ12为唇口表面形貌y方向波长；A2为转轴表面粗糙峰半幅值；λ21为转轴表面形貌x方向波长；λ22为转轴表面形貌y方向波长；h0为膜厚中值。

结合上述基本参数构建的初始膜厚如图6 所示。

图6 初始膜厚分布Fig.6 Initial film thickness distribution

将初始有限元仿真结果和基本参数代入到流固耦合模型中， 可求得无温度场条件下唇封摩擦力等密封性能， 其中摩擦力矩为3.35 N·m， 泄漏率为-2.73 g／mL （代表无泄漏产生）。 利用摩擦力计算唇口热通量， 结合有限元软件ABAQUS 进行热力耦合分析，最终得到的温度场分布， 如图7 所示， 可以看出唇口温度大约为55 ℃。

图7 温度场分布（℃）Fig.7 Temperature field distribution （℃）： （a） overall temperature field of lip seal and shaft； （b） local temperature field of lip

多次迭代后流固热耦合仿真模型精度满足收敛要求， 得到考虑温度场后摩擦力矩、 泄漏率等密封性能参数， 其中摩擦力矩为3.14 N·m， 泄漏率为-2.47 g／mL （即泵送率2.47 g／mL）。 考虑温度场后， 摩擦力矩下降， 泵送率降低， 主要因为唇口温度提高使得材料软化， 导致唇封径向力减小。

图8 所示为考虑温度场前后唇封接触压力分布。 可以看出考虑温度场后， 接触压力最大值由4.06 MPa 减小到3.52 MPa， 但接触宽度由3.55 mm增加到3.61 mm， 考虑温度的唇封接触情况更符合实际。

图8 考虑温度场前后接触压力分布Fig.8 Distribution of static contact pressure before and after considering temperature field

图9 展示了考虑温度场前后唇封的应力应变状态。 可见应力应变均有显著变化， 其中最大Mises应力值由2.3 MPa 减小到2.1 MPa， 最大主应变由0.21 减小到0.19。 最大Mises 应力和最大主应变都集中于唇口部分， 受温度场影响， 该区域材料变软， 是导致最大Mises 应力和最大主应变减小的主要原因。

图9 考虑温度场前后唇封应力应变状态Fig.9 Stress and strain state of lip seal before and after considering temperature field： （a） initial lip seal stress distribution （MPa）；（b） final lip seal stress distribution （MPa）； （c） initial lip seal strain distribution； （d） final lip seal strain distribution

3 结论

基于脂润滑条件下唇形密封流固耦合模型， 结合有限元软件实现温度场与唇封力场的完全耦合， 建立盾构机主驱动唇形密封流固热耦合仿真模型。 主要结论如下：

（1） 考虑温度场后， 唇封的摩擦力矩下降， 泵送率降低， 主要因为唇口温度提高使得材料软化， 导致唇封径向力减小。

（2） 考虑温度场后， 唇封接触压力最大值由4.06 MPa 减小到3.52 MPa， 但接触宽度由3.55 mm增加到3.61mm， 考虑温度的唇封接触情况更符合实际。

（3） 考虑温度场后， 最大Mises 应力值由2.3 MPa 减小到2.1 MPa， 最大主应变由0.21 减小到0.19。 最大Mises 应力和最大主应变都集中于唇口部分， 受温度场影响， 该区域材料变软， 是导致最大Mises 应力和最大主应变减小的主要原因。

（4） 温度对唇封应力应变状态及密封性能产生较大影响， 考虑温度场后的预测结果更贴合实际， 这对盾构机主驱动唇形密封设计具有一定指导作用。