諸永定 王 威 王 軍 李 斌 尹國(guó)慶

（1.寧波方太廚具有限公司；2.華中科技大學(xué)能源與動(dòng)力工程學(xué)院）

0 引言

在葉輪模型的幾何表達(dá)時(shí)，通過參數(shù)化造型及擬合，可以將一個(gè)復(fù)雜的葉輪幾何模型用若干個(gè)簡(jiǎn)單的控制參數(shù)來表達(dá)，這樣避免了CAD 模型中復(fù)雜的三維曲面和離散點(diǎn)模型中龐大的離散點(diǎn)數(shù)據(jù)，同樣也便于改變控制參數(shù)來實(shí)現(xiàn)葉輪的改型。

CFD 數(shù)值分析方法手段是現(xiàn)代葉輪機(jī)械設(shè)計(jì)的特征，其分析精度與效果評(píng)價(jià)與其內(nèi)流流道網(wǎng)格劃分密切相關(guān)。

由于葉輪機(jī)流道及葉片形狀具有高度三維性以及內(nèi)部流動(dòng)的復(fù)雜性，網(wǎng)格質(zhì)量的正交性以及網(wǎng)格在流場(chǎng)中大梯度區(qū)域的密度對(duì)數(shù)值計(jì)算結(jié)果影響很大，因此高質(zhì)量的網(wǎng)格生成技術(shù)往往是提高求解效率并正確模擬復(fù)雜流動(dòng)的關(guān)鍵。

本文選取NASA Rotor67 作為研究對(duì)象，分別采用CST 和3 次B 樣條參數(shù)化方法來擬合Rotor67 二維葉型型線、積疊線以及子午通道型線。通過對(duì)控制點(diǎn)的反算求解，最終實(shí)現(xiàn)對(duì)Rotor67 的參數(shù)化擬合。對(duì)于參數(shù)化擬合后的葉輪機(jī)葉片，采用代數(shù)方法生成H型網(wǎng)格，通過給定的邊界上的節(jié)點(diǎn)分布，插值得到計(jì)算區(qū)域內(nèi)網(wǎng)格節(jié)點(diǎn)分布，并采用雅各比比率（Jacoby ratio）對(duì)生成的網(wǎng)格質(zhì)量進(jìn)行評(píng)價(jià)。為了驗(yàn)證參數(shù)化擬合的精確性以及網(wǎng)格的質(zhì)量，本文計(jì)算了設(shè)計(jì)工況下Rotor67 特性曲線并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比對(duì)，論證了方法的可行性。

1 葉輪的參數(shù)化造型

通常原始葉片型面是由多個(gè)截面二維葉型型線按一定的規(guī)律積疊而成的空間曲面[1]，因此本文參數(shù)化擬合主要思路是先參數(shù)化二維葉型，再通過積疊線構(gòu)成一個(gè)三維葉片方式完成。

1.1 葉片成型的積疊線的參數(shù)化擬合

積疊線的參數(shù)化是為了更好的控制彎、掠規(guī)律，并探索新型的三維復(fù)雜造型。為計(jì)算方便，本文選擇二維葉型的前緣點(diǎn)作為積疊點(diǎn)，葉片的彎掠規(guī)律的控制便是通過控制積疊點(diǎn)形成的積疊線的形狀來實(shí)現(xiàn)的。葉片積疊線屬于三維空間曲線，將其分為切向彎和前后掠兩條二維曲線分別進(jìn)行參數(shù)化擬合，同樣采用3 次B 樣條進(jìn)行參數(shù)化擬合。本文采用14 個(gè)型值點(diǎn)擬合出Rotor67 的參數(shù)化彎、掠積疊線，如圖1 所示。

圖1 B樣條參數(shù)化彎、掠積疊線Fig.1 Parameterised shape of stacking line by B-Spline method

1.2 二維葉型型線的參數(shù)化擬合

