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

李子強(qiáng) ，羅會信 ，左兵權(quán) ，張 希

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

1 引言

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

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

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

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

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

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

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

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

2.2 弱形式的雷諾方程

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

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

2.3 剛度矩陣推導(dǎo)

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

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

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

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

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

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

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

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

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

在前文描述的雷諾方程式中，物理域D已被歸一化為［0，1］×［0，1］。所以幾何域映射定義為F這里的參數(shù)域D0=［0，1］2，矩陣DF（u，v）為雅可比矩陣。對式（15）進(jìn)行適當(dāng)?shù)耐茖?dǎo)后，剛度矩陣A可以計(jì)算如下：

將剛度矩陣進(jìn)行組裝，定義：

推導(dǎo)出的剛度矩陣為：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

具體細(xì)分算法參考文獻(xiàn)[10]。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 數(shù)值算例

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

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

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

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

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

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

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

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

假設(shè)v0=0，vh=v（其正方向?yàn)樯闲校瑵櫥瑒o側(cè)漏，因此u0=uh=0。那么有：

將求解域D歸一化為［0，1］×［0，1］，進(jìn)行如下轉(zhuǎn)換：

右側(cè)（因?yàn)閷ΨQ，只考慮前面部分）

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

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

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

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

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

4.3 徑向滑動軸承模型

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

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

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

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

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

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

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

將雅可比迭代、高斯塞德爾迭代、超松弛迭代三種代法與等幾何多重網(wǎng)格法的求解效率一起進(jìn)行對比。超松弛迭代松弛因子的值依據(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 可以看出可以得到與缸套活塞潤滑模型相似的結(jié)論。不同的是兩種模型里面高斯塞德爾迭代和雅可比迭代的收斂速度有較小差別，而且兩種模型的多重網(wǎng)格法迭代收斂曲線拐點(diǎn)出現(xiàn)的位置不同。這兩種區(qū)別的形成是由于兩種模型不同參數(shù)形成了不同的剛度矩陣的特征。

5 結(jié)論

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

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