一种机械密封润滑液膜空化预测方法.pdf

一种机械密封润滑液膜空化的预测方法，是一种用于润滑液膜空化模拟且遵循质量守恒的预测方法，采用流线迎风SUPG有限元法求解雷诺方程并进行高效的数值迭代。首先分别计算液膜压力p和液膜密度比θ的刚度矩阵Kp和Kθ，然后引入压力和密度的统一综合变量Φ，并对原有线性方程进行修改，最后利用数值迭代技术求解上述获得的线性互补方程，最终同时求得润滑液膜压力p和密度θ。本发明为机械密封等润滑液膜空化问题的数值模拟提供了一种高效精确的预测方法，能够对具有复杂动压槽的元件性能进行准确的预测，相比于其它方法，本发明采用的方法具有计算速度更快、计算精度更高、满足液膜质量守恒要求的优点。

1.  一种机械密封润滑液膜空化预测方法，其特征在于：
所述的预测方法是将机械密封润滑液膜空化问题的所涉及的液膜压力和液膜密度两未知变量综合到一个偏微分方程中，对该方程采用流线迎风SUPG有限元法进行处理，将人工粘度沿速度方向引入偏微分方程中的对流项中，得到其弱积分形式，引入润滑液膜压力和密度比的统一综合变量，形成关于综合变量的互补线性方程组，采用数值迭代技术获得综合变量的解，最终得到端面的压力分布和密度比分布，则密度比θ<1的区域为空化区，而密度比θ＝1的区域为液膜区，从而实现了空化问题的处理；
详细包括如下步骤：
(1)对计算域进行有限元网格划分；
(2)计算刚度矩阵K^p和K^θ，设置初始条件k＝0时(k为迭代次数)，F^k＝1：刚度矩阵K^p和K^θ采用流线迎风格式获得，其计算表达式分别为
Kijp=&Integral;Ωh36μ(&PartialD;Ni&PartialD;x&PartialD;Nj&PartialD;x+&PartialD;Ni&PartialD;y&PartialD;Nj&PartialD;y)dΩKijθ=&Integral;ΩhNi(U&PartialD;Nj&PartialD;x+V&PartialD;Nj&PartialD;y)dΩ-&Integral;Ω12τSUPGh(U&PartialD;Ni&PartialD;x+V&PartialD;Ni&PartialD;y)(U-&PartialD;(hNj)&PartialD;x+V&PartialD;(hNj)&PartialD;y)dΩ]]>
其中，N为有限元形函数，h是液膜厚度，μ是润滑液动力粘度，U、V是滑动表面沿x、y坐标轴的分速度，是稳定系数并且h_T是在速度U(U_x,U_y)方向上有限元网格的特征长度，Ω为积分单元；
(3)计算得到矩阵A＝K^pC-K^θ(I-C)；
其中，F是数值为1或0的空化指示向量(1代表液膜区，0代表空化区)；I为单位矩阵；C为一个对角线元素为空化指示向量F的稀疏矩阵；
(4)采用GMRES迭代方法求解线性方程组AΦ＝K^θF，得到统一综合变量Φ，并且初始计算时假设润滑液膜全部为非空化区，采用块迭代技术更新和检查统一变量Φ和空化指示向量F的一致性；
(5)采用如下规则检查并更新F^k：如果Φ_i<0和F_i^k＝1，则令F_i^k+1＝0；如果Φ_i>1和F_i^k＝0，则令F_i^k+1＝1；
(6)如果F_i^k+1-F_i^k＝0，则迭代终止，否则k＝k+1，重复步骤(3)～(6)；
(7)润滑液膜压力p和密度θ通过如下表达式求解获得:
p＝FΦ
θ＝(1-F)Φ+F
由此同时获得了机械密封润滑液膜压力p和密度比θ；润滑液膜空化区和非空化区可通过密度比θ的大小加以区分：密度比θ<1的润滑区域为空化区，密度比θ＝1的区域为液膜完整区；机械密封的相关密封性能则可由获得液膜压力进行计算获取。

