考虑空穴效应的广义特征分解法解雷诺方程

耿煜,郭振,陈渭

(西安交通大学现代设计及转子轴承系统教育部重点实验室,710049,西安)

铰接副作为必要的连接部件在机械系统中广泛应用,由于受到加工误差、材料变形、磨损等因素的影响,铰接副大多存在间隙,导致机构连接处发生碰撞、振动,加快铰接副的磨损,因此间隙铰接副通常依靠润滑将不利影响降到最低。

滑动轴承非线性油膜力模型是含润滑间隙铰接副多体系统动力学特性分析与研究的基础和关键,油膜力的计算精度和速度将直接影响非线性动力学特性分析的效率。目前关于油膜力模型的研究已取得很大进展,涌现出了多种求解非线性油膜力的理论与方法,其中应用较为普遍的为有限差分法(FDM)[1-4]和有限元法(FEM)[5-8]。虽然数值方法精度高,但是计算效率很低,因此运用近似解析法求解油膜力的研究日渐增多,目前滑动轴承非线性油膜力近似解表达式主要有无限短轴承假设下的Capone模型[9]和Sommerfeld提出的无限长轴承油膜力模型[10],但是考虑到油膜的温升和端泄,滑动轴承的宽径比通常在0.8～1.0之间[11],现实中较少存在无限短或无限长轴承,实际应用受到限制,近年来,很多学者针对有限长轴承做了相当多的工作。Chasalevris等基于油膜假设,利用变分原理求出有限长轴承非线性油膜力的近似解析解,但是并没有考虑空穴效应[12-14]。在流体润滑问题中,由于发散间隙的存在会导致油膜出现负压,当油膜达到最大可承受负压时会发生破裂现象,进而产生空穴,所以不考虑空穴效应所获得的油膜力是不可靠的[14]。张永芳等基于下游雷诺边界条件,运用变分原理和分离变量法提出了一种针对有限长轴承的近似解析解方法,计算精度较高,但是计算速度相比传统有限差分法并没有得到很大的提升[11,15]。除此之外,还有诸如数据库法和精确解析法等油膜力计算模型,但是由于通用性不强、计算效率较低等问题,并不适合于一般性径向滑动轴承油膜力的求解。因此,建立考虑空穴效应、通用性较高、计算速度较快的滑动轴承非线性油膜力模型对于提高含润滑间隙铰接副多体系统动力学分析效率和满足实际工程应用至关重要。

广义特征分解法(PGD)主要是通过将求解域转变为张量积函数之和来减少待求未知量个数,从而有效提高计算效率。将待求函数u分解为N维Q个张量积之和的形式,表达式如下

(1)

使用M个自由度对每个张量进行离散,待求未知量有NQM个,而如果使用传统的FEM法或FDM法,待求未知量则为MN个。Amine等第一次使用PGD法对偏微分方程进行求解[16];Dumon等将PGD法成功应用到流体力学领域,并证明了PGD法的高效性[17];Tamellini等使用PGD法对稳态不可压缩Navier-Stokes方程进行求解[18-19];Cherabi等使用PGD法求解稳态工况下的雷诺方程[20],但是并没有考虑空穴对油膜力分布的影响,实际应用具有一定的局限性。

本文针对含润滑间隙铰接副中油膜力的计算问题,在文献[20]的基础上,建立了动载滑动轴承非线性油膜力的PGD计算模型,并引入雷诺边界条件,充分考虑空穴效应对油膜力分布的影响。利用Matlab数值计算验证了模型的高效性和稳定性,明晰了PGD法的适用范围。将模型引入到含润滑间隙铰接副多体系统的计算中去,解决了模型的收敛问题,证明了模型的可靠性和工程应用价值。

1 计算模型的建立

1.1 动载雷诺方程

动载雷诺方程如下式

(2)

可以将雷诺方程转换为以下的形式

(3)

即

(4)

将式(5)代入式(4),得

∂x=R∂θ;U=Rω

(5)

可得到式(4)圆柱坐标系下的表达式

(6)

油膜厚度表达式如下

h=c(1+εcosθ)

(7)

式中:θ为周向坐标;z为轴向坐标;h为油膜厚度;R为轴承半径;μ为润滑油黏度;ω为旋转角速度;c为半径间隙;ε为偏心率。

1.2 PGD求解雷诺方程

将雷诺方程式(6)转换为等效积分的弱形式,若对于任意的权函数p′下式成立

(8)

则式(6)成立。

将油膜力函数p进行分解,获得其近似解

(9)

式中:i为当前维数;Xi和Yi为第i维度下周向和轴向方向向量。一般情况下,维数越多,近似解越逼近精确解。

