等幾何多重網(wǎng)格法在雷諾方程中的應用

李子強 ，羅會信 ，左兵權 ，張 希

（1.武漢科技大學機械自動化學院，冶金裝備及其控制教育部重點實驗室，湖北 武漢 430081；2.武漢科技大學機械自動化學院，機械傳動與制造工程湖北省重點實驗室，湖北 武漢 430081）

1 引言

對于雷諾方程的求解，有限差分法[3]和有限體積法[4]對于區(qū)域的連續(xù)性要求較為嚴格，結果也容易受到網(wǎng)格插值步長、迭代方法的影響，從而導致求解精度不高；有限元法[2]雖在精度上有提高，但在離散化過程中所花費的時間較長，效率低。等幾何法采用CAD 系統(tǒng)中的精確幾何模型來執(zhí)行物理場分析，避免了數(shù)值計算的二次建模，節(jié)省了求解域的離散過程所花的時間，在工程應用中有較大的優(yōu)勢。這些優(yōu)勢包括精確的幾何建模、高階連續(xù)性、簡單的網(wǎng)格劃分和網(wǎng)格細化等。

等幾何法可以在較少的自由度下得到同比于有限元法的求解精度。但隨著精度要求的進一步提高，高階的等幾何分析模型依然面臨大矩陣，低求解效率的困難。現(xiàn)今大部分的矩陣求解多采用迭代方法計算，多重網(wǎng)格法是普遍采用的方法之一。因此，為了加快高精度要求下等幾何法求解雷諾方程效率，將多重網(wǎng)格法與等幾何法相結合來對雷諾方程進行求解。

本研究選擇雷諾方程作為求解對象，將等幾何法，多重網(wǎng)格法與雷諾方程進行適配性研究：首先對雷諾方程進行歸一化處理，包括求解域簡化為矩形區(qū)域，并對雷諾方程進行進一步推導和簡化，構建等幾何分析模型；結合幾何多重網(wǎng)格方法要求，對IGA 法中的細分方法的研究，建立了基于h 細化的多重網(wǎng)格求解模型，提出了基于h細化的多重網(wǎng)格映射矩陣算法；最后引入活塞-缸套潤滑系統(tǒng)和滑動軸承兩個數(shù)值算例，選擇W循環(huán)和V循環(huán)兩種循環(huán)方式，用高斯-塞德爾迭代法進行求解，并與簡單IGA法計算效率進行比較。

2 求解雷諾方程的等幾何法

2.1 B 樣條與 NURBS 基函數(shù)

與有限元法類似，等幾何法里面也有“形函數(shù)”，即為樣條幾何自帶的NURBS 基函數(shù)。但等幾何法與有限元法也有明顯區(qū)別，NURBS 包括節(jié)點矢量、控制點向量和權重等要素，而有限元里面只有節(jié)點。有限元法是在節(jié)點上進行數(shù)值計算，等幾何法是在控制點上進行數(shù)值計算，等幾何法與有限元法的基函數(shù)的形式也有很大差別。這里首先對樣條基函數(shù)進行簡要介紹。首先給出節(jié)點向量：

式中：p—基函數(shù)的階數(shù)；n—基函數(shù)和控制點的個數(shù)。B樣條基函數(shù)可根據(jù)節(jié)點矢量定義為[6]：

與權重結合得到非均勻有理B樣條（NURBS），使用B樣條基函數(shù)Ni，p（ξ）、Ni，p（η）和權重wi，j可以獲得二維 NURBS 基函數(shù)其表達式為：

2.2 弱形式的雷諾方程

基于Navier-Stokes 方程和質(zhì)量連續(xù)性方程，可以在做出一些假設后建立雷諾方程[7，8]如下：

對方程（3）進行變量替換和歸一化后，在域D=［0，1］×［0，1］上運用伽遼金法并進行推導獲得雷諾方程的弱形式：

2.3 剛度矩陣推導

在獲得雷諾方程的弱形式之后，對雷諾方程的剛度矩陣進行推導。令C=diag（k1，k2），有：