葉型的位置和范圍確定后，即可以構(gòu)造葉型型線并確定其在流面上的位置。對(duì)于一般的二維葉型，應(yīng)根據(jù)所給離散點(diǎn)，擬合出二維葉型型線。葉型型線的擬合，越接近實(shí)際葉型形狀，對(duì)后面的流場(chǎng)計(jì)算產(chǎn)生的誤差就越小，同時(shí)為了減少氣流經(jīng)過葉柵時(shí)的損失，要求葉柵型面光滑。由于CST 參數(shù)化方法具有設(shè)計(jì)變量少，可調(diào)節(jié)，設(shè)計(jì)空間廣等優(yōu)點(diǎn)[2]，本文采用CST 參數(shù)化方法對(duì)二維葉型進(jìn)行擬合。

CST參數(shù)化下的通用葉型可以表示如下[3]：

N1和N2定義了所表示幾何外形具有的基本幾何特性，對(duì)于雙圓弧翼型取N1=1，N2=1。形狀函數(shù)S(ψ)通常采用n階Bernstein 多項(xiàng)式的加權(quán)和作為表達(dá)式，即：

式中，n為多項(xiàng)式中的指數(shù)；i為多項(xiàng)式階數(shù)；bi(i=0,1,2,…n)為Bernstein系數(shù)。

對(duì)于三維葉片可由多個(gè)葉型型面組成，如圖2 所示。因此可定義三維葉片的方程為：

圖2 CST 葉型截面定義Fig.2 Blade sectional definition by CST method

由以上分析可知，葉型型線的CST 參數(shù)化方法的控制參數(shù)為類函數(shù)中的N1,N2，形狀函數(shù)中的Bernstein系數(shù)b0,b1…bn以及翼型后緣厚度ζT。為了計(jì)算Bernstein 系數(shù)，需要得到一組控制點(diǎn)，由于Bernstein多項(xiàng)式的極值點(diǎn)均勻分布在翼型的弦上：

代入式（4）即是所需要的一組控制點(diǎn)[ψi,ζ(ψi)]，因此有：

從而可以求解出一組CST參數(shù)b0,b1…bn。

圖3展示了Rotor67不同截面位置葉型的CST翼型擬合效果。

圖3 Rotor 67葉型的CST擬合效果Fig.3 Rotor 67 profile parameterised CST method

1.3 子午通道型線參數(shù)化擬合

子午通道型線包括輪轂線和機(jī)匣線，其形狀對(duì)葉輪機(jī)的工作性能有重要影響，合理的子午面流道變化不僅能有效改變?nèi)~輪機(jī)的流量與壓比，還能在一定程度上抑制泄漏渦發(fā)展，減小其造成的損失，使得總體性能獲得提升。

圖4 B樣條參數(shù)化子午通道型線Fig.4 Parameterised shape of meridional plane by BSpline method

NASA Rotor 67 的子午通道型線如圖4 所示，可以看出輪轂沿流向逐漸增大，機(jī)匣沿流向逐漸減小，由于影響葉輪機(jī)械性能的主要因素在輪轂和機(jī)匣的下壓與上抬[4]，因此將控制參數(shù)選擇在徑向變化量較大的地方，其它數(shù)據(jù)點(diǎn)根據(jù)插值得到。為了比較好的擬合端壁型線，并且比較少地采用擬合參數(shù)，本文選擇5 個(gè)控制點(diǎn)的三次B樣條曲線分別擬合輪轂線和機(jī)匣線。圖4為Rotor 67子午通道型線參數(shù)化擬合效果。

通過分別對(duì)葉輪葉片和子午流道型線的參數(shù)化擬合，最終實(shí)現(xiàn)對(duì)Rotor 67的參數(shù)化擬合，如圖5所示。

圖5 Rotor 67參數(shù)化擬合效果Fig.5 Parameterised 3D model of Rotor 67

2 葉片網(wǎng)格生成

本文采用代數(shù)方法生成葉片計(jì)算域的H 型結(jié)構(gòu)網(wǎng)格，根據(jù)邊界節(jié)點(diǎn)通過直接的代數(shù)計(jì)算與插值得到通道中計(jì)算區(qū)域內(nèi)結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)分布，生成分塊結(jié)構(gòu)化網(wǎng)格。

