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

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

0 引言

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

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

圖1 計(jì)算流程圖

1 網(wǎng)格生成

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

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

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

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

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

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

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

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

其中p是壓強(qiáng)，γ為比熱比，對(duì)于空氣，γ=1.4。

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

2)有限體積法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4)物面邊界處理

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

5)基本耗散模型

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

6)旋翼壓力盤模型

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

圖3 旋翼槳葉

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

B點(diǎn)的有效攻角αe:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 截面流場(chǎng)圖

4 結(jié)論

通過(guò)對(duì)縱列式雙旋翼直升機(jī)氣動(dòng)力和流場(chǎng)的數(shù)值模擬，得到如下結(jié)論:

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

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

3)動(dòng)力艙周圍的渦流場(chǎng)是紊亂的，動(dòng)力艙周圍機(jī)身表面的壓力系數(shù)分布也驗(yàn)證了空間流場(chǎng)的紊亂特性。

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