假設物理域D由幾何域映射過IGA中NURBS 函數(shù)定義為：

對于在參數(shù)域D0=［0，1］2上的雙變量 NURBS 基函數(shù)Rij，以及cij∈R2上的控制點，物理域D上的被積函數(shù)g（X，Y）的積分可以轉換為參數(shù)域上的積分，其積分形式為：

這里，DF是 2×2 雅可比矩陣，基于P（X，Y）=P（F（u，v））的鏈式規(guī)則，微分可以寫為：

應用上述變換規(guī)則，式（4）左右兩端可以分別轉化為：

數(shù)值解P可以通過NURBS 基函數(shù)以及控制點qij∈R近似表示成如下形式：

通過將式（11）代入式（9）中，并用v=Rlk將v進行替換，對于L=1，…，m；K=1，…，n有：

其中，De是參數(shù)域的單元，并且獲得線性方程組Aq=b，其中這里，單元剛度矩陣Ae=（a（RI，RJ））表示為：

對于I，J=1，…，mn右側向量b=（＜f，RI＞）的組成部分為：

在前文描述的雷諾方程式中，物理域D已被歸一化為［0，1］×［0，1］。所以幾何域映射定義為F這里的參數(shù)域D0=［0，1］2，矩陣DF（u，v）為雅可比矩陣。對式（15）進行適當?shù)耐茖Ш螅瑒偠染仃嘇可以計算如下：

將剛度矩陣進行組裝，定義：

推導出的剛度矩陣為：

式（18）中we，p，q的是單元De上的高斯積分點處的權重，式（17）中的Ue，p，q和Ve，p，q依賴于 NURBS 基函數(shù)和對應的高斯積分點，c11，e，p，q和c22，e，p，q是從高斯積分點處根據(jù)油膜厚度 h計算的系數(shù)。

3 適于等幾何的多重網(wǎng)格法

3.1 多重網(wǎng)格法基本流程和方法

3.1.1 多重網(wǎng)格法的基本思想

上文的推導知，與有限元法類似，等幾何法離散之后會得到形如線性方程組Au=f的求解形式，只是剛度矩陣A和右邊項f的計算方法不同。求解線性方程組的方法主要有兩類一類是直接法，另一類是迭代法，對于雷諾方程來說，雖然等幾何法能以比有限元法更少的自由度來表達油膜壓力，但有時候為了達到更高的求解精度，自由度數(shù)也會比較多，求解的剛度矩陣也會隨之變大，當剛度矩陣A比較小且矩陣里面非零元素比較少時，適合用直接法進行求解。對于剛度矩陣維度很大的線性方程組，直接法的求解效率就會很低，這種條件下就需要用迭代法進行求解。在實際的工程問題偏微分方程求解中，大部分都采用迭代法進行求解。其求解公式如下：

收斂速度是數(shù)值計算方法的一個重要因素，等幾何法也不例外。為了滿足更高的精度要求，需要較多的自由度來表達油膜壓力，但這樣需要更多的迭代法收斂比較慢。為了加快收斂速度，在IGA 求解雷諾方程過程中引入多重網(wǎng)格法。迭代過程中產(chǎn)生的誤差分為高頻誤差和低頻誤差，當網(wǎng)格密度確定后，高頻誤差會在迭代過程中快速降低，而低頻誤差的消除較慢。而在網(wǎng)格密度變粗以后，細網(wǎng)格的低頻誤差在粗網(wǎng)格里面頻率就顯得相對較高，所以粗網(wǎng)格誤差收斂速度相對更快。多重網(wǎng)格法再細網(wǎng)格上滿足精度要求，用粗網(wǎng)格實現(xiàn)快速收斂。

3.1.2 多重網(wǎng)格法的基本流程

重網(wǎng)格法通過構建不同網(wǎng)格之間的映射關系進而加速計算，其中包含由粗到細的延拓矩陣R和細到粗的限制矩陣P。依據(jù)粗網(wǎng)格剛度矩陣的來源不同分為兩類，一類是幾何多重網(wǎng)格法，一類是代數(shù)多重網(wǎng)格法，它們的求解流程是基本相同。其中代數(shù)多重網(wǎng)格法對映射矩陣給出了如下限制：

