基于層次模型的雷諾方程自適應(yīng)變步長差分算法研究*

(1.武漢科技大學(xué)冶金裝備及其控制教育部重點(diǎn)實(shí)驗(yàn)室 湖北武漢 430081；2.武漢科技大學(xué)機(jī)械傳動與制造工程湖北省重點(diǎn)實(shí)驗(yàn)室 湖北武漢 430081)

雷諾方程實(shí)際上是由不可壓縮黏性流體Navier-Stokes動量方程和連續(xù)性方程變形而得到的一個關(guān)于壓力的二階偏微分方程。在諸多工程問題中，流體潤滑問題的求解首先需要對雷諾方程進(jìn)行求解[1-2]。當(dāng)今，求解雷諾方程使用最為普遍的是有限差分法[3]，其通常用等分網(wǎng)格對整個模型區(qū)域進(jìn)行離散，并采用差商代替導(dǎo)數(shù)對方程進(jìn)行離散求解，計(jì)算機(jī)比較容易實(shí)現(xiàn)。對雷諾方程差分求解而言，網(wǎng)格的質(zhì)量對計(jì)算的結(jié)果有重要影響，王羽[4]在網(wǎng)格單元與節(jié)點(diǎn)相鄰的半步長上插入若干節(jié)點(diǎn)，在提高計(jì)算精度的同時也增大了計(jì)算量。針對無織構(gòu)軸承，于如飛和陳渭[5]用有限差分法對雷諾方程進(jìn)行了網(wǎng)格無關(guān)性分析，精度會有所提高，但與低精度相差甚微，且計(jì)算機(jī)耗時較高。采用傳統(tǒng)有限差分法求解雷諾方程，隨著步長的減小會大大增加計(jì)算量，對計(jì)算機(jī)要求較高并且浪費(fèi)資源。

針對傳統(tǒng)有限差分法求解雷諾方程的限制，BERGERA和OLIGERB提出了網(wǎng)格自適應(yīng)的思想[6]，即將網(wǎng)格點(diǎn)進(jìn)行最優(yōu)分布，對物理解變化較大的地方用小步長的網(wǎng)格，對需要重點(diǎn)模擬的區(qū)域也使用小步長網(wǎng)格，從而使網(wǎng)格步長根據(jù)待解物理量的特性和對應(yīng)模型的形狀而改變。JASTRAM和TESSMER[7]將該思想運(yùn)用到彈性波數(shù)值模擬中，提出了縱向網(wǎng)格步長逐漸變化的算法。劉志強(qiáng)等[8]將基于變分原理的自適應(yīng)網(wǎng)格生成技術(shù)，引入到了起伏海底速度模型的網(wǎng)格剖分中。宓鐵良等[9]采用自適應(yīng)網(wǎng)格，利用小波函數(shù)的多尺度分解的特性模擬了地震波的傳播，該方法通過判定在初始網(wǎng)格點(diǎn)處的小波系數(shù)與閾值之間的大小決定是否對網(wǎng)格細(xì)化。

本文作者基于早期學(xué)者的研究成果，提出了一種基于網(wǎng)格層次化思想的雷諾方程差分算子對網(wǎng)格進(jìn)行優(yōu)化的方法，并在此基礎(chǔ)開發(fā)了以壓力梯度變化與區(qū)域積分值為基準(zhǔn)的自適應(yīng)細(xì)化區(qū)域算法；以活塞-缸套潤滑系統(tǒng)和徑向滑動軸承為研究對象，詳細(xì)研究對比了不同算法下不同自由度對求解結(jié)果的影響，驗(yàn)證了基于層次化思想的自適應(yīng)變步長差分算法的優(yōu)越性及適應(yīng)性。

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

1.1 雷諾方程的形式及邊界條件