采用变方向迭代策略,式(9)转换为

(10)

为了方便计算,将权函数p′设为

(11)

将式(10)(11)代入式(8),得

(12)

将式(12)进行整理,得

(13)

其中

将式(13)进一步化简,得

(14)

式(10)、式(11)转换为

(15)

(16)

将式(15)(16)代入式(8),得

(17)

同理,将式(17)进行整理,得

(18)

其中

将式(18)进一步化简,得

(19)

1.3 油膜承载力的求解

采用中差分五点差分形式对式(14)(19)进行求解,将求解得到的Xi(θ)、Yi(z)代入式(9),算出油膜力分布函数p。

轴承承载力计算公式如下

(20)

其中

Fx=∬pRsinθdθdz;Fy=∬pRcosθdθdz

1.4 边界条件

轴承的端泄边、周期边界和雷诺边界的边界条件分别如下

(21)

雷诺边界条件假定在油膜完整区和油膜破裂区的交界线上压力的法向梯度为0,且压力等于空穴压力pc,本文空穴压力设为0。

(22)

1.5 收敛条件

ξ1和ξ2的计算表达式如下

(23)

(24)

1.6 计算程序

计算程序的基本流程如下:

步骤1 输入3个轴承参数,润滑油黏度、最大允许误差λ1、最大允许维度数N=1;

步骤2 划分网格节点,计算油膜厚度分布;

步骤8 由式(9)计算当前油膜力分布,若

步骤9 由式(21)计算误差ξ1,如果ξ1<λ1,则转到步骤10,否则q=q+1,重复步骤4～9。

步骤10 由式(20)计算轴承承载力,输出油膜压力分布。

2 计算模型的验证

2.1 计算精度的验证

Kumar和Booker[4]基于JFO边界条件,采用有限元法对有限长轴承进行了油膜力求解,已有研究[21]表明,在偏心率小于0.02的情况下,雷诺边界条件和JFO边界条件所得到油膜力误差在5%以内;当偏心率高于0.02时,两种边界条件所得结果几乎一致。因此本文使用文献[4]的计算参数和结果,来验证基于雷诺边界条件的PGD法的准确性,具体的轴承参数意义及数值见表1。

表1 文献[4]中的轴承参数

将有限差分法、有限元法和PGD法计算结果与文献[4]的计算结果进行对比,结果如表2所示,文中不作特殊说明FDM和PGD法均采用中差分五点差分格式进行数值求解。采用PGD法得到的油膜承载力和最高压力的结果误差均在10%以内。值得注意的是,随着网格密度的加密,PGD法的计算精度将得到提高,但是相比FDM法计算时间得到了大幅缩减,这将在2.2节进行验证。另外,将使用PGD、FDM、FEM法得出的油膜压力分布图进行比较,如图1所示,可以看到图像曲线具有较好的一致性,再一次验证了基于雷诺边界条件的PGD法的准确性。

表2 PGD、FDM、FEM法与文献[4]方法的结果对比

图1 轴向中心处油膜压力分布图

2.2 计算效率的验证

相比FDM法,FEM法计算精度更高,处理复杂几何结构时适应性强,但是FDM法使用简单,在处理一般性滑动轴承方面应用更为广泛,计算效率相比FEM法有本质上的提高,且其计算精度已经能够满足一般性工程需要,因此在这里将采用新模型的FGD法与FDM法进行比较,验证FGD法的高效性。

设置相同的最大允许误差λ1和SOR迭代因子分别为10-5和1.41,比较PGD和FDM法在不同网格密度下耗费的计算时间。由于PGD法受到初始随机赋值的影响,为了严谨,每种网格密度下的PGD计算时间和承载力取3次平均值。如表3所示,在低网格密度下,PGD与FDM法的相对误差较低,且随着网格密度的增加,误差控制在1%以内,计算效率最大可以提高3.59倍。在网格密度达到300×300以上时,两种方法的计算误差开始突然增大,通过承载力可以看出,FDM法算出的承载力发生明显突变,而PGD法算出的承载力变化较小,这说明在相同的λ1条件下,PGD法具有更好的计算稳定性,随着网格密度的提高,依旧保持了较高的计算精度,而FDM法由于网格密度过高,原有λ1已经无法满足其计算精度的要求,导致结果过早收敛进而使得计算结果发生突变。

为了验证上述结论,在600×600的固定网格密度情况下,采用不同的λ1进行了实验,实验结果如表4所示。在降低最大允许误差后,PGD与FDM法的相对误差始终控制在1%以内,最大加速比达到了4.42,表明PGD法具有相当高的计算精度和效率。结合表2和表3可以看出,与FDM法相比,PGD法不仅保持了较高的计算精度,而且具有更好的计算稳定性与更快的计算速度,在高网格密度和低允许误差条件下,计算优势将会更加明显。

