基于层次模型的雷诺方程自适应变步长差分算法研究*

2019-08-02

润滑与密封 2019年7期

(1.武汉科技大学冶金装备及其控制教育部重点实验室湖北武汉 430081；2.武汉科技大学机械传动与制造工程湖北省重点实验室湖北武汉 430081)

雷诺方程实际上是由不可压缩黏性流体Navier-Stokes动量方程和连续性方程变形而得到的一个关于压力的二阶偏微分方程。在诸多工程问题中，流体润滑问题的求解首先需要对雷诺方程进行求解[1-2]。当今，求解雷诺方程使用最为普遍的是有限差分法[3]，其通常用等分网格对整个模型区域进行离散，并采用差商代替导数对方程进行离散求解，计算机比较容易实现。对雷诺方程差分求解而言，网格的质量对计算的结果有重要影响，王羽[4]在网格单元与节点相邻的半步长上插入若干节点，在提高计算精度的同时也增大了计算量。针对无织构轴承，于如飞和陈渭[5]用有限差分法对雷诺方程进行了网格无关性分析，精度会有所提高，但与低精度相差甚微，且计算机耗时较高。采用传统有限差分法求解雷诺方程，随着步长的减小会大大增加计算量，对计算机要求较高并且浪费资源。

针对传统有限差分法求解雷诺方程的限制，BERGERA和OLIGERB提出了网格自适应的思想[6]，即将网格点进行最优分布，对物理解变化较大的地方用小步长的网格，对需要重点模拟的区域也使用小步长网格，从而使网格步长根据待解物理量的特性和对应模型的形状而改变。JASTRAM和TESSMER[7]将该思想运用到弹性波数值模拟中，提出了纵向网格步长逐渐变化的算法。刘志强等[8]将基于变分原理的自适应网格生成技术，引入到了起伏海底速度模型的网格剖分中。宓铁良等[9]采用自适应网格，利用小波函数的多尺度分解的特性模拟了地震波的传播，该方法通过判定在初始网格点处的小波系数与阈值之间的大小决定是否对网格细化。

本文作者基于早期学者的研究成果，提出了一种基于网格层次化思想的雷诺方程差分算子对网格进行优化的方法，并在此基础开发了以压力梯度变化与区域积分值为基准的自适应细化区域算法；以活塞-缸套润滑系统和径向滑动轴承为研究对象，详细研究对比了不同算法下不同自由度对求解结果的影响，验证了基于层次化思想的自适应变步长差分算法的优越性及适应性。

1 有限差分法求解雷诺方程

1.1 雷诺方程的形式及边界条件

雷诺方程的一般形式是三维的，为了方便求解，根据不同的假设条件可以将雷诺方程简化为二维形式和一维形式，在相关理论下都可以得到较精确的雷诺方程的解析解，但一维形式的假设往往不符合实际的工程设计，其求解结果必然与实际的工况有较大的出入。文中选取工程中最常见的雷诺方程的二维简化形式进行数值计算。以径向滑动轴承为例，由Navier-Stokes方程和质量连续性方程导出三维雷诺方程[10-12]的一般形式，如下式：

(1)

上式右边三项依次称为楔形项、伸张项及挤压模项。对于一般径向轴承往往是轴承固定，轴颈转动，且外载荷为静载荷。即u1=0,u2=u。另外，v2-v1可以用油膜厚度及其导数的形式来表示：

(2)

根据上述分析可以得到径向滑动轴承雷诺方程的二维形式：

(3)

将求解域D归一化到区间[0,1]×[0,1]后得到

(4)

其他形式的雷诺方程也可以用上述形式进行表示，如缸套活塞-缸套模型对三维雷诺方程进行降维、变量替换和归一化得到雷诺方程后相对应的系数项和右边项：

求解域D的边界为Γ=ΓD∪ΓN∪ΓR，其中ΓD为Dirichlet边界条件，ΓN为Neumann边界条件,ΓR为Robin边界条件。

1.2 变步长差分格式及迭代矩阵的建立

将所求的油膜润滑区域网格化，形成若干单元区域，如图1所示，用各个节点上的压力值p构成一阶、二阶差商来逼近导数，将二维偏微分方程和其边界条件转化成节点的差分方程。非均匀网格划分下的油膜区域在X方向用j来表示列数，Y方向用i来表示行数。图中总共有m×n个节点，其中在X方向有n个节点，在Y方向上有m个节点，并且任意节点的一阶、二阶导数可以由其周围节点的压力值进行近似表示：

得到最终的差分格式为

pi,j=kr·pi+1,j+kl·pi-1,j+ku·pi,j+1+kd·pi,j-1-ff

(5)

式中各系数值随节点位置的改变而改变，其中：

图1 油膜网格划分

