考慮空穴效應的廣義特征分解法解雷諾方程

耿煜,郭振,陳渭

(西安交通大學現代設計及轉子軸承系統教育部重點實驗室,710049,西安)

鉸接副作為必要的連接部件在機械系統中廣泛應用,由于受到加工誤差、材料變形、磨損等因素的影響,鉸接副大多存在間隙,導致機構連接處發生碰撞、振動,加快鉸接副的磨損,因此間隙鉸接副通常依靠潤滑將不利影響降到最低。

滑動軸承非線性油膜力模型是含潤滑間隙鉸接副多體系統動力學特性分析與研究的基礎和關鍵,油膜力的計算精度和速度將直接影響非線性動力學特性分析的效率。目前關于油膜力模型的研究已取得很大進展,涌現出了多種求解非線性油膜力的理論與方法,其中應用較為普遍的為有限差分法(FDM)[1-4]和有限元法(FEM)[5-8]。雖然數值方法精度高,但是計算效率很低,因此運用近似解析法求解油膜力的研究日漸增多,目前滑動軸承非線性油膜力近似解表達式主要有無限短軸承假設下的Capone模型[9]和Sommerfeld提出的無限長軸承油膜力模型[10],但是考慮到油膜的溫升和端泄,滑動軸承的寬徑比通常在0.8～1.0之間[11],現實中較少存在無限短或無限長軸承,實際應用受到限制,近年來,很多學者針對有限長軸承做了相當多的工作。Chasalevris等基于油膜假設,利用變分原理求出有限長軸承非線性油膜力的近似解析解,但是并沒有考慮空穴效應[12-14]。在流體潤滑問題中,由于發散間隙的存在會導致油膜出現負壓,當油膜達到最大可承受負壓時會發生破裂現象,進而產生空穴,所以不考慮空穴效應所獲得的油膜力是不可靠的[14]。張永芳等基于下游雷諾邊界條件,運用變分原理和分離變量法提出了一種針對有限長軸承的近似解析解方法,計算精度較高,但是計算速度相比傳統有限差分法并沒有得到很大的提升[11,15]。除此之外,還有諸如數據庫法和精確解析法等油膜力計算模型,但是由于通用性不強、計算效率較低等問題,并不適合于一般性徑向滑動軸承油膜力的求解。因此,建立考慮空穴效應、通用性較高、計算速度較快的滑動軸承非線性油膜力模型對于提高含潤滑間隙鉸接副多體系統動力學分析效率和滿足實際工程應用至關重要。

廣義特征分解法(PGD)主要是通過將求解域轉變為張量積函數之和來減少待求未知量個數,從而有效提高計算效率。將待求函數u分解為N維Q個張量積之和的形式,表達式如下

(1)

使用M個自由度對每個張量進行離散,待求未知量有NQM個,而如果使用傳統的FEM法或FDM法,待求未知量則為MN個。Amine等第一次使用PGD法對偏微分方程進行求解[16];Dumon等將PGD法成功應用到流體力學領域,并證明了PGD法的高效性[17];Tamellini等使用PGD法對穩態不可壓縮Navier-Stokes方程進行求解[18-19];Cherabi等使用PGD法求解穩態工況下的雷諾方程[20],但是并沒有考慮空穴對油膜力分布的影響,實際應用具有一定的局限性。

本文針對含潤滑間隙鉸接副中油膜力的計算問題,在文獻[20]的基礎上,建立了動載滑動軸承非線性油膜力的PGD計算模型,并引入雷諾邊界條件,充分考慮空穴效應對油膜力分布的影響。利用Matlab數值計算驗證了模型的高效性和穩定性,明晰了PGD法的適用范圍。將模型引入到含潤滑間隙鉸接副多體系統的計算中去,解決了模型的收斂問題,證明了模型的可靠性和工程應用價值。

1 計算模型的建立

1.1 動載雷諾方程

動載雷諾方程如下式

(2)

可以將雷諾方程轉換為以下的形式

(3)

即

(4)

將式(5)代入式(4),得

?x=R?θ;U=Rω

(5)

可得到式(4)圓柱坐標系下的表達式

(6)

油膜厚度表達式如下

h=c(1+εcosθ)

(7)

式中:θ為周向坐標;z為軸向坐標;h為油膜厚度;R為軸承半徑;μ為潤滑油黏度;ω為旋轉角速度;c為半徑間隙;ε為偏心率。

1.2 PGD求解雷諾方程

將雷諾方程式(6)轉換為等效積分的弱形式,若對于任意的權函數p′下式成立

(8)

則式(6)成立。

將油膜力函數p進行分解,獲得其近似解

(9)

式中:i為當前維數;Xi和Yi為第i維度下周向和軸向方向向量。一般情況下,維數越多,近似解越逼近精確解。

采用變方向迭代策略,式(9)轉換為