雷諾方程的一般形式是三維的，為了方便求解，根據(jù)不同的假設(shè)條件可以將雷諾方程簡化為二維形式和一維形式，在相關(guān)理論下都可以得到較精確的雷諾方程的解析解，但一維形式的假設(shè)往往不符合實(shí)際的工程設(shè)計(jì)，其求解結(jié)果必然與實(shí)際的工況有較大的出入。文中選取工程中最常見的雷諾方程的二維簡化形式進(jìn)行數(shù)值計(jì)算。以徑向滑動軸承為例，由Navier-Stokes方程和質(zhì)量連續(xù)性方程導(dǎo)出三維雷諾方程[10-12]的一般形式，如下式：

(1)

上式右邊三項(xiàng)依次稱為楔形項(xiàng)、伸張項(xiàng)及擠壓模項(xiàng)。對于一般徑向軸承往往是軸承固定，軸頸轉(zhuǎn)動，且外載荷為靜載荷。即u1=0,u2=u。另外，v2-v1可以用油膜厚度及其導(dǎo)數(shù)的形式來表示：

(2)

根據(jù)上述分析可以得到徑向滑動軸承雷諾方程的二維形式：

(3)

將求解域D歸一化到區(qū)間[0,1]×[0,1]后得到

(4)

其他形式的雷諾方程也可以用上述形式進(jìn)行表示，如缸套活塞-缸套模型對三維雷諾方程進(jìn)行降維、變量替換和歸一化得到雷諾方程后相對應(yīng)的系數(shù)項(xiàng)和右邊項(xiàng)：

求解域D的邊界為Γ=ΓD∪ΓN∪ΓR， 其中ΓD為Dirichlet邊界條件，ΓN為Neumann邊界條件,ΓR為Robin邊界條件。

1.2 變步長差分格式及迭代矩陣的建立

將所求的油膜潤滑區(qū)域網(wǎng)格化，形成若干單元區(qū)域，如圖1所示，用各個節(jié)點(diǎn)上的壓力值p構(gòu)成一階、二階差商來逼近導(dǎo)數(shù)，將二維偏微分方程和其邊界條件轉(zhuǎn)化成節(jié)點(diǎn)的差分方程。非均勻網(wǎng)格劃分下的油膜區(qū)域在X方向用j來表示列數(shù)，Y方向用i來表示行數(shù)。圖中總共有m×n個節(jié)點(diǎn)，其中在X方向有n個節(jié)點(diǎn)，在Y方向上有m個節(jié)點(diǎn)，并且任意節(jié)點(diǎn)的一階、二階導(dǎo)數(shù)可以由其周圍節(jié)點(diǎn)的壓力值進(jìn)行近似表示：

得到最終的差分格式為

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

(5)

式中各系數(shù)值隨節(jié)點(diǎn)位置的改變而改變，其中：

圖1 油膜網(wǎng)格劃分

非均勻網(wǎng)格劃分的求解域總共有m×n個節(jié)點(diǎn)，編碼順序?yàn)閺淖笸摇⒆韵露希来螢?,2,3,......,m×n。方程(5)是有限差分法的計(jì)算方程，在離散的網(wǎng)格化區(qū)域內(nèi)，每一個節(jié)點(diǎn)都對應(yīng)一個差分方程，根據(jù)邊界條件確定差分格式的各項(xiàng)系數(shù)并通過編碼的形式建立系數(shù)迭代矩陣。這樣便于調(diào)試程序時迭代矩陣系數(shù)項(xiàng)的修改，減小了編程的工作量并且保證了程序的穩(wěn)定性，適應(yīng)于不同模型的系數(shù)迭代矩陣的生成。其中第一類邊界的加載是索引處差分格式系數(shù)所在的編碼位置1，其余置0。第二類邊界條件加載將相應(yīng)的差分格式的系數(shù)替換成原來的兩倍。

2 自適應(yīng)變步長差分算法

2.1 自適應(yīng)變步長網(wǎng)格思想

