顧春興 戴 黎

(上海理工大學(xué)機(jī)械工程學(xué)院 上海 200093)

在大多數(shù)機(jī)械系統(tǒng)中，流體動壓潤滑是廣泛存在的。這些機(jī)械系統(tǒng)包括徑向滑動軸承系統(tǒng)、止推軸承系統(tǒng)、活塞環(huán)/缸套系統(tǒng)和機(jī)械密封系統(tǒng)等。一般情況下，當(dāng)機(jī)械系統(tǒng)的2個接觸面之間存在收斂和發(fā)散間隙時，2個接觸面間將會產(chǎn)生流體動壓力。此時，潤滑油能有效地分離接觸表面以減少摩擦[1]。

近年來，表面織構(gòu)作為一種降低摩擦磨損的有效手段，在學(xué)術(shù)界和工業(yè)界得到了廣泛的應(yīng)用和研究[2-7]。然而，在涉及流體動壓潤滑的數(shù)值模擬中，流體潤滑區(qū)域可能會出現(xiàn)空穴現(xiàn)象[8-10]。同時，由于表面織構(gòu)的存在，潤滑油膜破裂后會發(fā)生潤滑油膜的重構(gòu)。因此，在計算流體壓力時，需要對空穴邊界進(jìn)行質(zhì)量守恒處理，以準(zhǔn)確地完成摩擦學(xué)系統(tǒng)潤滑性能的預(yù)測[11]。在過去的幾十年里，人們提出了不同的質(zhì)量守恒空化算法。其中，由JAKOBSSON和FLOBERG[12]、OLSSON[13]提出的JFO(Jacobsson-Floberg-Olsson)空化邊界條件具有特殊的地位，已經(jīng)在很多研究中得到了成功應(yīng)用。

然而，采用質(zhì)量守恒的空化邊界條件將會使流體動壓潤滑模型的求解變得不穩(wěn)定，計算時間增加。此外，潤滑模型的求解時間還與網(wǎng)格密度有關(guān)。在表面織構(gòu)的模擬中，常需要采用精密的網(wǎng)格才能精確地表征表面織構(gòu)。網(wǎng)格數(shù)的增加將導(dǎo)致數(shù)值計算的時間成倍增加。因此，在求解含有質(zhì)量守恒空化邊界條件的雷諾方程時，高效和穩(wěn)定的求解方法是迫切需要的。

為了解決由網(wǎng)格細(xì)化和采用質(zhì)量守恒空化邊界條件引起的潤滑模型不易收斂和計算時間過長的問題，本文作者通過結(jié)合Fischer-Burmister-Newton-Schur (FBNS)方法和網(wǎng)格細(xì)化(GR)策略，提出一種高效的數(shù)值計算方法。采用這種數(shù)值計算方法將顯著減少流體動壓潤滑模型的求解時間。

1 數(shù)學(xué)模型

1.1 控制方程

為了預(yù)測油膜壓力，文中將采用具有JFO空化邊界條件的二維雷諾方程[12-13]。方程式如下：

(1)

(2)

式中：p為油膜壓力；h為油膜厚度；?為流體飽和度；ρ為潤滑油密度；μ為潤滑油黏度；t為時間步長。

流體飽和度?和油膜壓力p之間的關(guān)系可由公式(2)揭示。如公式(2)所示，當(dāng)油膜壓力p大于空化壓力(pca=0)時，?設(shè)為1。當(dāng)油膜壓力p等于空化壓力時，流體飽和度?將小于1。

1.2 加速收斂的數(shù)值計算方法

在含油潤滑機(jī)械系統(tǒng)的數(shù)值模擬中，求解流體潤滑模型往往需要多次迭代。因此，流體潤滑模型的求解速度是非常關(guān)鍵的。數(shù)值模擬的運(yùn)行時間取決于流體潤滑模型的收斂速度。因此，迫切需要加快流體潤滑模型計算速度的方法。

