諸永定 王 威 王 軍 李 斌 尹國慶

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

0 引言

在葉輪模型的幾何表達時，通過參數化造型及擬合，可以將一個復雜的葉輪幾何模型用若干個簡單的控制參數來表達，這樣避免了CAD 模型中復雜的三維曲面和離散點模型中龐大的離散點數據，同樣也便于改變控制參數來實現葉輪的改型。

CFD 數值分析方法手段是現代葉輪機械設計的特征，其分析精度與效果評價與其內流流道網格劃分密切相關。

由于葉輪機流道及葉片形狀具有高度三維性以及內部流動的復雜性，網格質量的正交性以及網格在流場中大梯度區(qū)域的密度對數值計算結果影響很大，因此高質量的網格生成技術往往是提高求解效率并正確模擬復雜流動的關鍵。

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

1 葉輪的參數化造型

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

1.1 葉片成型的積疊線的參數化擬合

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

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

1.2 二維葉型型線的參數化擬合

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

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

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

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

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

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

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

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

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

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

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

1.3 子午通道型線參數化擬合

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

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

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

通過分別對葉輪葉片和子午流道型線的參數化擬合，最終實現對Rotor 67的參數化擬合，如圖5所示。

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

2 葉片網格生成

本文采用代數方法生成葉片計算域的H 型結構網格，根據邊界節(jié)點通過直接的代數計算與插值得到通道中計算區(qū)域內結構網格節(jié)點分布，生成分塊結構化網格。

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

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

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

圖7 葉片表面網格分布Fig.7 Grid in the blade surface

將生成的子午流面上的節(jié)點根據設定的周向節(jié)點分布進行周向旋轉，可以得到整個流道網格點的軸向坐標、周向坐標和徑向坐標()z,r,x，圖8 給出了三個流面上的網格生成情況。

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

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

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

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

為了保證連接處網格的銜接質量，讓控制參數α滿足給定第一個網格間距的已知條件：

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

圖9 為加密后葉片流道的H 型網格。采用網格自動生成算法生成網格后，需要對所生成的網格進行質量檢測，評價其有效性。本文選用雅各比比率，即最小雅各比矩陣與最大雅各比矩陣行列式的比值作為六面體網格單元的質量評價標準[6]。單元的雅各比矩陣反映了網格的基本質量指標，包括角度、長度和面積等。

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

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

圖9 Rotor 67計算網格Fig.9 Computional grid of Rotor 67

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

3 計算結果與討論

為了驗證參數化擬合精度以及網格生成質量，生成的網格數據文件導入Fluent 中進行數值計算。Strazisar[7]等人對Rotor 67 進行了詳細、高精度的熱態(tài)幾何、流場結構和整機性能試驗測試，因而作為流動機理研究和葉輪機CFD校核的標準案例而廣為使用。Rotor 67的基本設計參數如表1所示。

表1 Rotor 67基本設計參數Tab.1 Rotor 67 design parameters

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

表2 網格節(jié)點分布系數值Tab.2 Coefficients of node distribution

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

圖10 網格質量檢查Fig.10 Grid quality checking

采用密度基隱式求解，二階離散格式。湍流模型采用標準SSTk-ω模型，壁面為無滑移邊界條件。計算中，與試驗條件一致，設置進口總壓為101325Pa、進口總溫為288.15K。

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

圖11 Rotor 67設計轉速下的特性曲線與實驗對比Fig.11 Comparison of simulated and measured overall performance curves for Rotor 67

圖12給出了最高效率點附近不同葉高截面相對馬赫數的等值線圖以及相應的實驗結果。經過對比發(fā)現，流場中的相對馬赫數與實驗結果符合較好，能夠很好地描述流場細節(jié)。

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

4 結論

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

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

3）用Rotor 67 轉子實例的參數化擬合以及網格劃分的計算結果驗證本文方法的可行性，擬合精度以及網格生成質量滿足計算分析要求。