而幾何多重網(wǎng)格則是基于網(wǎng)格層次關系的不同而構建出來的。這里采用的是幾何多重網(wǎng)格法，其求解方法將在下文進行詳細說明。多重網(wǎng)格法的基本思路是利用粗網(wǎng)格消除細網(wǎng)格上的低頻誤差，以達到更快收斂的目的。這里以V循環(huán)為例對多重網(wǎng)格法基本流程進行表述：

該方法用三層控制點網(wǎng)格對雷諾方程進行求解，第一層較細網(wǎng)格上的求解結果作為最終計算結果，用第二層較粗控制點網(wǎng)格對第一層細網(wǎng)格的求解進行加速，用第三層最粗網(wǎng)格對第二層網(wǎng)格上的計算進行進一步的加速。W循環(huán)類似。

從上面可以看出，求解過程需要兩個映射關系，一個是將rh從Ωh映射到Ω2h得到r2h的映射關系（即限制矩陣P），另一個是將e2h從Ω2h映射到Ωh上得到eh的映射關系（延拓矩陣R）。這是多重網(wǎng)格方法計算的關鍵，能否得到這個映射關系直接決定了方法的可行性。

3.2 多重網(wǎng)格的生成

基于有限元的多重網(wǎng)格法的映射關系局限在單元內(nèi)，細網(wǎng)格上節(jié)點的待解物理量可由粗網(wǎng)格上相應單元上的相關節(jié)點線性地表達，粗網(wǎng)格上要求解的節(jié)點的物理量也是由細網(wǎng)格上相應單元節(jié)點線性表示。而等幾何分析方法中的基函數(shù)對應的作用范圍是跨單元的，所以對應的多重網(wǎng)格間的映射關系也需要是跨單元的。由于等幾何是在控制點上進行數(shù)值計算，所以網(wǎng)格上對應的網(wǎng)格點是NURBS 里面的控制點。控制點的坐標由節(jié)點的坐標進行表達，所以每一次細化都是先在原來的節(jié)點向量上插入一些節(jié)點形成新的節(jié)點向量，然后再根據(jù)新的節(jié)點向量和原來的控制點和節(jié)點向量的值計算得到新的控制點。等幾何法中網(wǎng)格細化方法有h細化、p細化和k細化，由于經(jīng)p細化和k細化都需要對基函數(shù)進行升階，而傳統(tǒng)多重網(wǎng)格法中不同網(wǎng)格上沒有階次的差別，因此文中選用基于h細化的方式對網(wǎng)格進行細化，故有必要分析一下等幾何分析方法中的h細分算法。這里先以一維為例，介紹一下h細分方法。設已給一條k次B樣條曲線：

其中，B 樣條基函數(shù)只由節(jié)點向量 τ=［τ0，τ1，…，τn，τn+k+1］中的節(jié)點進行表達。如果在定義域內(nèi)某個節(jié)點區(qū)間內(nèi)插入一個節(jié)點t∈［τi，τi+1］?［τk，τn+1］以得到新的節(jié)點矢量T=［τ0，τ1，…，τi，t，τi+1，…，τn，τn+k+1］重新編號后成為T=［t0，t1，…，ti，ti+1，ti+2，…，tn+k+2］。然后根據(jù)新的節(jié)點矢量表達一組新的B樣條基函數(shù)0，1，…，n+1）。然后原樣條曲線可以用新的B樣條基函數(shù)與未知新頂點Di（j=0，1，…，n+1）共同表示為：

具體細分算法參考文獻[10]。

式中：Di—細網(wǎng)格控制點；Pi—粗網(wǎng)格控制點；i—粗節(jié)點向量上ti所在位置左邊的節(jié)點編號；k—基函數(shù)階數(shù)；t—細節(jié)點向量上節(jié)點值；τ—粗節(jié)點向量上節(jié)點值，r的值從k-1 開始隨著公式的遞推減1，直至r=1。