先計(jì)算葉片區(qū)子午面上網(wǎng)格點(diǎn)坐標(biāo)()z,r。通過上節(jié)參數(shù)化得到的葉型截面數(shù)據(jù)和輪轂、機(jī)匣的數(shù)據(jù)，用代數(shù)方法求得葉片前、尾緣與輪轂、機(jī)匣的交點(diǎn)，然后根據(jù)預(yù)先設(shè)置的葉片區(qū)軸向與徑向的節(jié)點(diǎn)分布設(shè)置，可以得到葉片區(qū)子午面上邊界的節(jié)點(diǎn)分布，進(jìn)而采用代數(shù)插值得到葉片區(qū)子午面上網(wǎng)格分布[5]，如圖6所示。同樣的方法可以得到葉片區(qū)上游和下游子午面上的網(wǎng)格。

圖6 葉片區(qū)子午面網(wǎng)格分布Fig.6 Grid in the meridional plane

將葉片區(qū)子午面上的節(jié)點(diǎn)信息與葉型截面數(shù)據(jù)聯(lián)立求解，可以分別得到葉片吸力面與壓力面的網(wǎng)格節(jié)點(diǎn)分布，如圖7所示：

圖7 葉片表面網(wǎng)格分布Fig.7 Grid in the blade surface

將生成的子午流面上的節(jié)點(diǎn)根據(jù)設(shè)定的周向節(jié)點(diǎn)分布進(jìn)行周向旋轉(zhuǎn)，可以得到整個(gè)流道網(wǎng)格點(diǎn)的軸向坐標(biāo)、周向坐標(biāo)和徑向坐標(biāo)()z,r,x，圖8 給出了三個(gè)流面上的網(wǎng)格生成情況。

圖8 S1，S2，S3流面的生成網(wǎng)格Fig.8 Grid in the S1,S2,S3 stream surface

在流動(dòng)計(jì)算的網(wǎng)格生成中，在保持某一方向上網(wǎng)格總數(shù)不變的前提下，為了捕捉局部地區(qū)高梯度的物理量變化，須要加密局部網(wǎng)格，同時(shí)為了保證計(jì)算的精度，還須使網(wǎng)格的跨度光滑過渡。對(duì)于周向網(wǎng)格和徑向網(wǎng)格節(jié)點(diǎn)分布要求兩端密而中間稀，本文采用如下雙曲正切關(guān)系式來控制柵距內(nèi)周向網(wǎng)格及輪轂和機(jī)匣的附面層徑向網(wǎng)格的疏密，即

式中，M為節(jié)點(diǎn)個(gè)數(shù)；s1，sM分別為起始與終止點(diǎn)；α，β分別為阻尼和伸展因子；當(dāng)α=1,β=1 時(shí)，網(wǎng)格節(jié)點(diǎn)均勻分布；當(dāng)α ＞1,β ＜1 時(shí)，網(wǎng)格節(jié)點(diǎn)在兩端對(duì)稱加密；α越大，β越小，則兩端節(jié)點(diǎn)密度越大。

同樣，對(duì)于葉輪進(jìn)口截面至葉柵通道進(jìn)口和葉柵通道出口至葉輪出口截面的網(wǎng)格線上的節(jié)點(diǎn)分布要求從一段向另一端節(jié)點(diǎn)分布由密到疏，本文采用下列形式來進(jìn)行節(jié)點(diǎn)分布，即

為了保證連接處網(wǎng)格的銜接質(zhì)量，讓控制參數(shù)α滿足給定第一個(gè)網(wǎng)格間距的已知條件：

其中，Δs=s2-s1；Δq=

