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

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

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

0 引言

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

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

圖1 計算流程圖

1 網格生成

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

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

圖2 縱列式直升機的最終表面網格

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

1)數學模型

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

其中xj是n維實數空間Rn的坐標分量，1≤j≤n，解域 Ω?Rn×［0，t)，n=3。W是守恒變量，Fj為流通量矢量，其表達式分別為:

理想氣體的狀態方程為:

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

其中α是來流攻角，M∞是來流馬赫數。

2)有限體積法

應用Jameson的有限體積法及四步Runge－Kutta時間推進，此外，還引入了當地時間步長，焓阻擬等技術，提高收斂速度，減少了計算時間。在遠場邊界處理上采用無反射邊界條件。

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

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

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

從該網格中流出的凈通量為:

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

其中W左、W右分別為面的左右兩側網格中心處的守恒量。

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

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

3)遠場邊界處理

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

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

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

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

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

4)物面邊界處理

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

5)基本耗散模型

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

6)旋翼壓力盤模型

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

圖3 旋翼槳葉

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

B點的有效攻角αe:

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

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

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

壓力跳躍邊界條件變為:

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

3 計算狀態及結果分析

圖5 縱列式雙旋翼直升機機身表面網格

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

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

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

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

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

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

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

圖9 截面流場圖

4 結論

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

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

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

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

［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.