文中所提出的高效的數(shù)值計算方法，采用如下2種方法，以提高流體潤滑模型的計算效率。一方面，計算中將采用FBNS方法[14]，以提高流體潤滑模型求解的穩(wěn)定性和收斂性。FBNS方法的采用，可以將流體潤滑模型(一個約束問題)的求解轉(zhuǎn)化為一個無約束問題的求解，而無約束問題的求解則可以直接采用現(xiàn)有的商用求解器。文獻(xiàn)[14-15]中已經(jīng)給出了FBNS方法的詳細(xì)步驟，文中不再詳述，下面僅簡要介紹FBNS方法的求解步驟。按照FBNS方法，首先需要對考慮JFO空化邊界的雷諾方程進(jìn)行修改[12-13]。質(zhì)量守恒的空穴邊界條件將采用一個互補(bǔ)約束來代替。相應(yīng)的表達(dá)式如下：

(3)

(4)

其中θ=1-?。當(dāng)且僅當(dāng)p、θ≥0和pθ=0時，公式(4)所涉及的約束條件才能滿足。根據(jù)文獻(xiàn)[14]可知，采用FBNS方法求解流體潤滑模型，其計算速度將明顯快于傳統(tǒng)方法(2個數(shù)量級)。

如圖1所示，沿x和y方向?qū)ρ芯繉ο筮M(jìn)行網(wǎng)格劃分。網(wǎng)格點(diǎn)沿x方向的位置用i進(jìn)行標(biāo)記，沿y方向的位置用j進(jìn)行標(biāo)記。網(wǎng)格點(diǎn)(i,j)處的油膜壓力可以用p(i,j)表示，網(wǎng)格點(diǎn)(i,j)處的油膜厚度可以用h(i,j)表示。流體潤滑模型數(shù)值方程可采用有限體積法進(jìn)行離散。如圖1所示，分別使用下標(biāo)P、E、W、N、S、e、w、n和s來替換下標(biāo)(i,j)、(i+1,j)、(i-1,j)、(i,j+1)、(i,j-1)、(i+0.5,j)、(i-0.5,j)、(i,j+0.5)和(i,j-0.5)。

圖1 流體動壓潤滑模擬中的網(wǎng)格點(diǎn)分布情況

基于有限體積法，公式(3)可以分解成[15]：

6ΔyU[(1-θW)ρWhW-(1-θP)ρPhP]+

(5)

其中he=0.5(hE+hP),hw=0.5(hW+hP),hn=0.5(hN+hP),hs=0.5(hS+hP),ρe=0.5(ρE+ρP),ρw=0.5(ρW+ρP),ρn=0.5(ρN+ρP),ρs=0.5(ρS+ρP),μe=0.5(μE+μP),μw=0.5(μW+μP),μn=0.5(μN(yùn)+μP),μs=0.5(μS+μP)。另外，帶上標(biāo)pre的項(xiàng)代表前一時間步的項(xiàng),Δt則表示時間步長。

因此，對公式(5)進(jìn)行重新排列，可得如下表達(dá)式：

ANPN+ASPS+APPP+AWPW+AEPE+BPθP+BWθW+CP=0

(6)

公式(6)揭示了油膜壓力和流體飽和度之間的線性關(guān)系，因此離散化結(jié)果可以表示為

G(p,θ)=Ap+Bθ+c=0

(7)

式中：A是包含AN、As、AP、AE和AW的信息的矩陣；矩陣B包含BW和BP的信息；c是包含邊界條件和CP信息的向量。

此外，公式(4)可離散為如下的表達(dá)式：

(8)

公式(8)所示的離散化結(jié)果可用如下表達(dá)式表示：

(9)

另一方面，如圖2所示，可采用GR(網(wǎng)格細(xì)化)策略以進(jìn)一步加快流體潤滑模型的計算速度。在GR策略中，粗網(wǎng)格上的p和θ的最終解將作為細(xì)網(wǎng)格上p和θ的初始值。通過這種方式，可保證細(xì)網(wǎng)格上p和θ的初始值接近于最終解。根據(jù)牛頓迭代法的特點(diǎn)，當(dāng)初始值足夠接近最終解時，牛頓迭代法具有二次收斂。另外，由于網(wǎng)格數(shù)目較少，粗網(wǎng)格上數(shù)值計算的收斂速度相對較快。當(dāng)合理選擇網(wǎng)格層數(shù)時，GR策略可以有效地縮短計算時間。

圖2 網(wǎng)格細(xì)化策略示意