圖9 為加密后葉片流道的H 型網(wǎng)格。采用網(wǎng)格自動(dòng)生成算法生成網(wǎng)格后，需要對(duì)所生成的網(wǎng)格進(jìn)行質(zhì)量檢測(cè)，評(píng)價(jià)其有效性。本文選用雅各比比率，即最小雅各比矩陣與最大雅各比矩陣行列式的比值作為六面體網(wǎng)格單元的質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)[6]。單元的雅各比矩陣反映了網(wǎng)格的基本質(zhì)量指標(biāo)，包括角度、長(zhǎng)度和面積等。

假設(shè)xk∈?3為任意六面體單元的第k個(gè)節(jié)點(diǎn)的方向矢量且xk=[xk,yk,zk]T,k=1,2,…8。xk,i∈?3(i=1,2,3) 表示與該節(jié)點(diǎn)相鄰的所有節(jié)點(diǎn)。ek,i=xk,i-xk(i=1,2,3) 為節(jié)點(diǎn)k和其相鄰節(jié)點(diǎn)的方向矢量按順序構(gòu)成三個(gè)邊矢。節(jié)點(diǎn)k的雅各比矩陣行列式的值Jk=det[ek,1ek,2ek,3]，將六面體單元上8 個(gè)節(jié)點(diǎn)的單位雅各比矩陣行列式的值Jk的最小值和最大值的比重定義為雅各比比率JR，即

其中，Jmin=min(Jk),Jmax=max(Jk)，如果一個(gè)單元Jmin≤0 則該單元不滿足計(jì)算要求，在0～+1 之間，Jmin越大，說明單元的質(zhì)量越好。

圖9 Rotor 67計(jì)算網(wǎng)格Fig.9 Computional grid of Rotor 67

為了便于氣動(dòng)熱力計(jì)算，需要將網(wǎng)格數(shù)據(jù)轉(zhuǎn)化為計(jì)算網(wǎng)格文件，對(duì)于一個(gè)網(wǎng)格文件，須包含網(wǎng)格節(jié)點(diǎn)、網(wǎng)格單元、網(wǎng)格域、網(wǎng)格邊界等網(wǎng)格信息信息。ANSYS求解器可以讀取多種形式的網(wǎng)格文件，本文采用擴(kuò)展名為cfx5 格式為網(wǎng)格文件，*.cfx5 網(wǎng)格文件是按照非結(jié)構(gòu)化網(wǎng)格的數(shù)據(jù)形式存儲(chǔ)計(jì)算網(wǎng)格的，即網(wǎng)格文件僅存儲(chǔ)一系列離散的網(wǎng)格單元，不直接存儲(chǔ)網(wǎng)格鄰接關(guān)系。它可以直接采用ASC-II字碼存儲(chǔ)網(wǎng)格，方便閱讀，同時(shí)也易于實(shí)現(xiàn)程序輸出。

3 計(jì)算結(jié)果與討論

為了驗(yàn)證參數(shù)化擬合精度以及網(wǎng)格生成質(zhì)量，生成的網(wǎng)格數(shù)據(jù)文件導(dǎo)入Fluent 中進(jìn)行數(shù)值計(jì)算。Strazisar[7]等人對(duì)Rotor 67 進(jìn)行了詳細(xì)、高精度的熱態(tài)幾何、流場(chǎng)結(jié)構(gòu)和整機(jī)性能試驗(yàn)測(cè)試，因而作為流動(dòng)機(jī)理研究和葉輪機(jī)CFD校核的標(biāo)準(zhǔn)案例而廣為使用。Rotor 67的基本設(shè)計(jì)參數(shù)如表1所示。

表1 Rotor 67基本設(shè)計(jì)參數(shù)Tab.1 Rotor 67 design parameters

將采用上述方法生成的網(wǎng)格導(dǎo)入Fluent 中進(jìn)行數(shù)值計(jì)算。本文設(shè)置單流道網(wǎng)格節(jié)點(diǎn)數(shù)為50×50×150（切向×徑向×軸向），表2為設(shè)置網(wǎng)格節(jié)點(diǎn)分布規(guī)律的系數(shù)值。α,β分別為式(7)中的阻尼和伸展因子，而進(jìn)出口段流向節(jié)點(diǎn)分布規(guī)律的系數(shù)α0值根據(jù)第一個(gè)網(wǎng)格間距迭代求得。