(10)

為了方便計算,將權函數p′設為

(11)

將式(10)(11)代入式(8),得

(12)

將式(12)進行整理,得

(13)

其中

將式(13)進一步化簡,得

(14)

式(10)、式(11)轉換為

(15)

(16)

將式(15)(16)代入式(8),得

(17)

同理,將式(17)進行整理,得

(18)

其中

將式(18)進一步化簡,得

(19)

1.3 油膜承載力的求解

采用中差分五點差分形式對式(14)(19)進行求解,將求解得到的Xi(θ)、Yi(z)代入式(9),算出油膜力分布函數p。

軸承承載力計算公式如下

(20)

其中

Fx=?pRsinθdθdz;Fy=?pRcosθdθdz

1.4 邊界條件

軸承的端泄邊、周期邊界和雷諾邊界的邊界條件分別如下

(21)

雷諾邊界條件假定在油膜完整區和油膜破裂區的交界線上壓力的法向梯度為0,且壓力等于空穴壓力pc,本文空穴壓力設為0。

(22)

1.5 收斂條件

ξ1和ξ2的計算表達式如下

(23)

(24)

1.6 計算程序

計算程序的基本流程如下:

步驟1 輸入3個軸承參數,潤滑油黏度、最大允許誤差λ1、最大允許維度數N=1;

步驟2 劃分網格節點,計算油膜厚度分布;

步驟8 由式(9)計算當前油膜力分布,若

步驟9 由式(21)計算誤差ξ1,如果ξ1<λ1,則轉到步驟10,否則q=q+1,重復步驟4～9。

步驟10 由式(20)計算軸承承載力,輸出油膜壓力分布。

2 計算模型的驗證

2.1 計算精度的驗證

Kumar和Booker[4]基于JFO邊界條件,采用有限元法對有限長軸承進行了油膜力求解,已有研究[21]表明,在偏心率小于0.02的情況下,雷諾邊界條件和JFO邊界條件所得到油膜力誤差在5%以內;當偏心率高于0.02時,兩種邊界條件所得結果幾乎一致。因此本文使用文獻[4]的計算參數和結果,來驗證基于雷諾邊界條件的PGD法的準確性,具體的軸承參數意義及數值見表1。

表1 文獻[4]中的軸承參數

將有限差分法、有限元法和PGD法計算結果與文獻[4]的計算結果進行對比,結果如表2所示,文中不作特殊說明FDM和PGD法均采用中差分五點差分格式進行數值求解。采用PGD法得到的油膜承載力和最高壓力的結果誤差均在10%以內。值得注意的是,隨著網格密度的加密,PGD法的計算精度將得到提高,但是相比FDM法計算時間得到了大幅縮減,這將在2.2節進行驗證。另外,將使用PGD、FDM、FEM法得出的油膜壓力分布圖進行比較,如圖1所示,可以看到圖像曲線具有較好的一致性,再一次驗證了基于雷諾邊界條件的PGD法的準確性。

表2 PGD、FDM、FEM法與文獻[4]方法的結果對比

圖1 軸向中心處油膜壓力分布圖

2.2 計算效率的驗證

相比FDM法,FEM法計算精度更高,處理復雜幾何結構時適應性強,但是FDM法使用簡單,在處理一般性滑動軸承方面應用更為廣泛,計算效率相比FEM法有本質上的提高,且其計算精度已經能夠滿足一般性工程需要,因此在這里將采用新模型的FGD法與FDM法進行比較,驗證FGD法的高效性。

設置相同的最大允許誤差λ1和SOR迭代因子分別為10-5和1.41,比較PGD和FDM法在不同網格密度下耗費的計算時間。由于PGD法受到初始隨機賦值的影響,為了嚴謹,每種網格密度下的PGD計算時間和承載力取3次平均值。如表3所示,在低網格密度下,PGD與FDM法的相對誤差較低,且隨著網格密度的增加,誤差控制在1%以內,計算效率最大可以提高3.59倍。在網格密度達到300×300以上時,兩種方法的計算誤差開始突然增大,通過承載力可以看出,FDM法算出的承載力發生明顯突變,而PGD法算出的承載力變化較小,這說明在相同的λ1條件下,PGD法具有更好的計算穩定性,隨著網格密度的提高,依舊保持了較高的計算精度,而FDM法由于網格密度過高,原有λ1已經無法滿足其計算精度的要求,導致結果過早收斂進而使得計算結果發生突變。

為了驗證上述結論,在600×600的固定網格密度情況下,采用不同的λ1進行了實驗,實驗結果如表4所示。在降低最大允許誤差后,PGD與FDM法的相對誤差始終控制在1%以內,最大加速比達到了4.42,表明PGD法具有相當高的計算精度和效率。結合表2和表3可以看出,與FDM法相比,PGD法不僅保持了較高的計算精度,而且具有更好的計算穩定性與更快的計算速度,在高網格密度和低允許誤差條件下,計算優勢將會更加明顯。