由于選擇的NURBS 基函數(shù)是三階的，故k=3，按上述公式可得，如式（24）所示。

這樣即可得到新的控制點向量D=［D1，D2，…，Dn+1］。

由于計算對象的求解域是二維的，所以在進行網(wǎng)格細化時分別要在兩個方向插入節(jié)點來進行網(wǎng)格細化（三維方向同理）。若需要l層網(wǎng)格，在最初網(wǎng)格上進行l(wèi)-1 次網(wǎng)格細化即可。

對于二維求解域，依次對兩個維度節(jié)點向量分別插入節(jié)點，得到的細化算法如下：

式中：V，V′—第二方向上的粗、細節(jié)點向量。

然后用式（25）計算方法對網(wǎng)格進行細化，一粗網(wǎng)格控制點網(wǎng)格和經(jīng)過一次細化后的細控制點網(wǎng)格圖，如圖1 所示。

圖1 控制點網(wǎng)格圖Fig.1 Control Point Grid Diagram

等幾何的控制點都是由求解域內(nèi)初始控制點進行細化得到。這種方法使多重網(wǎng)格的生成更加簡單，需要三層控制點網(wǎng)格，當進行l(wèi)次細化時，保存第l-2 次細化得到的控制點作為第三層，這一層控制點網(wǎng)格最粗，然后在進行兩次細化，分別得到的二層相對較粗的控制點網(wǎng)格和第一層最細控制點網(wǎng)格。得到三層控制點網(wǎng)格后，根據(jù)等幾何法求出各層網(wǎng)格上的剛度矩陣A1，A2，A3和最細網(wǎng)格上的右邊項f1。

3.3 多重網(wǎng)格法映射矩陣

與節(jié)點的細化方式對應，映射矩陣也是按照文獻[10]中的方法進行求解。根據(jù)式（23）推導出的式（24）的計算方法，當插入一個節(jié)點時，新的控制點向量和原控制點向量就會出現(xiàn)一個對應的映射關系R，此時新的控制點有n+1 個，而原控制點有n個。則新控制點向量和原控制點向量基于R的映射關系為：

以一個一維三次B樣條為例，其原有節(jié)點向量為τ=［0 0 0 0 1 1 1 1］，原有控制點向量為P=［0 0.333 0.6667 1.0000］。

當在原始節(jié)點向量中插入一個節(jié)點t=0.5，就會得到一個新的節(jié)點向量t=［0 0 0 0 0.5 1 1 1 1］，和一個新的控制點向量D＝［D1D2D3D4D5］，按照式（24）即可得到如下映射關系：D=D5×4P，此時的映射矩陣為：

而在h細化中，往往要插入的節(jié)點數(shù)目不只一個。在這種情況下，可依次單獨插入相應的節(jié)點來得到新的節(jié)點向量和控制點。例如要在節(jié)點向量 τ=［τ0，τ1，…，τn，τn+k+1］，中插入節(jié)點t1，t2，…，tm，此時依然可以按照式（24）得到新的控制點：

其中R（n+m）Xn即為粗網(wǎng)格到細網(wǎng)格的映射矩陣。

這里雷諾方程的求解域是二維的。其粗細網(wǎng)格間的映射矩陣可根據(jù)式（25）計算得到。

對細網(wǎng)格上每個控制點根據(jù)式（27）建立線性映射關系，即可組裝得到粗網(wǎng)格到細網(wǎng)格之間的整體映射矩陣R。由于細網(wǎng)格到粗網(wǎng)格的映射矩陣（即限制矩陣P）的求解比較困難，這里采用文獻ade 的方法求該映射矩陣，即映射矩陣P取R的轉置。