表2 網(wǎng)格節(jié)點(diǎn)分布系數(shù)值Tab.2 Coefficients of node distribution

網(wǎng)格質(zhì)量如圖10所示，由圖可知，計(jì)算域的網(wǎng)格雅各比比率在0.3 以上，70%以上的網(wǎng)格其雅各比比率均在0.7，網(wǎng)格正交性滿足數(shù)值計(jì)算要求。

圖10 網(wǎng)格質(zhì)量檢查Fig.10 Grid quality checking

采用密度基隱式求解，二階離散格式。湍流模型采用標(biāo)準(zhǔn)SSTk-ω模型，壁面為無滑移邊界條件。計(jì)算中，與試驗(yàn)條件一致，設(shè)置進(jìn)口總壓為101325Pa、進(jìn)口總溫為288.15K。

Rotor 67 在設(shè)計(jì)轉(zhuǎn)速下通過實(shí)驗(yàn)測(cè)得的堵塞流量為34.96kg/s，本文數(shù)值計(jì)算得到的堵塞流量為34.48kg/s，在實(shí)驗(yàn)測(cè)量結(jié)果的誤差范圍之內(nèi)。圖11是計(jì)算得到的Rotor 67 在設(shè)計(jì)轉(zhuǎn)速下的特性曲線與實(shí)驗(yàn)值的對(duì)比圖。從圖中可以看出，數(shù)值計(jì)算得到的壓氣機(jī)設(shè)計(jì)轉(zhuǎn)速下的效率與實(shí)驗(yàn)值相比具有較為一致的趨勢(shì)，但數(shù)值偏低，偏差基本在5%以內(nèi)；數(shù)值計(jì)算得到的不同流量下的總壓比和實(shí)驗(yàn)值相接近，且變化趨勢(shì)與實(shí)驗(yàn)相一致，計(jì)算具有滿意的可信度。

圖11 Rotor 67設(shè)計(jì)轉(zhuǎn)速下的特性曲線與實(shí)驗(yàn)對(duì)比Fig.11 Comparison of simulated and measured overall performance curves for Rotor 67

圖12給出了最高效率點(diǎn)附近不同葉高截面相對(duì)馬赫數(shù)的等值線圖以及相應(yīng)的實(shí)驗(yàn)結(jié)果。經(jīng)過對(duì)比發(fā)現(xiàn)，流場(chǎng)中的相對(duì)馬赫數(shù)與實(shí)驗(yàn)結(jié)果符合較好，能夠很好地描述流場(chǎng)細(xì)節(jié)。

圖12 最高效率點(diǎn)附近不同葉高處相對(duì)馬赫數(shù)分布等值線圖Fig.12 Relative Mach number contours line of different blade height near peak efficiency

4 結(jié)論

1）建立葉輪葉片三維實(shí)體的參數(shù)化造型一種方法，采用CST和3次B樣條參數(shù)化方法完成葉片的參數(shù)化擬合，運(yùn)用代數(shù)方法建立葉輪葉片六面體網(wǎng)格生成。

2）開發(fā)前處理計(jì)算程序，可以高效、快速地對(duì)葉片進(jìn)行參數(shù)化擬合并對(duì)流道進(jìn)行網(wǎng)格劃分，通過邊界曲線和邊界曲面上的網(wǎng)格節(jié)點(diǎn)分布控制整各流道網(wǎng)格，輸出網(wǎng)格數(shù)據(jù)能夠被調(diào)用，從而能快速完成數(shù)值模擬計(jì)算，有效縮短CFD前處理時(shí)間，為葉片的參數(shù)優(yōu)化設(shè)計(jì)打下基礎(chǔ)。

3）用Rotor 67 轉(zhuǎn)子實(shí)例的參數(shù)化擬合以及網(wǎng)格劃分的計(jì)算結(jié)果驗(yàn)證本文方法的可行性，擬合精度以及網(wǎng)格生成質(zhì)量滿足計(jì)算分析要求。