考虑油液流变时变特性和深度特征尺寸的静压转台流动分析

金玉峰, 段鸾芳, 吴伟伟

(1.江苏省常州技师学院 智能制造学院, 江苏 常州 213032; 2. 南京工业职业技术大学 电气工程学院, 江苏 南京 210023; 3.扬州大学 机械工程学院, 江苏 扬州 225127)

引言

液压油是数控机床加工过程中不可或缺的材料,主要起到润滑与冷却的作用[1-2]。当液压油用于导轨、丝杠及轴承(如静压导轨、静压丝杠及静压轴承)中时,还需要具有一定的承载能力[3-7]。静压转台是数控机床中的常见部件,其以静压轴承为支承,液压油作为工作介质,在压力作用下形成压力油膜,具有较好的承载刚性[8-9]。

静压转台的结构是影响刚度和承载能力等的主要因素。张瑞涛等[10]利用Fluent软件对节流器孔径及轴承间隙对刚度的影响进行了分析,基于分析结果设计了具有较高刚度的静压转台。仝志伟等[11]分析了支承布局对平面度误差的影响,并对支承布局进行重构,可减小平面度误差,提高转台的承载能力。王京等[12]发现双环形油腔较传统油腔具有更好的承载能力。王永柱等[13]对扇形转台中油腔形状对承载能力的影响进行了分析对比,结果表明口字型油腔具有更好的承载能力。张庆锋[14]分析研究了位移率和间隙对油膜承载能力的影响,根据结果优化了转台结构。申峰等[15]具体分析了封油边形状特征参数对转台承载能力的影响,利用分析结果优化封油边参数可显著改善承载能力。MICHALEC M等[16]利用CFD模型和MATLAB提出一种双参数优化方法,对静压轴承的油腔形状进行了优化,可有效降低能量损失。ZHANG Yanqin等[17]针对双矩形腔的静压转台中载荷对油膜性能的影响展开了研究,研究结果表明,载荷大小与油膜的平均温升之间存在非线性关系。YADAV S K等[18]利用有限元分析对比分析了圆形、椭圆形、方形及矩形油腔对油膜刚度系数的影响,结果表明油腔结构对油膜承载能力具有重要影响。

在静压转台的实际运行工作过程中,通过实验发现液压油的油液品质(油液的黏度)对机床零部件运作时的性能也具有一定影响[19-21]。液压油在使用过程中,随着时间的推移,油液的黏度会发生变化,需要了解该变化对承载能力的影响。本研究将以承载能力较好的双环形油腔结构为研究对象,利用格子玻尔兹曼方法(LBM)分析讨论封油边厚度对油液黏度的时变特性的影响,为后续选择合适的油液黏度和封油边厚度提供理论依据。

1 流变测试

表1 两种油品流变数据拟合结果

图1 两种油品实验数据与拟合结果对比

根据表1结果,将其与流变实验的数据相比较,结果显示出较好的一致性,且表1中的相关系数无限接近于或等于1,因此拟合获得的流变模型是可信的。由表1结果可知,1#和2#油品均呈现出非牛顿特性,具体为带屈服特征的Herschel-Bulkley型流体,仅当剪切应力大于初始屈服应力时,流体才会发生剪切运动。

2 适应Herschel-Bulkley流体的LBM方法

2.1 虚外力项LBM方法

LBM的演化方程为:

(1)

为了解决由非牛顿特性引起的仿真离散及精度差的问题,引入“虚外力项”用以描述Herschel-Bulkley特征,对于一般含外力项的LBM演化方程为[22]:

feq(r,t)]+δtF′

(2)

式(2)中包含了LBM的两个主要过程——粒子迁移和粒子碰撞,其中粒子迁移具体可描述为:

f(r+eiδt,t+δt)=f+(r,t)

(3)

粒子碰撞具体可描述为:

(4)

式中,r—— 位置矢量

t—— 某时刻

f—— 粒子分布函数

feq—— 平衡态分布函数

其具体表达式为:

(5)

式中,ρ—— 格子密度,通常取1

u—— 速度矢量

ωi—— 权重系数

对于LBM的D2Q9模型,其具体取值如下:

(6)

其中,ei为D2Q9模型的速度配置,具体取值为:

(7)

(8)

黏度一般可根据流变模型计算获得,对于牛顿流体,黏度不随剪切率的变化而变化,其为恒定常数,则松弛时间也为恒定值;而对于非牛顿流体,黏度随剪切率的变化而变化,因此松弛时间为一变化值。对于液压油所属的Herschel-Bulkley流体,其流变模型为:

(9)

式中,τh—— Herschel-Bulkley流体的剪切应力

τh0—— 初始屈服应力

μh0—— 黏度系数

n—— 幂律指数

为了方便计算,前人针对Herschel-Bulkley流体提出了修正黏度-剪切率方程为:

(10)