網(wǎng)格劃分的質(zhì)量對雷諾方程數(shù)值求解精度和效率起著至關(guān)重要的作用，人們在網(wǎng)格的構(gòu)造過程中往往傾向于選擇大小均勻，長寬適當(dāng)、相互正交的網(wǎng)格；為了獲得高精度的計(jì)算結(jié)果往往需要大規(guī)模的網(wǎng)格來離散整個模型，這種方法計(jì)算靈活性較差，并造成了計(jì)算資源的浪費(fèi)。為了改善這一弊端，將自適應(yīng)變步長網(wǎng)格思想引入到雷諾方程的差分計(jì)算模型中。自適應(yīng)變步長網(wǎng)格方法將網(wǎng)格點(diǎn)進(jìn)行了最優(yōu)分布，做到了網(wǎng)格點(diǎn)的分布與物理解的耦合，從而最終達(dá)到在物理解變動較大的區(qū)域網(wǎng)格自動密集，在物理解變化平緩區(qū)域網(wǎng)格相對稀疏的目的，摒棄了傳統(tǒng)有限差分法所使用的均分網(wǎng)格的做法，在保證了精度的同時又提高了計(jì)算速率。

2.2 自適應(yīng)變步長網(wǎng)格構(gòu)造過程

自適應(yīng)變步長網(wǎng)格[13]構(gòu)造過程中最主要的問題便是根據(jù)網(wǎng)格品質(zhì)或者數(shù)值計(jì)算的特點(diǎn)來確定自適應(yīng)細(xì)分區(qū)域。以壓力梯度變化與單元區(qū)域積分值為基準(zhǔn)的細(xì)化區(qū)域算法來確定自適應(yīng)細(xì)分區(qū)域，并通過歸一化后灰度圖來對細(xì)分區(qū)域進(jìn)行直觀的表達(dá)。自適應(yīng)細(xì)分區(qū)域采用矩形網(wǎng)格進(jìn)行局部的細(xì)化，這樣要比采用非構(gòu)造網(wǎng)格計(jì)算速度快。

采用矩形網(wǎng)格對求解域進(jìn)行均勻劃分，步長為h0，這種矩形網(wǎng)格定義為基網(wǎng)格。所有細(xì)分的網(wǎng)格區(qū)域就是由基網(wǎng)格確定。基網(wǎng)格用Ω0表示。Ω0的下標(biāo)0表示尺度。根據(jù)劃分的網(wǎng)格的粗細(xì)來定義網(wǎng)格的尺度。

使用有限差分法進(jìn)行數(shù)值計(jì)算得到基網(wǎng)格上的數(shù)值解。首先以一階密度梯度為判定準(zhǔn)則來進(jìn)行網(wǎng)格細(xì)化區(qū)域的確定。如圖1計(jì)算i,j點(diǎn)的一階密度梯度，過程如下：

(6)

將求解域中的網(wǎng)格節(jié)點(diǎn)的一階密度梯度進(jìn)行歸一化處理。以圖1中節(jié)點(diǎn)pij為中心，將網(wǎng)格節(jié)點(diǎn)的一階梯度均分給周圍四塊網(wǎng)格區(qū)域，計(jì)算出每一個網(wǎng)格點(diǎn)僅僅在梯度方面對周圍的四塊單元網(wǎng)格區(qū)域貢獻(xiàn)值。由于每一個網(wǎng)格區(qū)域含有4個網(wǎng)格點(diǎn)，則這4個網(wǎng)格點(diǎn)的梯度密度都會影響到自適應(yīng)細(xì)分區(qū)域的確定。將每個網(wǎng)格區(qū)域的4個節(jié)點(diǎn)的貢獻(xiàn)值疊加可以得到單元網(wǎng)格區(qū)域的一階密度梯度。將網(wǎng)格節(jié)點(diǎn)梯度的變化轉(zhuǎn)移到網(wǎng)格區(qū)域上以此來更好地確定細(xì)分區(qū)域。然后將滿足下式中的網(wǎng)格區(qū)域過濾出來：

W≥ηWmax

(7)

其中η根據(jù)具體的求解和需要的效果確定。

按照上述步驟過濾出來的網(wǎng)格細(xì)分區(qū)域分布可能比較集中，也可能比較分散，細(xì)分區(qū)域分別用Ω1來表示，則Ω1可以表示為