在GR策略中，較低網(wǎng)格層(亦稱粗網(wǎng)格)上的計算結(jié)果應(yīng)該首先進(jìn)行插值，然后才能用作較高網(wǎng)絡(luò)層(細(xì)網(wǎng)格)上的計算初始值?；诹黧w潤滑壓力和流體飽和度的特點(diǎn)，文中所采用插值算法為非線性插值運(yùn)算，相應(yīng)的表達(dá)式為

uk(i,j)=uk-1(I,J),i=2I,j=2J

(10)

9uk-1(I+1,J)-uk-1(I+2,J)),i=2I+1,j=2J

(11)

9uk-1(I,J+1)-uk-1(I,J+2)),i=2I,j=2J+1

(12)

她的高燒退了一些，但仍然燙著，青辰繼續(xù)給她煎藥，喂藥，喂米湯或者面湯，然后再給她擦洗四肢，一套忙活下來，大半天的時間便已過去。

uk-1(I,J-1)+uk-1(I,J+2)+uk-1(I+1,J-1)+

(13)

其中uk(i,j)是第k網(wǎng)格層中的參數(shù)，而uk-1(I,J)位于第k-1網(wǎng)格層中。

圖3總結(jié)了流體潤滑模型的快速數(shù)值計算流程，其中Gr為所采用的網(wǎng)格層數(shù)。在每一網(wǎng)格層的計算中，p和θ值需要不斷更新，直到更新的值(在當(dāng)前迭代步)和前面的值(在前面的迭代步)之間的偏差小于或者等于一個允許的精度(10-6)。

圖3 流體潤滑模型的快速數(shù)值計算流程

2 結(jié)果和討論

2.1 油膜厚度的表征

油膜厚度的表征是流體潤滑數(shù)值計算的關(guān)鍵。摩擦學(xué)系統(tǒng)中油膜厚度是最小油膜厚度和輪廓形狀之和。其表達(dá)式如下：

h(x,y)=h0+hprof(x,y)

(14)

式中：hprof(x,y)是輪廓形狀。

織構(gòu)化的表面輪廓的表達(dá)式[16]如下:

(15)

式中：Λ=|x′cosα±y′sinα|cosα[16]；α代表凹槽角度；dg是平均凹槽深度；wg表示平均凹槽寬度；x′和y′都在局部坐標(biāo)軸上，而局部坐標(biāo)位于織構(gòu)單元的中心。

圖4給出了織構(gòu)化輪廓的數(shù)值重建云圖。其中，量綱一X可由X=(2x-b)/b，而量綱一Y可由Y=(2y-l)/l計算，b表示接觸面的寬度，l表示接觸面的長度。

圖4 織構(gòu)輪廓的數(shù)值重建云圖

2.2 數(shù)值模擬條件

為了評估所提出的快速數(shù)值計算方法的有效性，進(jìn)行了一系列穩(wěn)態(tài)仿真試驗(yàn)。在這些算例中，假定滑動速度為10 m/s，潤滑油黏度為0.004 Pa·s，潤滑油密度為860 kg/m3。同時，織構(gòu)輪廓尺寸如下：溝槽深度為3 μm，溝槽寬度為40 μm，凹槽角度設(shè)為45°?？棙?gòu)單元尺寸為200 μm×200 μm，接觸面長度為800 μm，接觸面寬度為400 μm。文中對3種算例進(jìn)行數(shù)值計算，以評估所提出的快速數(shù)值計算方法的求解效率。3種算例所采用的最小油膜厚度分別為0.141、0.283和0.566 μm。即，3種算例的最小油膜厚度值分別為綜合表面粗糙度的0.5倍、1倍和2倍。數(shù)值模擬條件見表1中。

表1 數(shù)值模擬條件

2.3 網(wǎng)格收斂性研究

在評估所提出的快速數(shù)值計算方法的計算效率之前，首先應(yīng)該進(jìn)行網(wǎng)格收斂性研究，以評價網(wǎng)格數(shù)量對數(shù)值計算結(jié)果的影響。在流體動壓潤滑模擬中，網(wǎng)格必須足夠精密，才能精確地表征表面輪廓。如圖5所示，當(dāng)網(wǎng)格大于1 200×600時，最大流體壓力的計算結(jié)果將趨于穩(wěn)定，逐漸接近于精確值。因此，在下面的數(shù)值計算中，將采用1 200×600網(wǎng)格。

