縱列式雙旋翼直升機氣動力和流場的數(shù)值模擬

王彩楓，楊永飛，林永峰，黃明恪

(1.中國直升機設(shè)計研究所，江西 景德鎮(zhèn) 333001;2.南京航空航天大學(xué)航空宇航學(xué)院，江蘇南京 210016)

0 引言

縱列式直升機旋翼的尾跡受前飛自由流、機身和旋翼之間的干擾影響，旋翼尾跡不僅影響機身的空氣動力，也影響旋翼自身的空氣動力特性。因此，旋翼之間以及旋翼和機身的相互干擾一直是直升機空氣動力學(xué)中的重要課題。經(jīng)典的旋翼渦流理論，無論是固定尾跡還是自由尾跡模型，都不能考慮機身的影響，因為它們都是以無旋位流理論為基礎(chǔ)，用線渦系來模擬尾跡。歐拉方程是精確的無粘氣體動力學(xué)方程組，可以模擬有旋流動。用歐拉方程計算無需專門考慮翼面后的尾跡旋渦，它們會自動地包含在計算結(jié)果之中。對直升機來說，能比較容易地考慮旋翼之間以及旋翼和機身的相互干擾。

本文應(yīng)用計算流體力學(xué)(CFD)方法，采用自適應(yīng)非結(jié)構(gòu)直角網(wǎng)格和歐拉方程計算雙旋翼直升機全機繞流。直升機繞流多為非定常流，但是，如果采用非定常來處理，計算時間很長，難以在工程上應(yīng)用。本文采用的方法是對直升機槳盤平面進(jìn)行時間和空間的平均［1-4］，從而將非定常問題轉(zhuǎn)化為定常問題。將槳盤平面看作是氣流可穿透的壓力不連續(xù)面［5］，采用CFD方法計算流場。由計算得到的槳盤平面位置的流動速度，再用葉素理論計算作用在旋翼上的氣動力，然后將旋翼上的氣動力進(jìn)行時間和空間的平均，轉(zhuǎn)化為槳盤平面上的壓力差，再轉(zhuǎn)回到CFD計算，從而構(gòu)成迭代循環(huán)(見圖1)。

圖1 計算流程圖

1 網(wǎng)格生成

以某一種重型縱列式直升機模型為例，說明網(wǎng)格生成方法。

網(wǎng)格生成時將直升機的幾何外形表面劃分為許許多多的三角形，并認(rèn)為直升機表面的三角形為平面。確定流場計算區(qū)域，構(gòu)成較粗的第一層直角網(wǎng)格。找出與物面相碰的單元，將相碰的單元及其鄰居單元等分為八個單元，如此形成第二層加密網(wǎng)格。再對第二層網(wǎng)格做上述自適應(yīng)加密，形成第三層加密網(wǎng)格。繼續(xù)上述過程，從而生成空間自適應(yīng)八叉樹非結(jié)構(gòu)直角網(wǎng)格。將上面生成的靠近物面的網(wǎng)格進(jìn)行最后的光順，形成最終的貼體非結(jié)構(gòu)網(wǎng)格(見圖2)。

圖2 縱列式直升機的最終表面網(wǎng)格

2 基于歐拉方程的CFD分析方法

1)數(shù)學(xué)模型

歐拉方程是無粘氣體動力學(xué)方程組，它包含連續(xù)方程、動量方程和能量方程，可寫為:

其中xj是n維實數(shù)空間Rn的坐標(biāo)分量，1≤j≤n，解域 Ω?Rn×［0，t)，n=3。W是守恒變量，F(xiàn)j為流通量矢量，其表達(dá)式分別為:

理想氣體的狀態(tài)方程為:

其中p是壓強，γ為比熱比，對于空氣，γ=1.4。

其中α是來流攻角，M∞是來流馬赫數(shù)。

2)有限體積法

應(yīng)用Jameson的有限體積法及四步Runge－Kutta時間推進(jìn)，此外，還引入了當(dāng)?shù)貢r間步長，焓阻擬等技術(shù)，提高收斂速度，減少了計算時間。在遠(yuǎn)場邊界處理上采用無反射邊界條件。

對于三維情形，將Euler方程在任一有限單元上積分，并應(yīng)用格林定理得到積分型歐拉方程，有:

其中Ωk表示編號為K的網(wǎng)格單元，?Ωk為其邊界，則在網(wǎng)格單元面積足夠小的情況下，近似地，有:

其中Ak為網(wǎng)格單元面積，Wk可視作網(wǎng)格中心處的守恒變量。

從該網(wǎng)格中流出的凈通量為:

式中為面積向量，其大小為面積，方向指向面的外法線，M是多面體的面數(shù)。網(wǎng)格交界處的通量Fi的值，用W在相鄰網(wǎng)格中心處的值的平均值來計算:

其中W左、W右分別為面的左右兩側(cè)網(wǎng)格中心處的守恒量。

最終得到半離散化方程為:

后面將要對該式沿時間方向積分，最終得到不隨時間變化的定常解。

3)遠(yuǎn)場邊界處理

遠(yuǎn)場邊界條件數(shù)目需根據(jù)邊界附近信息傳播的性質(zhì)確定，取一維Riemann不變量型無反射邊界條件，可分為四種情形。由于流場外部自由流信息是沿特征線傳入流場內(nèi)的，所以有:

①亞音速入流邊界:有四條入流特征線，需規(guī)定四個條件;

②亞音速出流邊界:只有一條入流特征線，只需規(guī)定一個條件;

③超音速入流邊界:有五條入流特征線，需規(guī)定全部五個條件，即邊界上取自由流值;

④超音速出流邊界:無入流特征線，不需在邊界上規(guī)定條件，邊界上的值由流場內(nèi)向邊界上外插確定。

4)物面邊界處理

在鄰近物面的有限面積元上，考慮到物面不可穿透，從而有:

5)基本耗散模型

前述的方法中，相鄰單元交接面處的數(shù)值通量采取平均計算。對均勻網(wǎng)格是二級精度的中心差分格式。這種格式有奇偶不關(guān)連的缺陷，需附加四級耗散項防止高頻誤差不阻尼。但這有引起激波附近解的震蕩，又需在此處附加二級耗散而關(guān)閉四級耗散。這就是我們要用的Jameson的二級和四級融合耗散。

6)旋翼壓力盤模型

如圖3，旋翼槳葉的掃掠軌跡可以劃分成槳盤的網(wǎng)格，旋翼槳葉B點處的速度可以確定。

圖3 旋翼槳葉

依據(jù)垂直和平行于槳葉軌跡的速度分量，可以確定誘導(dǎo)攻角αi:

B點的有效攻角αe:

這里，αg是由槳葉幾何形狀和槳葉變距確定的幾何攻角，αg、αe和 αi之間的關(guān)系見圖 4。

圖4 旋翼槳葉剖面的相對速度和受力示意圖

確定了槳葉的相對速度，就能計算作用在旋翼槳葉上的升力和阻力，得到旋翼拉力和扭矩分量:

壓力跳躍邊界條件變?yōu)?

扭矩邊界條件是模擬切向速度跳躍，切向速度跳躍的量值跟誘導(dǎo)速度vi及扭矩拉力比相關(guān)。誘導(dǎo)速度確定后，切向速度跳躍的量值為:

3 計算狀態(tài)及結(jié)果分析

圖5 縱列式雙旋翼直升機機身表面網(wǎng)格

本文對縱列式雙旋翼直升機機身繞流進(jìn)行了數(shù)值計算，通過對雙旋翼不同交疊度δ和直升機不同前飛速度V的計算(交疊度是兩個旋翼相互重疊的最大長度與旋翼半徑的比值)，比較分析機身氣動特性變化。計算包括了三種交疊度:δ=25.5%，δ=26.7%，δ =39%;兩種前飛速度:V=82.8km/h(μ=0.1)，V=245km/h(μ =0.3)。

圖7 截面壓力系數(shù)分布(X=14.5m，δ=39%)

計算表明，飛行速度對機身表面的壓力系數(shù)分布影響較大，而旋翼的交疊度變化對機身表面的壓力系數(shù)分布影響不大。圖6給出了機身表面的壓力系數(shù)分布。

計算表明前后雙旋翼對機身產(chǎn)生的干擾影響，在機身的兩側(cè)產(chǎn)生了不對稱的流場，表現(xiàn)在機身兩側(cè)產(chǎn)生了不對稱的壓力分布。在后機身截面(X=14.5m)處的壓力分布出現(xiàn)了較大的波動，主要是此處受到了動力艙干擾的影響(見圖7);

從機身對稱截面(Y=0)的壓力系數(shù)分布(圖8)可以看到，機頭處為正的壓力系數(shù)，機身中段的壓力很平穩(wěn)，壓力數(shù)值接近0;前后旋翼塔附近的壓力波動較大，特別是動力艙附近的壓力波動更大。

圖9給出了典型截面的流場圖，在X=15.5m截面，由于受到旋翼及動力艙的影響，其流場分布比較復(fù)雜，在動力艙的后部，流場是紊亂的。

圖8 截面壓力系數(shù)分布(Y=0m，δ=39%)

圖9 截面流場圖

4 結(jié)論

通過對縱列式雙旋翼直升機氣動力和流場的數(shù)值模擬，得到如下結(jié)論:

1)飛行速度對機身表面的壓力系數(shù)分布影響較大，表明前飛速度的變化會影響旋翼/旋翼/機身的氣動干擾程度;

2)前后旋翼的交疊度變化對機身表面的壓力系數(shù)分布影響不大;

3)動力艙周圍的渦流場是紊亂的，動力艙周圍機身表面的壓力系數(shù)分布也驗證了空間流場的紊亂特性。

［1］Strawn R C，Barth T J.A finite－volume Euler solver for computing rotary－wing aerodynamics on unstructured meshes［J］.Journal of American Helicopter Society，1993，38(2).

［2］Fejtek I，Roberts L.Navier － Stokes computation of wing/rotor interaction for a tilt rotor in hover［J］.Journal of American Helicopter Society，1992，30(11).

［3］Rajagopalan R G，Mathur S R.Three dimensional analysis of a rotor in forward flight［J］.Journal of American Helicopter Society，1993，38(3).

［4］Zori L A J，Rajagopalan，R G.Navier－Stokes calculation of rotor － airframe interaction in forward flight［J］.Journal of American Helicopter Society，1995，40(2).

［5］Chaffin M S，Berry J.Navier－Stokes Simulation of a Rotor Using a Distributed Pressure Disk Methed［C］.51st American Helicopter Society，May 1995.