質量守恒潤滑模型的快速數值計算方法*

顧春興 戴 黎

(上海理工大學機械工程學院 上海 200093)

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

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

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

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

1 數學模型

1.1 控制方程

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

(1)

(2)

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

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

1.2 加速收斂的數值計算方法

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

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

(3)

(4)

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

如圖1所示，沿x和y方向對研究對象進行網格劃分。網格點沿x方向的位置用i進行標記，沿y方向的位置用j進行標記。網格點(i,j)處的油膜壓力可以用p(i,j)表示，網格點(i,j)處的油膜厚度可以用h(i,j)表示。流體潤滑模型數值方程可采用有限體積法進行離散。如圖1所示，分別使用下標P、E、W、N、S、e、w、n和s來替換下標(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 流體動壓潤滑模擬中的網格點分布情況

基于有限體積法，公式(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+μP),μs=0.5(μS+μP)。另外，帶上標pre的項代表前一時間步的項,Δt則表示時間步長。

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

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

(6)

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

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

(7)

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

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

(8)

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

(9)

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

圖2 網格細化策略示意

在GR策略中，較低網格層(亦稱粗網格)上的計算結果應該首先進行插值，然后才能用作較高網絡層(細網格)上的計算初始值。基于流體潤滑壓力和流體飽和度的特點，文中所采用插值算法為非線性插值運算，相應的表達式為

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)

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

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

(13)

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

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

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

2 結果和討論

2.1 油膜厚度的表征

油膜厚度的表征是流體潤滑數值計算的關鍵。摩擦學系統中油膜厚度是最小油膜厚度和輪廓形狀之和。其表達式如下：

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

(14)

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

織構化的表面輪廓的表達式[16]如下:

(15)

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

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

圖4 織構輪廓的數值重建云圖

2.2 數值模擬條件

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

表1 數值模擬條件

2.3 網格收斂性研究

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

圖5 采用不同網格時的最大流體壓力的計算結果

2.4 標準方法和改進方法結果比較

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

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

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

圖7示出了3種算例采用不同方法進行求解所需的計算時間。計算時間是通過CPU時間衡量的，可以發現，與標準方法相比，考慮多層網格的快速數值計算方法可以顯著減少計算時間。當采用快速數值計算方法時，每個非目標網格層的計算都將有助于目標網格層上的精確初始值的獲取。研究結果表明，對于不同的算例，與傳統方法相比，快速數值計算方法可縮短73%～81%的數值計算時間。在快速數值計算方法中，由于初始值更加接近精確值，目標網格上的計算能夠迅速收斂。事實上，在一個完整的計算流程中，流體壓力分布的更新需要進行多次，即需要多次求解流體潤滑模型。因此，加快求解流體動壓潤滑模型，就能顯著減少整個計算流程的時間?？梢姡捎每焖贁抵涤嬎惴椒▉砬蠼饬黧w潤滑模型可明顯提高求解效率。

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

3 結論

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

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