非均匀网格划分的求解域总共有m×n个节点，编码顺序为从左往右、自下而上，依次为1,2,3,......,m×n。方程(5)是有限差分法的计算方程，在离散的网格化区域内，每一个节点都对应一个差分方程，根据边界条件确定差分格式的各项系数并通过编码的形式建立系数迭代矩阵。这样便于调试程序时迭代矩阵系数项的修改，减小了编程的工作量并且保证了程序的稳定性，适应于不同模型的系数迭代矩阵的生成。其中第一类边界的加载是索引处差分格式系数所在的编码位置1，其余置0。第二类边界条件加载将相应的差分格式的系数替换成原来的两倍。

2 自适应变步长差分算法

2.1 自适应变步长网格思想

网格划分的质量对雷诺方程数值求解精度和效率起着至关重要的作用，人们在网格的构造过程中往往倾向于选择大小均匀，长宽适当、相互正交的网格；为了获得高精度的计算结果往往需要大规模的网格来离散整个模型，这种方法计算灵活性较差，并造成了计算资源的浪费。为了改善这一弊端，将自适应变步长网格思想引入到雷诺方程的差分计算模型中。自适应变步长网格方法将网格点进行了最优分布，做到了网格点的分布与物理解的耦合，从而最终达到在物理解变动较大的区域网格自动密集，在物理解变化平缓区域网格相对稀疏的目的，摒弃了传统有限差分法所使用的均分网格的做法，在保证了精度的同时又提高了计算速率。

2.2 自适应变步长网格构造过程

自适应变步长网格[13]构造过程中最主要的问题便是根据网格品质或者数值计算的特点来确定自适应细分区域。以压力梯度变化与单元区域积分值为基准的细化区域算法来确定自适应细分区域，并通过归一化后灰度图来对细分区域进行直观的表达。自适应细分区域采用矩形网格进行局部的细化，这样要比采用非构造网格计算速度快。

采用矩形网格对求解域进行均匀划分，步长为h0，这种矩形网格定义为基网格。所有细分的网格区域就是由基网格确定。基网格用Ω0表示。Ω0的下标0表示尺度。根据划分的网格的粗细来定义网格的尺度。

使用有限差分法进行数值计算得到基网格上的数值解。首先以一阶密度梯度为判定准则来进行网格细化区域的确定。如图1计算i,j点的一阶密度梯度，过程如下：

(6)

将求解域中的网格节点的一阶密度梯度进行归一化处理。以图1中节点pij为中心，将网格节点的一阶梯度均分给周围四块网格区域，计算出每一个网格点仅仅在梯度方面对周围的四块单元网格区域贡献值。由于每一个网格区域含有4个网格点，则这4个网格点的梯度密度都会影响到自适应细分区域的确定。将每个网格区域的4个节点的贡献值叠加可以得到单元网格区域的一阶密度梯度。将网格节点梯度的变化转移到网格区域上以此来更好地确定细分区域。然后将满足下式中的网格区域过滤出来：

W≥ηWmax

(7)

其中η根据具体的求解和需要的效果确定。

按照上述步骤过滤出来的网格细分区域分布可能比较集中，也可能比较分散，细分区域分别用Ω1来表示，则Ω1可以表示为

Ω1=∪Ω1，j

(8)

j表示需要进行网格细分的区域。例如在判定准则1下的网格灰度图如图2所示。

图2 一阶梯度密度灰度图

然后通过单元网格的油膜压力积分来二次确定细化区域。通常采用高斯积分法[14]对0尺度各个网格的油膜进行积分计算，得到细分区域的油膜压力积分S0。

S0=∬p(x,y)cosφdA

(9)

其中p(x,y)为积分点的函数值。

求解时候在每个网格单元上选择一些点作为积分点，求出其压力函数在网格节点上的数值，最后对这些点进行加权求和，求出积分点的油膜压力函数值，二维高斯积分的表达式如下：

∬f(x,y)dA=∑∑HiHjf(xi,yi)

(10)

其中积分点xi、yi和权函数Hi、Hj按相应的算法选取。

这样将网格细分的区域弱化成点，对单元网格的油膜压力积分再进行归一化对比后排序，就可以清晰明了地看到每块网格的压力积分的大小，从而将满足式(11)中的网格区域过滤出来。

S1≥γSmax

(11)

其中γ根据具体的求解和需要的效果确定。

按照油膜压力积分过滤出来的网格细分区域用Ω2来表示，Ω2来表示需要进行网格细分的区域。则Ω2可以表示为

Ω2=∪Ω2，j

(12)

j表示在判定准则2下确定网格细分的各个区域。例如在判定准则2下的灰度图如图3所示。

图3 单元网格压力积分灰度图

判定准则1是通过将点的一阶梯度密度映射到网格区域范围内，而判定准则2是将油膜区域弱化成点比较油膜压力积分的大小，各自划分的Ω1和Ω2区域可能集中也可能比较分散。