Ω1=∪Ω1，j

(8)

j表示需要進(jìn)行網(wǎng)格細(xì)分的區(qū)域。例如在判定準(zhǔn)則1下的網(wǎng)格灰度圖如圖2所示。

圖2 一階梯度密度灰度圖

然后通過單元網(wǎng)格的油膜壓力積分來二次確定細(xì)化區(qū)域。通常采用高斯積分法[14]對0尺度各個網(wǎng)格的油膜進(jìn)行積分計(jì)算，得到細(xì)分區(qū)域的油膜壓力積分S0。

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

(9)

其中p(x,y)為積分點(diǎn)的函數(shù)值。

求解時候在每個網(wǎng)格單元上選擇一些點(diǎn)作為積分點(diǎn)，求出其壓力函數(shù)在網(wǎng)格節(jié)點(diǎn)上的數(shù)值，最后對這些點(diǎn)進(jìn)行加權(quán)求和，求出積分點(diǎn)的油膜壓力函數(shù)值，二維高斯積分的表達(dá)式如下：

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

(10)

其中積分點(diǎn)xi、yi和權(quán)函數(shù)Hi、Hj按相應(yīng)的算法選取。

這樣將網(wǎng)格細(xì)分的區(qū)域弱化成點(diǎn)，對單元網(wǎng)格的油膜壓力積分再進(jìn)行歸一化對比后排序，就可以清晰明了地看到每塊網(wǎng)格的壓力積分的大小，從而將滿足式(11)中的網(wǎng)格區(qū)域過濾出來。

S1≥γSmax

(11)

其中γ根據(jù)具體的求解和需要的效果確定。

按照油膜壓力積分過濾出來的網(wǎng)格細(xì)分區(qū)域用Ω2來表示，Ω2來表示需要進(jìn)行網(wǎng)格細(xì)分的區(qū)域。則Ω2可以表示為

Ω2=∪Ω2，j

(12)

j表示在判定準(zhǔn)則2下確定網(wǎng)格細(xì)分的各個區(qū)域。例如在判定準(zhǔn)則2下的灰度圖如圖3所示。

圖3 單元網(wǎng)格壓力積分灰度圖

判定準(zhǔn)則1是通過將點(diǎn)的一階梯度密度映射到網(wǎng)格區(qū)域范圍內(nèi)，而判定準(zhǔn)則2是將油膜區(qū)域弱化成點(diǎn)比較油膜壓力積分的大小，各自劃分的Ω1和Ω2區(qū)域可能集中也可能比較分散。

假設(shè)1：當(dāng)區(qū)域都比較集中時，則最后細(xì)分的區(qū)域?yàn)?/p>

Ω=Ω1∪Ω2

(13)

假設(shè)2：當(dāng)Ω1區(qū)域較集中時，則以Ω1區(qū)域?yàn)橹饕獏⒖紖^(qū)域，根據(jù)Ω2區(qū)域的分布進(jìn)行適當(dāng)?shù)臄U(kuò)展,由此得到最終的細(xì)分區(qū)域Ω。

假設(shè)3：當(dāng)區(qū)域都比較分散時，則將細(xì)分區(qū)域進(jìn)行分開處理，但數(shù)值實(shí)驗(yàn)證明這種情況很少。

進(jìn)行變步長計(jì)算時，將區(qū)域網(wǎng)格按照假設(shè)區(qū)域的灰度值大小(假設(shè)1就按照各自的區(qū)域灰度值排序，假設(shè)2就按照Ω1區(qū)域灰度值進(jìn)行排序)分為3塊區(qū)域，步長依次為h1,h2,h3，h1=2h2=3h3=1/3h0，隨著細(xì)分的進(jìn)行將步長縮小為原來的1/2。

3 自適應(yīng)變步長差分算法求解流程

在MATLAB 中開發(fā)出求解程序，其求解流程圖如圖4所示，步驟如下：