表3 不同網格密度下PGD和FDM法的計算結果對比

表4 不同允許誤差下PGD與FDM法的計算時間對比

3 PGD法在含潤滑間隙鉸接副多體系統中的應用

為了提高含潤滑間隙鉸接副多體系統的計算效率,符合工業應用所需要的實用性特點,引入模型對含潤滑鉸接副非線性油膜力進行計算。本文以曲柄滑塊為例,曲柄滑塊機構見圖2,相關參數見表5。

圖2 曲柄滑塊機構示意圖

表5 曲柄滑塊機構參數

采用歐拉-拉格朗日法對多體系統進行建模,引入Flores過渡模型來解決潤滑間隙鉸接副的工作狀態過渡問題,具體參數詳見表6,其中E與μ表示材料的彈性模量與泊松比。接觸力模型與摩擦力模型需要設定恢復系數與摩擦系數,本文分別取為0.9和0.2。Flores模型過渡臨界值設為2 μm,潤滑油黏度為0.08 Pa·s。具體的公式與計算步驟在此不作贅述,詳見文獻[1]。

表6 間隙鉸接副參數

在多體系統計算中采用FDM法求解雷諾方程時,往往使用上一時刻計算得出的壓力分布解作為下一時刻的初始值進行計算加速,使用PGD法時同樣使用上一次計算出的周向和軸向向量Xi、Yi來作為初始值進行計算加速,但是Xi、Yi不僅參與迭代,還會影響到迭代系數Aθ～Fθ,Az～Hz的計算,因此連續使用上一次計算解作初始解會產生累計誤差,導致算法難以收斂,因此需要在一定數量的時間步長后,對初始向量進行新一輪的隨機賦值,以消除累計誤差對計算所帶來的影響。

曲柄轉動3周,設置相同的λ1為10-5,SOR迭代因子為1.41,驗證PGD法在含潤滑間隙鉸接副多體系統計算中的可靠性與實用性。圖3為PGD與FDM法在30×30網格密度下的軸心軌跡、最小油膜厚度、滑塊速度隨時間變化的對比。從圖3可以看出,PGD與FEM法的計算結果相差不大,滑塊速度曲線近乎重合,可見PGD法具有較高的計算精度。

(a)軸心軌跡圖

(b)最小油膜厚度隨時間變化

(c)滑塊速度隨時間變化圖3 30×30網格下PGD與FDM法的計算結果對比

表7為兩種方法在不同網格密度下運轉3個周期的計算時間對比,可以看出PGD法對于計算的加速效果是非常明顯的,在30×30和40×40的網格密度下,PGD法平均一個周期的計算時間分別縮短了6 min和20 min,在實際計算中,為了觀察間隙鉸接副的潤滑效果,至少需要進行5個周期的運算,隨著機構復雜性的提高,達到穩定潤滑狀態所需要的周期數會進一步增加,因此采用新模型所節省的時間將是巨大的。

表7 在含潤滑間隙鉸接副多體系統中的計算時間對比

4 結 論

本文基于廣義特征分解法建立了動載有限長滑動軸承非線性油膜力模型,引入雷諾邊界條件充分考慮空穴對油膜壓力分布的影響,與現有文獻數據進行對比,證明了新模型的準確性與高效性,并研究了采用新模型方法的計算精度與計算效率隨網格密度的變化趨勢。以曲柄滑塊機構為例,采用新模型的方法對含潤滑間隙鉸接副進行動力學分析,證明了新模型的實際應用價值。可以得出以下結論:

(1)PGD法計算編程較為復雜,相較傳統計算方法更難理解,但是在較低網格密度下,PGD法與FDM法的承載力計算相對誤差最高不超過7%,在較高網格密度下,相對誤差控制在1%以內,具有較高的計算精度且穩定性相比FDM法更高。

(2)PGD法計算效率較高,在不損失精度的前提下,在不同網格密度下優勢明顯,計算速度相比FEM法提高2～5倍。

(3)使用PGD法對潤滑間隙鉸接副動力學分析,在30×30和40×40網格密度下,單個運算周期節省時間分別高達6 min和20 min,隨著網格密度和機械系統復雜性的提高,新模型的計算優勢將會更加顯著,證明了模型的可靠性與實用價值。

本文只是將PGD法應用于普通雷諾方程的求解,考慮表面粗糙度及溫度的影響,例如將PGD法引入平均流模型和熱彈流計算中是PGD法未來的發展方向。本文考慮空穴效應的PGD法只是針對普通滑動軸承而言,而對于織構化軸承而言,由于存在大量的收斂、發散間隙,油膜會反復破裂并再形成,這對算法的收斂問題會產生巨大挑戰,這是PGD法目前所面臨的難題。