假设1：当区域都比较集中时，则最后细分的区域为

Ω=Ω1∪Ω2

(13)

假设2：当Ω1区域较集中时，则以Ω1区域为主要参考区域，根据Ω2区域的分布进行适当的扩展,由此得到最终的细分区域Ω。

假设3：当区域都比较分散时，则将细分区域进行分开处理，但数值实验证明这种情况很少。

进行变步长计算时，将区域网格按照假设区域的灰度值大小(假设1就按照各自的区域灰度值排序，假设2就按照Ω1区域灰度值进行排序)分为3块区域，步长依次为h1,h2,h3，h1=2h2=3h3=1/3h0，随着细分的进行将步长缩小为原来的1/2。

3 自适应变步长差分算法求解流程

在MATLAB 中开发出求解程序，其求解流程图如图4所示，步骤如下：

(1)输入方程的常值参数，即参数设置，包括活塞裙部半径、缸套半径、活塞裙部长度、活塞裙部上下端半径以及润滑油黏度等工况参数。

(2)初始化网格，按照图1所示的网格示意图进行编码，根据公式(5)求出系数迭代矩阵编码位置对应的系数，从而生成系数迭代矩阵和右边项，计算出初始网格的油膜压力。

(3)计算单元网格的油膜压力一阶密度梯度和油膜压力积分，根据自适应网格构造过程中的准则1和准则2来确定是否存在细分区域。若存在细分区域，继续细分生成第2层网格，然后按照准则1和准则2进一步确定细分区域，细化网格，最终得出自适应变步长网格下的系数矩阵和右边项。

(4)进行迭代计算得到油膜压力p的数值，并与初始化网格的油膜压力进行对比，若符合条件就输出，不符合条件就调整初始网格。迭代可以采用SOR迭代，要注意松弛因子的选取。

(5)设定如下的收敛判定准则：

(14)

式中：σ即为允许的相对误差值，取为1×10-6。

4 数值算例

4.1 内燃机活塞-缸套润滑模型

活塞缸套之间的膜润滑主要涉及通过求解雷诺方程得到油膜压力。图5(a)所示为活塞-缸套系统，其油膜润滑区域如图5(b)所示[15]。

因为左右润滑区域相对于平面φ=0都是对称的，所以只需考虑前面部分。通常，当活塞周期性地移动时，由于负载变化，密度和黏度的大小会发生变化。但为了简化数值模型，这里将密度和黏度的值作为常数。求解得出油膜压力的计算参数如表1所示。

图5 活塞-缸套润滑系统示意图

参数说明值单位r活塞裙部半径21.725×10-3mR缸套半径21.75×10-3mLsk活塞裙部长度22.5×10-3mη润滑油黏度0.012 95Pa·set活塞裙部上端偏心距0.785 65×10-5meb活塞裙部下端偏心距-0.249 24×10-4mdetet对时间的导数-0.243 5×10-2m·s-1debeb对时间的导数-0.844 61×10-2m·s-1v活塞往复运动速度-10.654 3m·s-1Cb活塞销到活塞冠部垂直距离6×10-3mCc活塞销到活塞中心线距离0m

假设v0=0,vh=v(其正方向为上行)，润滑剂无侧漏，因此u0=uh=0。那么有：

(15)

油膜厚度h(φ,y)采用椭圆法进行计算：

h(φ,y)=R-

将域D归一化为[0,1]×[0,1]，进行如下转换：右侧(因为对称，只考虑前面部分)

边界条件为

通过求解上述雷诺方程得到油膜面积划分区域中的油膜压力p。图6(a)所示为自适应网格下自由度为2 655时求得的油膜压力曲面。图6(b)所示为等距差分格式求解的自由度为2 025时的油膜压力曲面。因篇幅限制，通过其他自由度求得的油膜压力曲面文中并未给出，但表2中给出了2种数值解法在不同自由度下在求解域上的油膜压力和油膜压力在求解域上的压力积分值。

针对内燃机活塞缸套模型将2种算法下不同自由度数下的油膜压力积分做出对比，如图7所示绘制了2种算法下油膜压力积分曲线。

图6 自适应变步长差分法和传统有限差分法计算的油膜压力

表2自适应变步长差分法和传统有限差分法计算的油膜压力与压力积分值

Table 2 Oil film pressure and pressure integral valuescalculated by adaptive variable step differencealgorithm and traditional finite difference method

自由度数自适应变步长差分法油膜压力p/MPa油膜压力积分w/kN自由度数传统有限差分法油膜压力p/MPa油膜压力积分w/kN5721.1120.380 19001.110 30.380 41 4281.114 90.381 51 6001.113 10.381 22 6681.115 60.381 93 6001.114 80.381 85 2481.116 20.382 16 4001.115 30.3828 6921.116 40.382 18 1001.115 40.382 1