由式（24）可以發(fā)現(xiàn)，對于一維求解域，一個細網(wǎng)格控制點由4 個粗網(wǎng)格控制點線性表達，由式（25）可以看出，對于二維求解域，一個細網(wǎng)格控制點由16 個粗網(wǎng)格控制點線性表達，其映射關系，如圖2 所示。細網(wǎng)格向粗網(wǎng)格映射時，粗網(wǎng)格上的每個控制點也是由同樣數(shù)目的細網(wǎng)格控制點線性表達。粗、細網(wǎng)格上的映射關系是跨單元的，這樣就得到兩層控制點之間更好的線性關系。這種跨單元多對一的映射關系，有利于提高計算的穩(wěn)定性。

圖2 粗細網(wǎng)格控制點映射關系示意圖Fig.2 Schematic Diagram of Point Reflecting Relationship between Coarse and Fine Grid Control Points

4 數(shù)值算例

4.1 模型的數(shù)值求解

這里先分別用有限元法和等幾何法對活塞-缸套潤滑模型進行求解，通過各種自由度下兩種數(shù)值算法的油膜壓力積分值的對比來比較兩種解法的效果；然后分別用等幾何法和多重網(wǎng)格法多活塞-缸套和徑向滑動軸承兩種模型進行求解，并以全主元高斯消去法得到的值作為標準值，將得到的值與標準值相對誤差以10 為底求對數(shù)作為誤差參數(shù)e，求解與迭代次數(shù)對應的誤差值，繪制e-迭代次數(shù)折線圖。誤差參數(shù)e的表示方法，如式（28）所示。

其中，U=Ah＼fh，u是等幾何法、迭代法等得到的近似解。

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

通過對活塞缸之間的膜潤滑的雷諾方程進行求解得到的油膜壓力。活塞-缸套系統(tǒng)，如圖3（a）所示。其油膜潤滑區(qū)域[9]，如圖3（b）所示。

圖3 活塞--缸套潤滑系統(tǒng)示意圖Fig.3 Piston-Cylinder Liner Lubrication System

因為左右潤滑區(qū)域相對于平面φ=0 都是對稱的，所以只需考慮前面部分。通常，當活塞周期性地移動時，由于負載變化，密度和粘度的大小會發(fā)生變化。但為了簡化數(shù)值模型，這里將密度和粘度的值作為常數(shù)。給定油膜壓力的參數(shù)，如表1 所示。

表1 計算參數(shù)Tab.1 Calculate Parameters

假設v0=0，vh=v（其正方向為上行），潤滑劑無側漏，因此u0=uh=0。那么有：

將求解域D歸一化為［0，1］×［0，1］，進行如下轉換：

右側（因為對稱，只考慮前面部分）

（1）等幾何法與有限元法精度比較

分別用有限元法和等幾何法對活塞-缸套潤滑模型油膜壓力進行求解，并將其比較兩種算法得到的壓力積分隨自由度變化的趨勢。其結果，如圖4 所示。

圖4 不同自由度時有限元法和等幾何法壓力積分圖Fig.4 Pressure Integral of FEM and Isogeometric Under Different Degree of Freedom

根據(jù)圖4 可以看出，隨著自由度的增加，兩種數(shù)值算法都更趨近于收斂值，但等幾何法趨近的速度更快。因此可以證明在自由度相等的情況下，等幾何法的精度高于有限元法。這種特性綜合了有限差分法和有限元法的速度與精度的優(yōu)勢，非常適合流體仿真分析。

（2）將雅可比迭代、高斯塞德爾迭代、超松弛迭代三種代法與等幾何多重網(wǎng)格法的求解效率一起進行對比。超松弛迭代松弛因子的值根據(jù)最佳松弛因子選取。求解得到的油膜壓力圖和r-迭代次數(shù)趨勢圖，如圖5、圖6 所示。

圖5 缸套活塞油膜壓力Fig.5 Oil Film Pressure of Cylinder Piston

圖6 缸套活塞e-迭代次數(shù)圖Fig.6 e-Iteration Number Diagram For Cylinder Liner Piston