(11)

其中,DII为应变率张量第二不变量,其可通过下式获得:

(12)

式中,l—— 计算时的维数,此处l=2

Sαβ—— 应变率张量

其具体计算公式为:

(13)

对于式(4)中的虚外力项,具体计算方式如下:

(14)

(15)

结合式(14)和式(15),可知离散虚外力项具体可以描述为:

(16)

2.2 理论验证

为了验证上述提出的虚外力项方法的可靠性,利用Herschel-Bulkley流体在泊肃叶流中的理论解进行验证,其理论解为[23]:

(17)

式中,H为泊肃叶流中两平行平板之间的距离,yτ为一临界点,表征了剪切应力是否达到屈服应力,具体值为yτ=τh0/(∂p/∂x),∂p/∂x为在x方向上的压力梯度。具体参数设置如下:H=1,平板长度L=1,压力梯度∂p/∂x=-4×10-3,黏度系数μh0=0.1,幂律指数n=0.9,初始屈服应力τh0=6×10-4。在数值模拟计算时,格子数取为200×200,此处仅为了验证模型的可靠性,因此直接使用了无量纲数进行模拟与对比。经过仿真,对比结果如图2所示,对比结果表明提出的虚外力项LBM具有较好的可靠性,可用于Herschel-Bulkley流体的仿真中。

图2 泊肃叶流LBM数值解与理论解的对比

2.3 静压转台仿真步骤

将上述提出的虚外力项LBM用于液压油在静压转台中的流动分析,其具体过程如下。

(1) 初始化参数。设定初始参数,主要包括格子数(由转台的尺寸决定)、入口压力或入口速度、格子密度、初始粒子分布函数、初始黏度、初始松弛时间等。在此步骤中,需进行单位转换,LBM是一种无量纲方法,此处以雷诺数为准则数,由实际尺寸和格子大小可确定长度比例,自定义一较小格子速度,由实际速度值和格子速度可确定速度比例,基于式(18),保证实际雷诺数与格子雷诺数相等,即可获得运动黏度比例,以此类推可获得其他比例,进而计算出其他格子物理量,完成实际物理量与仿真格子参数之间的转换[24]:

(18)

其中,v为特征速度,d为特征长度,υ为运动黏度。

(2) 初始化平衡态分布函数。在初始阶段可以自定义一个函数作为初始平衡态分布函数,后续根据式(5)实际计算迭代获得新的平衡态分布函数。

(3) 执行碰撞步式(4)和迁移步式(5)。在执行碰撞步时,需根据式(13)～式(16)计算虚外力项,用以表征Herschel-Bulkley流体的非牛顿特性。

(4) 计算松弛时间。由式(13)计算应变率张量,结合式(11)、式(12)计算剪切率,由式(8)和式(10)计算下一次迭代所需的松弛时间。

(5) 边界处理。不同于传统有限元方法的边界处理,LBM的边界处理主要是对边界粒子的状态进行定义,目前比较成熟的方法有反弹格式、非平衡外推等,此处选用非平衡外推进行边界处理。

(6) 判断是否计算结束。设定一较小值为对照标准,以某一物理量为基准,计算上一循环与本次循环的结果差值的绝对值,当此值小于上述设定的对照标准时,即认为计算已满足条件,跳出迭代计算,否则进行下一轮迭代循环计算。

(7) 计算宏观物理量。参考式(19)进行物理量计算,并再次进行量纲转换获得宏观物理量。

(19)

式中,ρ′ —— 仿真中油液密度

u—— 仿真中获得的速度

3 不同静压转台结构中的流动分析

3.1 静压转台结构

基于双环形静压转台较传统静压转台具有更好的承载能力,采用如图3所示的双环形油腔结构的静压转台为研究对象。

图3 双环形静压转台结构

油液入口1的半径r1设置为4 mm,油液入口2和3的半径r2均设置为2 mm,内圆环的内圆半径R2设置为20 mm,外圆半径R3设置为24 mm,外圆环的内圆半径R4设置为46 mm,外圆半径R5设置为50 mm,油腔的总深度为H,封油边深度为h,为了研究油液黏度时变特性下深度对油液在腔体内部流动的影响,其具体设置如表2所示。由于双环形油腔是一典型的对称结构,为提高运算效率,在以下的仿真中均取一半结构进行分析,仿真时径向的格子数设置为1000,轴向的格子数由具体尺寸等比例设定。

表2 油腔深度的设置

3.2 1#油品流动分析