表3 不同网格密度下PGD和FDM法的计算结果对比

表4 不同允许误差下PGD与FDM法的计算时间对比

3 PGD法在含润滑间隙铰接副多体系统中的应用

为了提高含润滑间隙铰接副多体系统的计算效率,符合工业应用所需要的实用性特点,引入模型对含润滑铰接副非线性油膜力进行计算。本文以曲柄滑块为例,曲柄滑块机构见图2,相关参数见表5。

图2 曲柄滑块机构示意图

表5 曲柄滑块机构参数

采用欧拉-拉格朗日法对多体系统进行建模,引入Flores过渡模型来解决润滑间隙铰接副的工作状态过渡问题,具体参数详见表6,其中E与μ表示材料的弹性模量与泊松比。接触力模型与摩擦力模型需要设定恢复系数与摩擦系数,本文分别取为0.9和0.2。Flores模型过渡临界值设为2 μm,润滑油黏度为0.08 Pa·s。具体的公式与计算步骤在此不作赘述,详见文献[1]。

表6 间隙铰接副参数

在多体系统计算中采用FDM法求解雷诺方程时,往往使用上一时刻计算得出的压力分布解作为下一时刻的初始值进行计算加速,使用PGD法时同样使用上一次计算出的周向和轴向向量Xi、Yi来作为初始值进行计算加速,但是Xi、Yi不仅参与迭代,还会影响到迭代系数Aθ～Fθ,Az～Hz的计算,因此连续使用上一次计算解作初始解会产生累计误差,导致算法难以收敛,因此需要在一定数量的时间步长后,对初始向量进行新一轮的随机赋值,以消除累计误差对计算所带来的影响。

曲柄转动3周,设置相同的λ1为10-5,SOR迭代因子为1.41,验证PGD法在含润滑间隙铰接副多体系统计算中的可靠性与实用性。图3为PGD与FDM法在30×30网格密度下的轴心轨迹、最小油膜厚度、滑块速度随时间变化的对比。从图3可以看出,PGD与FEM法的计算结果相差不大,滑块速度曲线近乎重合,可见PGD法具有较高的计算精度。

(a)轴心轨迹图

(b)最小油膜厚度随时间变化

(c)滑块速度随时间变化图3 30×30网格下PGD与FDM法的计算结果对比

表7为两种方法在不同网格密度下运转3个周期的计算时间对比,可以看出PGD法对于计算的加速效果是非常明显的,在30×30和40×40的网格密度下,PGD法平均一个周期的计算时间分别缩短了6 min和20 min,在实际计算中,为了观察间隙铰接副的润滑效果,至少需要进行5个周期的运算,随着机构复杂性的提高,达到稳定润滑状态所需要的周期数会进一步增加,因此采用新模型所节省的时间将是巨大的。

表7 在含润滑间隙铰接副多体系统中的计算时间对比

4 结 论

本文基于广义特征分解法建立了动载有限长滑动轴承非线性油膜力模型,引入雷诺边界条件充分考虑空穴对油膜压力分布的影响,与现有文献数据进行对比,证明了新模型的准确性与高效性,并研究了采用新模型方法的计算精度与计算效率随网格密度的变化趋势。以曲柄滑块机构为例,采用新模型的方法对含润滑间隙铰接副进行动力学分析,证明了新模型的实际应用价值。可以得出以下结论:

(1)PGD法计算编程较为复杂,相较传统计算方法更难理解,但是在较低网格密度下,PGD法与FDM法的承载力计算相对误差最高不超过7%,在较高网格密度下,相对误差控制在1%以内,具有较高的计算精度且稳定性相比FDM法更高。

(2)PGD法计算效率较高,在不损失精度的前提下,在不同网格密度下优势明显,计算速度相比FEM法提高2～5倍。

(3)使用PGD法对润滑间隙铰接副动力学分析,在30×30和40×40网格密度下,单个运算周期节省时间分别高达6 min和20 min,随着网格密度和机械系统复杂性的提高,新模型的计算优势将会更加显著,证明了模型的可靠性与实用价值。

本文只是将PGD法应用于普通雷诺方程的求解,考虑表面粗糙度及温度的影响,例如将PGD法引入平均流模型和热弹流计算中是PGD法未来的发展方向。本文考虑空穴效应的PGD法只是针对普通滑动轴承而言,而对于织构化轴承而言,由于存在大量的收敛、发散间隙,油膜会反复破裂并再形成,这对算法的收敛问题会产生巨大挑战,这是PGD法目前所面临的难题。