從圖6 可以看出，在迭代次數(shù)從25 增加到400 的過程中，高斯塞德爾迭代和雅可比迭代收斂速度都最慢，SOR 的收斂效率相對較快，多重網(wǎng)格法收斂速度最快。在迭代次數(shù)為25 時，SOR、高斯塞德爾迭代和雅可比迭代的誤差相差不大，多重網(wǎng)格法誤差相對較小，在迭代次數(shù)為400 時，高斯塞德爾迭代和雅可比迭代的誤差接近e-2，SOR 迭代誤差在e-6 到e-7 之間，而多重網(wǎng)格法迭代誤差接近e-10。因此用高斯塞德爾迭代和雅可比迭代時，單純等幾何法收斂速度很慢，在改用SOR 進行迭代后，速度明顯加快，但仍然明顯比等幾何多重網(wǎng)格法收斂速度慢。多重網(wǎng)格法在迭代過程中收斂速度曲線出現(xiàn)了拐點，在拐點后比拐點前的趨勢更加平緩，此處是由于限制矩陣P是通過取延拓矩陣R的轉置得到的，而不是進行理論計算得到。這樣會使限制矩陣的映射有一定偏差，由于剛開始迭代時得到的數(shù)值解與真實值的誤差較大，這種偏差的影響較小，但當誤差較小時，這種偏差的影響就顯得很明顯。

4.3 徑向滑動軸承模型

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

圖7 徑向滑動軸承潤滑模型Fig.7 Radial Plain Bearing Lubrication Model

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

表2 計算參數(shù)Tab.2 Calculate Parameters

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

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

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

將雅可比迭代、高斯塞德爾迭代、超松弛迭代三種代法與等幾何多重網(wǎng)格法的求解效率一起進行對比。超松弛迭代松弛因子的值依據(jù)最佳松弛因子選取。求解得到的油膜壓力圖和r-迭代次數(shù)趨勢圖，如圖8、圖9 所示。

圖8 徑向滑動軸承油膜壓力Fig.8 Oil Film Pressure of Radial Sliding Bearing

圖9 徑向滑動軸承e-迭代次數(shù)圖Fig.9 e-Iteration Number Diagram of Radial Sliding Bearing

從圖9 可以看出可以得到與缸套活塞潤滑模型相似的結論。不同的是兩種模型里面高斯塞德爾迭代和雅可比迭代的收斂速度有較小差別，而且兩種模型的多重網(wǎng)格法迭代收斂曲線拐點出現(xiàn)的位置不同。這兩種區(qū)別的形成是由于兩種模型不同參數(shù)形成了不同的剛度矩陣的特征。

5 結論

在二次精度NURBS 樣條的基礎上，對雷諾方程進行推導，并對節(jié)點插入的細分算法進行研究，提出基于h細化的細分算法，并建立適于等幾何多重網(wǎng)格法的求解模型，提出基于h細化的映射矩陣求解方法。引入缸套活塞潤滑系統(tǒng)和徑向滑動軸承兩個數(shù)值算例，先用單純等幾何法和有限元法對活塞-缸套潤滑模型進行求解，比較不同自由度條件下的壓力積分。計算結果證明等幾何法的求解精度明顯高于有限元法。在此基礎之上，分別通過基于單純等幾何法的雅可比、高斯塞德爾和超松弛迭代三種代法和等幾何多重網(wǎng)格法對活塞-缸套潤滑模型和徑向滑動軸承潤滑模型油膜壓力進行求解，比較其收斂速度。

根據(jù)計算結果可以得出如下結論，相比于有限元法，等幾何法求解雷諾方程能以較小的自由度實現(xiàn)更高的精度。超松弛迭代的收斂效率高于雅克比迭代和高斯塞德爾迭代，但引入多重網(wǎng)格法進行求解后，收斂效率明顯高于超松弛迭代、雅克比迭代和高斯塞德爾迭代；由于限制矩陣P是通過取延拓矩陣R的轉置的方法得到而造成映射偏差，造成多重網(wǎng)格法收斂曲線出現(xiàn)拐點；對于相同的迭代方法，兩種模型的收斂趨勢有較小的區(qū)別，是由于不同模型剛度矩陣的特征不同。