一种机械密封润滑液膜空化预测方法
技术领域
本发明涉及旋转机械中的机械密封领域，特别是涉及一种机械密封润滑液膜空化问题的预测方法。
背景领域
液膜润滑机械密封广泛应用于旋转机械中，为提高承载能力和减小相对运动界面的摩擦磨损，往往在一个或两个运动界面上通常加工一些复杂的动压槽。正是由于复杂动压槽的存在以及泊肃叶流和库埃特流产生的空化效应，使得难以对润滑液膜流体动压效应进行准确的预测。造成采用传统方法求解润滑控制方程-Reynolds方程失败的主要原因有三：1.质量守恒条件下的空化效应；2.在高转速下(高Peclet数)，椭圆—双曲线型方程的对流占优带来的求解稳定性问题；3.复杂的动压槽结构导致液膜厚度的不连续和高的压力梯度。
如今，针对空化条件下雷诺方程的求解问题，已经提出了很多的数学模型和数值算法。其中，JFO空化边界条件是广泛应用的液膜质量守恒的边界条件其相应的数学模型为“ρ-θ”模型，其中JFO空化边界条件与方程组成为互补方程。Ausas等提出了一种有效的松弛方法来对满足互补性条件的解进行迭代更新。由于有限单元法能够建立复杂的边界，Kumar和Booker以及Boedo和Booker采用有限元法在JFO理论条件下对空化区域进行了预测，对液膜破裂和改善边界进行了区分。Hajjam和Bonneau提出了一种改进的ρ-θ有限元算法，在全膜区用高阶形函数，在空化区用低阶形函数来模拟唇形密封和瞬态滑动轴承。Schweizer采用拉格朗日-欧拉方程处理雷诺方程和JFO空化边界条件，并且为一些数值算例提供了有限元方程。但是这些算法普遍存在计算量大、算法实现较为复杂的缺点，在很大程度上影响了对润滑液膜空化问题的研究和发展。
发明内容
为了克服现有方法中空化边界条件液膜泄漏量不守恒、对具有复杂动压槽的液体润滑机械密封的性能无法进行准确预测以及计算量大、计算精度低、收敛速度慢等缺点，本发明提供了一种机械密封润滑液膜空化预测方法，该方法具有更快的计算速度、更高的计算精度，能满足液膜质量守恒的要求。
本发明的技术方案：
一种机械密封润滑液膜空化预测方法，其特征在于：
所述的预测方法是将机械密封润滑液膜空化问题的所涉及的液膜压力和液膜密度两未知变量综合到一个偏微分方程中，对该方程采用流线迎风SUPG有限元法进行处理，将人工粘度沿速度方向引入偏微分方程中的对流项中，得到其弱积分形式，引入润滑液膜压力和密度比的统一综合变量，形成关于综合变量的互补线性方程组，采用数值迭代技术获得综合变量的解，最终得到端面的压力分布和密度比分布，则密度比θ<1的区域为空化区，而密度比θ＝1的区域为液膜区，从而实现了空化问题的处理；
详细包括如下步骤：
(1)对计算域进行有限元网格划分；
(2)计算刚度矩阵K^p和K^θ，设置初始条件k＝0时(k为迭代次数)，F^k＝1：刚度矩阵K^p和K^θ采用流线迎风格式获得，其计算表达式分别为
Kijp=&Integral;Ωh36μ(&PartialD;Ni&PartialD;x&PartialD;Nj&PartialD;x+&PartialD;Ni&PartialD;y&PartialD;Nj&PartialD;y)dΩKijθ=&Integral;ΩhNi(U&PartialD;Nj&PartialD;x+V&PartialD;Nj&PartialD;y)dΩ-&Integral;Ω12τSUPGh(U&PartialD;Ni&PartialD;x+V&PartialD;Ni&PartialD;y)(U-&PartialD;(hNj)&PartialD;x+V&PartialD;(hNj)&PartialD;y)dΩ]]>
其中，N为有限元形函数，h是液膜厚度，μ是润滑液动力粘度，U、V是滑动表面沿x、y坐标轴的分速度，是稳定系数并且h_T是在速度U(U_x,U_y)方向上有限元网格的特征长度；Ω为积分单元。
(3)计算得到矩阵A＝K^pC-K^θ(I-C)；
其中，F是数值为1或0的空化指示向量(1代表液膜区，0代表空化区)；I为单位矩阵；C为一个对角线元素为空化指示向量F的稀疏矩阵；
(4)采用GMRES迭代方法求解线性方程组AΦ＝K^θF，得到统一综合变量Φ，并且初始计算时假设润滑液膜全部为非空化区，采用块迭代技术更新和检查统一变量Φ和空化指示向量F的一致性；
(5)采用如下规则检查并更新F^k：如果Φ_i<0和F_i^k＝1，则令F_i^k+1＝0；如果Φ_i>1和F_i^k＝0，则令F_i^k+1＝1；
(6)如果F_i^k+1-F_i^k＝0，则迭代终止，否则k＝k+1，重复步骤(3)～(6)；
(7)润滑液膜压力p和密度θ通过如下表达式求解获得
p＝FΦ
θ＝(1-F)Φ+F
由此同时获得了机械密封润滑液膜压力p和密度比θ；润滑液膜空化区和非空化区可通过密度比θ的大小加以区分：密度比θ<1的润滑区域为空化区，密度比θ＝1的区域为液膜完整区；机械密封的相关密封性能则可由获得液膜压力进行计算获取。
本发明所述的刚度矩阵K^p和K^θ用高斯积分法计算且只需计算一次，不需要参与循环迭代计算，因此计算效率较高。
本发明的有益效果主要表现在：
1.本发明采用流线迎风SUPG有限元法来求解雷诺方程，保证了解的稳定性，为润滑液膜空化问题的数值模拟提供了一种高效精确的计算方法。
2.用高斯积分法计算刚度矩阵K^p和K^θ且只需计算一次，不参与循环迭代计算，只有稀疏矩阵A在循环内部更新，因此大大减少了计算量。
3.所采用的有限单元网格可避免处理润滑液膜厚度不连续，因此解决了由此带来的计算误差问题。
4.相比于其它方法，本发明采用的算法具有计算速度更快、计算精度更高、能满足质量守恒要求的优点。
附图说明
图1为本发明的流程示意图。
图2为本发明采用的算法与其他两种算法的性能参数比较图。
具体实施方式
以下实施例用于说明本专利，但不用来限制本专利的范围。现通过附图和实施例对本发明做进一步的详细描述。
实施例1
图1为本发明实施例1提供的一种机械密封润滑液膜空化预测方法的流程示意图，本实施例主要包含以下步骤：
步骤1：对计算域进行有限元网格划分；
步骤2：在上述有限元网格上采用高斯积分法计算如下积分式，经整体组装后获得刚度矩阵K^p和K^θ
Kijp=&Integral;Ωh36μ(&PartialD;Ni&PartialD;x&PartialD;Nj&PartialD;x+&PartialD;Ni&PartialD;y&PartialD;Nj&PartialD;y)dΩKijθ=&Integral;ΩhNi(U&PartialD;Nj&PartialD;x+V&PartialD;Nj&PartialD;y)dΩ-&Integral;Ω12τSUPGh(U&PartialD;Ni&PartialD;x+V&PartialD;Ni&PartialD;y)(U-&PartialD;(hNj)&PartialD;x+V&PartialD;(hNj)&PartialD;y)dΩ]]>
步骤3：令k＝0，F^k＝1；
步骤4：通过下式计算稀疏矩阵A：
A＝K^pC-K^θ(I-C)
步骤5：用GMRES迭代方法求解如下线性方程，获得统一变量Φ；
AΦ＝K^θF
步骤6：检查并更新F^k：如果Φ_i<0和F_i^k＝1，则令F_i^k+1＝0；如果Φ_i>1和F_i^k＝0，则令F_i^k+1＝1；
步骤7：如果F_i^k+1-F_i^k＝0，则迭代终止，否则k＝k+1，重复步骤(4)～(7)；
步骤8：得到润滑液膜压力p和密度θ,
p＝FΦ，θ＝(1-F)Φ+F
实施例二
图2为本发明实施例二提供的本预测方法与采用JFO空化边界条件的有限差分法(JFO-FDM)、采用雷诺空化边界条件的有限单元法(REYNOLDS-FEM)的比较，参数包括时间t、开启力F_o、以及内外径处的泄漏率Q_i和Q_o。本预测方法得到的内外径 的泄漏量几乎相等，满足质量守恒要求；得到的开启力最准确并且所用的计算时间最短，计算效率最高。
通过实际计算证明，本发明得到的结果更加准确、可靠，验证了本发明的合理性和准确性。
本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。