(1)輸入方程的常值參數(shù)，即參數(shù)設(shè)置，包括活塞裙部半徑、缸套半徑、活塞裙部長度、活塞裙部上下端半徑以及潤滑油黏度等工況參數(shù)。

(2)初始化網(wǎng)格，按照圖1所示的網(wǎng)格示意圖進(jìn)行編碼，根據(jù)公式(5)求出系數(shù)迭代矩陣編碼位置對應(yīng)的系數(shù)，從而生成系數(shù)迭代矩陣和右邊項(xiàng)，計(jì)算出初始網(wǎng)格的油膜壓力。

(3)計(jì)算單元網(wǎng)格的油膜壓力一階密度梯度和油膜壓力積分，根據(jù)自適應(yīng)網(wǎng)格構(gòu)造過程中的準(zhǔn)則1和準(zhǔn)則2來確定是否存在細(xì)分區(qū)域。若存在細(xì)分區(qū)域，繼續(xù)細(xì)分生成第2層網(wǎng)格，然后按照準(zhǔn)則1和準(zhǔn)則2進(jìn)一步確定細(xì)分區(qū)域，細(xì)化網(wǎng)格，最終得出自適應(yīng)變步長網(wǎng)格下的系數(shù)矩陣和右邊項(xiàng)。

(4)進(jìn)行迭代計(jì)算得到油膜壓力p的數(shù)值，并與初始化網(wǎng)格的油膜壓力進(jìn)行對比，若符合條件就輸出，不符合條件就調(diào)整初始網(wǎng)格。迭代可以采用SOR迭代，要注意松弛因子的選取。

(5)設(shè)定如下的收斂判定準(zhǔn)則：

(14)

式中：σ即為允許的相對誤差值，取為1×10-6。

4 數(shù)值算例

4.1 內(nèi)燃機(jī)活塞-缸套潤滑模型

活塞缸套之間的膜潤滑主要涉及通過求解雷諾方程得到油膜壓力。圖5(a)所示為活塞-缸套系統(tǒng)，其油膜潤滑區(qū)域如圖5(b)所示[15]。

因?yàn)樽笥覞櫥瑓^(qū)域相對于平面φ=0都是對稱的，所以只需考慮前面部分。通常，當(dāng)活塞周期性地移動時，由于負(fù)載變化，密度和黏度的大小會發(fā)生變化。但為了簡化數(shù)值模型，這里將密度和黏度的值作為常數(shù)。求解得出油膜壓力的計(jì)算參數(shù)如表1所示。

圖5 活塞-缸套潤滑系統(tǒng)示意圖

參數(shù)說明值單位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對時間的導(dǎo)數(shù)-0.243 5×10-2m·s-1debeb對時間的導(dǎo)數(shù)-0.844 61×10-2m·s-1v活塞往復(fù)運(yùn)動速度-10.654 3m·s-1Cb活塞銷到活塞冠部垂直距離6×10-3mCc活塞銷到活塞中心線距離0m

假設(shè)v0=0,vh=v(其正方向?yàn)樯闲?，潤滑劑無側(cè)漏，因此u0=uh=0。那么有：

(15)

油膜厚度h(φ,y)采用橢圓法進(jìn)行計(jì)算：

h(φ,y)=R-

將域D歸一化為[0,1]×[0,1]，進(jìn)行如下轉(zhuǎn)換：右側(cè)(因?yàn)閷ΨQ，只考慮前面部分)

邊界條件為

通過求解上述雷諾方程得到油膜面積劃分區(qū)域中的油膜壓力p。圖6(a)所示為自適應(yīng)網(wǎng)格下自由度為2 655時求得的油膜壓力曲面。圖6(b)所示為等距差分格式求解的自由度為2 025時的油膜壓力曲面。因篇幅限制，通過其他自由度求得的油膜壓力曲面文中并未給出，但表2中給出了2種數(shù)值解法在不同自由度下在求解域上的油膜壓力和油膜壓力在求解域上的壓力積分值。