根据2.1节及2.3节中的描述,结合MATLAB 2017b软件将1#油品的流变方程代入LBM方法中进行仿真分析,设定油液在3个入口的速度均为0.1 m/s,可以得到如图4～图7的结果,其中图4为速度云图,图5为顶层油膜速度大小的分布情形,图6为顶层油膜速度矢量的角度情形,图7为顶层油膜的压力分布结果。通过对比图4a～图4c可知,在保持油腔深度h2(h2=H-h)不变的情形下,不断增加h,速度云图的变化不明显。根据图4d～图4f可以发现,在保持h不变的情形下,不断增加h2,速度云图差异较为明显,在入油口1处,由1个主涡逐渐发展成2个主涡相连,在入油口2处,主涡强度逐渐增强。在环形台阶上方,速度增强的现象逐渐减弱。

图4 1#油品在6种结构参数下的速度云图

图5 1#油品顶层油膜的速度大小分布

图6 1#油品顶层油膜的速度矢量角度

图7 1#油品顶层油膜的压力分布

图5a和图5b分别描述了Case 1～Case 3和Case 4～Case 6情形下顶层油膜的速度大小分布情况,两图具有相似的变化趋势,在两个环形台阶上方的顶层油膜具有较大的速度,且入油口2输入的油液增强了其左右两侧环形台阶上方顶层油膜的速度,造成同一台阶上方的速度大小有较明显的不同。图5b表明,相同的h下,较小的h2会导致环形台阶上方出现较大的速度,与图4d～图4f展现出的速度云图相吻合。图6中速度矢量与水平线的夹角展示了近乎一致的变化,该角度的变化与形成的涡流在顶层处的形状相关,由图4可看出在顶层油膜涡流具有相似形状,与图6量化的角度大小基本吻合。

图7的压力分布结果间接展示了转台的承载能力,6种情形的结果均表明,中心凹槽的压力远大于环形凹槽的压力,间接说明了入油口2处油液的输入形成了外围涡流大大增加了中心凹槽的压力,起到了增压的作用,提升了转台的承载能力。此外,通过对比发现H值越小,顶层油膜所受的压力越大,且呈现指数级增长的趋势,中心凹槽的压力与环形凹槽的压力差值也越大。

3.3 2#油品流动分析

图8～图11为2#油品仿真分析的结果,图8为速度云图,其中图8a～图8c中,保证h2不变,逐渐增大封油边h,速度云图未有明显差异。图8d～图8f中,保证h不变,h2较小时,中心凹槽内的涡流较明显,随着h2逐渐增大,中心凹槽内的涡流逐渐减弱,但均比图8a～图8c中的涡流明显。

图8 2#油品在6种结构参数下的速度云图

图9 2#油品顶层油膜的速度大小分布

图10 2#油品顶层油膜的速度矢量角度

图11 2#油品顶层油膜的压力分布

图9中展示了顶层油膜速度大小的分布情况,图9a表明随着h的逐渐增大,环形台阶上方的速度逐渐减小,图9b表明随着h2的逐渐增大,顶层油膜的速度在径向分布上均有一定程度的减小,在环形台阶上方的速度减幅较大。图10中所展示的速度矢量角度表明在6种情形下顶层油膜速度与水平线夹角的变化趋势几乎一致,与图8所示的速度云图基本吻合。

图10中展示了顶层油膜受压情况,由图可知,入油口2中油液的输入同样增加了中心凹槽内的压力,起到了增压作用,且Case 4～Case 6的结果表明,H较小时,压力的增幅较明显。

通过综合对比1#油品与2#油品的仿真结果可以发现,随着时间的推移,油液中混入的杂质等改变了油液的流变特性,增加了油液的黏度,对比图3与图7、图4与图8可知,油液黏度的增加使油液的流速明显减小,形成的主涡强度较弱。对比图6与图10的压力分布可知,油膜的承压能力也大大减弱。

4 结论

经过实际加工过程观察及流变实验分析发现,液压油在经过一段时间的使用后,油液的流变特性发生了明显变化,且刚使用与使用一段时间后的油液均呈现出典型的非牛顿特性,针对非牛顿油液流变随时间变化的特点,提出了一种虚外力项LBM方法,对油液在双环形转台的内部流动进行分析,通过仿真结果对比后,可以得出以下结论:

(1) 双环形转台的封油边深度h变化时,速度云图及顶部油膜速度变化主要集中在环形台阶上方,但压力分布差异较明显,随着封油边深度h的增加,入油口2油液的输入对中心凹槽的增压作用减弱;

(2) 双环形转台的油腔深度h2不断增大时,两种油品的顶层油膜速度沿径向分布减小,且在环形台阶上方的速度减幅较大。两种油品中心凹槽内的速度云图有明显不同,承载能力呈现“断崖式”减小;

(3) 液压油的流变特性随时间的推移发生变化,黏度不断增大,减弱了油液的流动,导致2#油品情形下,顶层油膜的速度减小,承载能力也大幅减弱。对于刚使用的1#油品,随着油腔深度增加,中心凹槽的主涡区域不断扩大,而2#油品的中心凹槽主涡区域不断减小,此现象也是由油液黏度的变化导致。为保证转台的正常工作,应定期更换转台用液压油。