圖5 采用不同網(wǎng)格時的最大流體壓力的計算結(jié)果

2.4 標(biāo)準(zhǔn)方法和改進(jìn)方法結(jié)果比較

針對3種算例，采用標(biāo)準(zhǔn)方法和文中所提出的快速數(shù)值計算方法進(jìn)行計算。在標(biāo)準(zhǔn)方法中，將采用FBNS方法來計算流體壓力。在快速數(shù)值計算方法中，將結(jié)合FBNS方法和GR策略來加速流體壓力和流體飽和度的計算。同時，還將分析網(wǎng)格層數(shù)對計算速度的影響。根據(jù)前文針對網(wǎng)格收斂問題的分析，以下模擬將采用1 200×600網(wǎng)格。數(shù)值計算是在一臺配備英特爾酷睿i7-7700 CPU的計算機(jī)(每個核心的時鐘頻率為3.60 GHz，內(nèi)存為24 GB)上進(jìn)行的。

圖6顯示了3種算例的流體壓力分布情況。需要指出，標(biāo)準(zhǔn)方法和快速數(shù)值計算方法所預(yù)測的流體壓力結(jié)果是完全一致的。這是因?yàn)椋喝鐖D3所示，標(biāo)準(zhǔn)方法和快速數(shù)值計算方法的主要差異在于是否采用GR策略。在快速數(shù)值計算方法中增加了GR策略，按照GR策略，為了使牛頓迭代法具有二次收斂，細(xì)網(wǎng)格(目標(biāo)網(wǎng)格)上的p和θ初始值是由粗網(wǎng)格上p和θ的最終解通過插值獲得的。然而，無論是標(biāo)準(zhǔn)方法還是改進(jìn)方法，p和θ最終解的計算總是在細(xì)網(wǎng)格上進(jìn)行的。因此，這2種方法所求的結(jié)果是完全一致的，其主要差異在于求解效率。

圖6 3種算例的流體壓力分布

圖7示出了3種算例采用不同方法進(jìn)行求解所需的計算時間。計算時間是通過CPU時間衡量的，可以發(fā)現(xiàn)，與標(biāo)準(zhǔn)方法相比，考慮多層網(wǎng)格的快速數(shù)值計算方法可以顯著減少計算時間。當(dāng)采用快速數(shù)值計算方法時，每個非目標(biāo)網(wǎng)格層的計算都將有助于目標(biāo)網(wǎng)格層上的精確初始值的獲取。研究結(jié)果表明，對于不同的算例，與傳統(tǒng)方法相比，快速數(shù)值計算方法可縮短73%～81%的數(shù)值計算時間。在快速數(shù)值計算方法中，由于初始值更加接近精確值，目標(biāo)網(wǎng)格上的計算能夠迅速收斂。事實(shí)上，在一個完整的計算流程中，流體壓力分布的更新需要進(jìn)行多次，即需要多次求解流體潤滑模型。因此，加快求解流體動壓潤滑模型，就能顯著減少整個計算流程的時間?？梢?，采用快速數(shù)值計算方法來求解流體潤滑模型可明顯提高求解效率。

圖7 3種算例采用不同方法所需的計算時間

3 結(jié)論

(1)提出一種高效精確的面向流體潤滑模型的快速數(shù)值計算方法，該方法結(jié)合了Fischer-Burmsister-Newton-Schur(FBNS)方法和網(wǎng)格細(xì)化(GR)策略，可以高效地求解流體潤滑模型。

(2)采用傳統(tǒng)數(shù)值方法和快速數(shù)值計算方法求解不同算例，結(jié)果表明，采用快速數(shù)值計算方法得到的流體壓力與傳統(tǒng)方法得到的結(jié)果完全一致；同時，采用快速數(shù)值計算方法能夠有效減少計算時間，計算時間下降可達(dá)73%～81%。在一個完整的計算流程中，會執(zhí)行多次潤滑模型的求解。因此，采用所提出的快速數(shù)值計算方法，將能大大減少整個計算流程的計算時間。