針對內(nèi)燃機(jī)活塞缸套模型將2種算法下不同自由度數(shù)下的油膜壓力積分做出對比，如圖7所示繪制了2種算法下油膜壓力積分曲線。

圖6 自適應(yīng)變步長差分法和傳統(tǒng)有限差分法計(jì)算的油膜壓力

表2自適應(yīng)變步長差分法和傳統(tǒng)有限差分法計(jì)算的油膜壓力與壓力積分值

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

自由度數(shù)自適應(yīng)變步長差分法油膜壓力p/MPa油膜壓力積分w/kN自由度數(shù)傳統(tǒng)有限差分法油膜壓力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

圖7 自適應(yīng)變步長差分算法與傳統(tǒng)有限差分法計(jì)算的油膜壓力積分值

可以看出，自適應(yīng)變步長差分算法下的收斂速度明顯較快，在較小的自由度下就能達(dá)到收斂后的解，而等分網(wǎng)格劃分下傳統(tǒng)有限差分法下求解達(dá)到收斂時需要更多的自由度，大大地降低了求解效率。

4.2 徑向滑動軸承模型

工作在穩(wěn)態(tài)條件下的徑向滑動軸承的潤滑模型如圖8所示。

圖8 徑向滑動軸承潤滑模型

徑向滑動軸承油膜壓力的計(jì)算參數(shù)如表3所示。

表3 計(jì)算參數(shù)

層流狀態(tài)下的二維不可壓縮流體動壓潤滑雷諾方程[16-17]為

(16)

將求解域D歸一化，得到雷諾方程的系數(shù)項(xiàng)和右邊項(xiàng)

邊界條件為

p(X,Y)|X=0=0；p(X,Y)|X=1=0;p(X,Y)|Y=0=0,p(X,Y)|Y=1=0

光滑滑動軸承的油膜厚度公式為

h=C+e·cosφ

(17)

針對徑向滑動軸承模型其求解結(jié)果如圖9所示，其中圖9(a)表示運(yùn)用自適應(yīng)變步長差分算法生成的自由度為3 876時的油膜壓力曲面，圖9(b)表示等分網(wǎng)格下傳統(tǒng)有限差分法求解的自由度為2 500時的油膜壓力曲面。

圖9 自適應(yīng)變步長差分法和傳統(tǒng)有限差分法計(jì)算的油膜壓力

表4中給出了兩種數(shù)值解法在不同自由度下求得的油膜壓力以及油膜壓力在求解域上的壓力積分值。

表4 自適應(yīng)變步長差分法和傳統(tǒng)有限差分法計(jì)算的油膜壓力與壓力積分值

同活塞缸套潤滑系統(tǒng)一樣，繪制2種計(jì)算方法下的油膜壓力積分圖，如圖10所示。可見，自適應(yīng)變步長差分算法能夠在低自由度下快速收斂，且能夠精確地表達(dá)油膜壓力，對不同雷諾方程模型具有很好的適應(yīng)性。

圖10 自適應(yīng)變步長差分算法與傳統(tǒng)有限差分法計(jì)算的油膜壓力積分值

5 結(jié)論

(1)自適應(yīng)變步長差分算法實(shí)現(xiàn)了用較小的自由度表達(dá)油膜壓力曲面，在保證精度的同時有效地降低了網(wǎng)格離散所需自由度數(shù)，收斂速度和計(jì)算效率較傳統(tǒng)均分網(wǎng)格得到了顯著的提升。并且在不同的計(jì)算模型中都得到了驗(yàn)證，具備較高的靈活性與適配性。

(2)但自適應(yīng)變步長網(wǎng)格算法采用SOR進(jìn)行單層網(wǎng)格迭代，并沒有充分有效地運(yùn)用細(xì)分過程中的所有細(xì)分網(wǎng)格，后期考慮將多重網(wǎng)格算法與自適應(yīng)變步長差分算法相結(jié)合，在自適應(yīng)差分算法松弛因子對多重網(wǎng)格算法的影響展開深入研究。