渦槳飛機滑流影響的非定常數值模擬驗證

馬 率,王建濤,邱 名,劉 剛

(中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000)

0 引 言

在螺旋槳飛機的氣動布局設計和性能分析中,存在著滑流和全機的氣動干擾,預測螺旋槳滑流對飛機氣動特性的影響一直是螺旋槳飛機設計中亟待解決又難度很大的問題[1]。由于螺旋槳的滑流是一種非管道控制的流動,不僅軸向速度較遠方的來流速度大,而且還有一定的回轉運動,流線呈螺旋狀[2],因此流動特性非常復雜。國內外學者長期以來不斷嘗試著用各種理論和方法來研究螺旋槳滑流的機理,目前螺旋槳氣動特性計算的方法主要有:動量理論、葉素理論、渦流理論、升力線理論、升力面理論以及各種CFD方法[3]。20世紀70年代以來,Strash[4]等將動量理論與葉素理論結合構造了等效盤模型用于螺旋槳/短倉/機翼組合部件的數值模擬;Sullivan[5]和陳澤民、潘杰元[6]提出了改進的升力線方法,并在此基礎上發展了螺旋槳的升力面方法;西工大的左歲寒等[7]采用基于動量理論的等效盤模型與CFD相結合的方法對螺旋槳/機翼的三維流場進行了求解;南航的李博等[8]發展了考慮螺旋槳滑流旋轉速度的等效盤模型,并與CFD相結合對某四發渦槳飛機的氣動性能進行了計算。以上的這些方法基于不同理論,對真實螺旋槳三維流場都有一定程度的簡化。隨著數值方法和計算機技術的發展,使得對螺旋槳三維流場直接進行Euler/N-S方程求解變得可能。近年來,國內許和勇[9]、楊小川[10]等用非定常Euler方程數值模擬了螺旋槳飛機的滑流影響并進行了分析。由于RANS方程考慮了氣體的黏性,可以更加深入地研究螺旋槳大側流、槳葉失速、滑流掃過機身部件的干擾等復雜流動現象的流場細節,采用URANS方程對單個螺旋槳的氣動特性開展數值模擬研究已有不少成果[11-12],但對螺旋槳/全機的URANS方程數值模擬及驗證的工作還很少。

本文為了研究螺旋槳滑流對某渦槳飛機氣動特性的影響,運用自主開發的大型三維可壓縮RANS方程解算器——PMBNS3D軟件,采用動態重疊網格技術實現槳葉的相對運動,計算采用雙時間步方法,SA一方程湍流模型和k-ωSST兩方程湍流模型,在內循環迭代當中,通過隱式求解、多重網格方法以及結合并行化算法加速收斂過程。通過以上方法對某螺旋槳飛機起飛構型下的低速特性流動進行了研究,并與試驗結果進行了對比,在確保計算結果總體可信的基礎上總結了螺旋槳滑流動對飛機全機及部件氣動特性和操控性能的影響規律,然后通過對螺旋槳滑流影響邊界區域定義及分析,對滑流干擾現象以及背后的物理流動機理進行了研究和探索。通過這些工作的開展,表明本文所發展的螺旋槳/機身流場的URANS方程數值方法可以為我國渦槳飛機操縱性及穩定性設計及分析提供技術支撐。

1 數值方法

1.1 控制方程

在慣性坐標下,針對任意運動控制體的的可壓縮N-S方程的積分形式為:

其中,Ω表示控制體的體積;S表示控制體邊界面的面積;VΩ為邊界面上的網格運動速度;Q為狀態變量:Q＝{ρ,ρu,ρv,ρw,e}T,這里ρ、u、v、w、e分別表示流體的密度、速度矢量在三個坐標方向(x、y、z)上的分量和單位體積的總能量。壓力p和內能ε以及和總能e之間滿足如下狀態方程:

無黏項的空間離散采用Roe平均迎風通量差分分裂格式,黏性項采用中心差分格式離散,計算分別采用了SA一方程湍流模型和k-ωSST二方程湍流模型[13-14]。非定常計算采用雙時間步方法[15],除了真實時間t之外,引入偽時間τ,好處是使原來只適用于定常計算的隱式時間迭代、局部時間步長、多重網格等方法都可以運用到非定常計算的偽時間迭代當中[16-17]。

1.2 多重網格方法

限制URANS方程應用的一個重要因素是對計算資源的要求過高。對雙時間步方法來說,提高子迭代過程的收斂速率是衡量方法的一個重要指標,而在加速方法中多重網格法被證明是最有效的手段之一[18-19]。多重網格方法的原理就是利用高頻誤差比低頻誤差更容易被迭代格式消除的特點,通過引入系列逐步粗化網格的辦法,將低頻誤差轉化為高頻誤差,消除各種誤差分量,再用套迭代技術將粗網格上的解返回到各級細網格上。多重網格法從理論上被證明至少對于線性橢圓型問題是一種最優化的數值方法,其計算工作量僅僅與網格節點數的一次方成正比,并且收斂速度與網格的尺度大小無關,從而特別適合于應用在超大型工程數值計算問題中。本文采用了多重網格法中非線性方程的全近似格式(FAS),套迭代技術為W循環。

1.3 動態重疊網格技術

在采用時間精確的數值方法求解中,網格運動是一個比較棘手的問題,較為成熟的方法是結構網格框架下的動態重疊網格方法。重疊網格[20-21]對相對運動的部件生成不同的網格,分別隨各自部件運動,相互之間構成動態的重疊關系。由于部件運動的關系,有必要在網格生成上對運動路徑上的網格進行互相適應,照顧好這個區域的網格密度以及均勻程度,以滿足重疊插值的要求;同時要求不同部件的網格在接合區域相互匹配,對它們之間流場的光滑過渡也是十分必要的[22-23]。

1.4 并行化計算

并行計算采用基于消息傳遞(MPI)的并行化算法,用區域分割方法將計算域映射到多個CPU上分別進行計算[24-25]。由于計算網格采用多塊結構網格方法,在網格塊數小于節點數或者現有網格分塊不足以使負載平衡的情況下,對最大塊進行剖分,剖分需要滿足兩個條件,一是沿最大維數方向進行剖分,這樣保證信息交換量最少;第二要求剖分后網格滿足多重網格計算條件。如果最大維數方向滿足不了該條件的話就沿次大維數方向進行剖分,依次進行。按這種方法劃分的網格可以在各計算節點上進行多重網格計算,提高整體計算效率[26]。

在發展了動態重疊網格的快速自動查詢方法基礎之上,實現了重疊網格挖洞、插值關系建立的完全并行化。其中找插值關系的具體過程可以描述如下:對任意插值點而言,首先在本節點查找貢獻單元,在沒有找到的情況下將該插值點信息發送到其余各個節點上,由這些節點同時查找該插值點的貢獻單元,如果找到,將該貢獻單元信息返回,并且建立兩個節點之間的消息傳遞關系。這種做法的最大好處就是去掉了所謂的根節點,即任何節點的重疊網格計算不須收集網格全局信息,降低了單節點的內存要求,同時避免了單個節點計算重疊插值關系,其余節點等待的局面,提高了計算效率。

2 螺旋槳單槳計算驗證

研究模型為單獨渦槳發動機吊艙帶六片槳葉模型,見圖1。螺旋槳直徑4 m,槳葉角為36°,轉速為7550 r/min,來流風速25 m/s,迎角為0°,模型縮比為1∶12。

計算網格采用重疊網格方法。分別對運動部件(槳葉及輪轂)和靜止部件(發動機吊艙)生成各自的網格。圖2中給出了槳葉網格和發動機吊艙背景網格的重疊關系,槳葉網格隨槳葉一起作剛體運動,動態地與發動機吊艙背景網格構成重疊關系。槳葉網格的拓撲結構為O型。為了模擬黏性附面層的需要,第一層網格距離物面的距離約為1×10-5m,網格單元總數約330萬。

圖1 帶六片槳葉的渦槳發動機外形Fig.1 The propeller engine with six blades

圖2 槳葉網格和發動機吊艙背景網格的重疊關系Fig.2 Overlap relationship between blades grids and engine background grids

流場計算從初場開始進行非定常計算,計算采用在慣性坐標系下求解URANS方程方法,在時間上計算采用雙時間步方法,每個真實時間步槳葉在周向運動3°,即每個旋轉周期包含真實時間步數為120步,每個真實時間步內的子迭代步數為50步,采用隱式LU-SGS方法進行迭代,空間格式采用Roe通量差分分裂方法,湍流模型分別采用了SA一方程湍流模型和k-ωSST二方程湍流模型。

并行計算在工作站集群上完成,計算結點為12核至強芯片(Intel(R)Xeon(R)CPU X5670,主頻2.93 GHz,內存48 GB),該單槳計算驗證采用了8核的并行計算,加速比為7,并行效率約為87%。

圖3給出了總的拉力系數CT和功率系數CP的收斂歷程,橫坐標表示槳葉周向角的變化,從圖中看到,拉力系數在槳葉運動1.5周以后基本達到收斂狀態。由于螺旋槳包含6片槳葉,每片槳葉的運動規律相同,所以流場變化的周期實際上是1/6周。

圖4給出的是螺旋槳效率計算結果與試驗結果的比較,可以看出采用SST湍流模型的計算結果與試驗值吻合的很好,誤差最大不超過2%,因此在以后的計算中均采用了SST湍流模型。渦槳發動機吊艙帶六片槳葉模型的計算結果證實了動態重疊網格及非定常方法解算器的可靠性,可用于螺旋槳飛機的初步氣動特性計算分析。

圖3 總的拉力系數和功率系數的收斂歷程Fig.3 Convergence course of propeller thrust coefficient and power coefficient

圖4 螺旋槳效率計算結果與試驗結果的比較Fig.4 Comparison of propeller efficiency between calculation and experimental data

3 某渦槳飛機起飛構型有無滑流氣動特性影響的計算驗證與分析

某渦槳飛機起飛構型無動力及帶動力外形如圖5所示,模型基本參數如下:

模型縮比: 1∶12

機翼面積: 0.4626 m2

機翼平均氣動弦長: 0.20125 m

試驗風速: 30 m/s

螺旋槳轉速: 4977 r/min(拉力TC＝0.15)

襟翼偏度: 35°

螺旋槳旋轉方向為雙發同向順時針旋轉(飛行員視角),對該起飛構型下的低速特性流動進行了計算分析。

3.1 計算結果的可靠性驗證

計算網格分為無動力全機和帶動力全機兩套網格,見圖5,無動力全機網格采用對接網格方法,網格單元總數約為2200萬,無動力全機單個定常計算狀態采用8核并行計算,迭代計算1萬步取收斂結果,耗時約為24 h。帶動力全機計算網格采用重疊網格方法,分別對運動部件(槳葉及輪轂)和靜止部件(全機其它部分)生成各自的網格,槳葉網格隨槳葉一起運動,動態地與全機背景網格構成重疊關系,為了模擬黏性附面層的需要,第一層網格距離物面的距離約為3×10-5m,y＋＝O(1),湍流模型采用了k-ωSST二方程湍流模型。全機網格規模約3050萬,帶動力全機單個非定常計算狀態采用64核并行計算,加速比為49,并行效率約為75%,迭代計算2520個真實時間步數(螺旋槳旋轉21轉)取流場變化的周期值,耗時約135個小時。

圖5 全機無動力及帶動力起飛構型物面網格Fig.5 Grid of takeoff airplane surface without and with propellers

圖6給出了無動力及帶動力不同情況下氣動特性隨迎角變化的計算值與試驗值的比較。需要特別說明的是為了和試驗結果對比,圖中的升力和阻力特性曲線乃是扣除了槳葉的直接拉力,而帶動力的俯仰力矩曲線則是考慮了槳葉的俯仰力矩。從試驗結果和計算結果的對比情況來看,升力系數在線性段幾乎重合,計算結果與試驗值在線性段的總體誤差不超過2%,在失速點附近因為機翼上出現的分離區導致數值模擬的難度較大,造成了計算與實驗的一些偏差;阻力系數計算結果相對試驗值整體偏大,但誤差基本控制在8%以內,偏大原因應該來自于數值計算解的耗散行為以及采用全湍流模型帶來的摩阻的增加;俯仰力矩曲線在線性段帶動力情況的計算值與試驗值吻合較好,而在勺型區域計算值與試驗值差別較大,只有帶動力情況的計算值與試驗值出現拐點的迎角一致;在滑流影響作用下,全機升力、阻力系數有明顯提升,且隨著迎角的變大,升阻系數的增量值也開始放大,而俯仰力矩系數斜率的絕對值減小,說明帶螺旋槳動力的全機俯仰安定性對比無動力有所降低。

圖6 全機無動力與帶動力的氣動特性變化Fig.6 Change of takeoff airplane’s aerodynamics characteristics(without and with propellers)

圖7 無動力起飛的全機各部件的氣動力特性(α＝8°)Fig.7 Parts contribution in takeoff airplane’s aerodynamics characteristics(without propellers,α＝8°)

3.2 滑流對全機及部件的氣動力影響分析

與巡航狀態不同的是,即使在小迎角狀態下,滑流對起飛狀態飛機的升力及阻力亦有明顯的放大作用,可見襟翼的打開導致流場阻塞效應進一步增大,這個阻塞作用通過壓力場施加給螺旋槳一個反作用力,從而使螺旋槳拉力及滑流效應進一步增強。由于低速起飛時迎角8°為典型飛行狀態,下面對該狀態進行細節上的分析,圖7是無動力巡航的全機各部件對氣動力特性的貢獻圖,圖8給出了帶動力巡航的全機各部件對氣動力特性的貢獻圖。在升力的貢獻圖對比中可以看到,有無動力全機的升力都主要來自于機翼的貢獻,它占到總升力的70%左右,其次機身貢獻占總升力的11%左右,另外襟翼貢獻能占到總升力的9%左右;受滑流影響,帶動力影響的機翼升力比無動力偏大了約10%,機身升力比無動力偏大了約19%,襟翼升力比無動力偏大約60%;由于螺旋槳同向旋轉的原因,左側(即關鍵發動機側)的機翼升力比右側機翼大2%;在滑流影響下左右平尾的升力發生明顯不對稱變化,值得注意的是機身和立尾受到的是右向側力,其中立尾的右向側力非常明顯,而左機翼和右機翼受到的側力值亦不小,但它們各自反向而近似抵消;俯仰力矩的貢獻圖中可以發現,機翼、襟翼、平尾和槳葉的低頭力矩有利于縱向的靜穩定性,特別是右平尾的低頭力矩相比無動力時明顯減弱,而此時機身有很大的抬頭力矩,輪轂也有一定的抬頭力矩,它們對全機的縱向安定性帶來了不利的影響。

圖8 帶動力起飛的全機各部件的氣動力特性(α＝8°)Fig.8 Parts contribution in takeoff airplane’s aerodynamics characteristics(with propellers,α＝8°)

圖9 有動力及無動力機翼弦向站位的壓力系數比(Ma＝0.088,α＝8°,T C＝0.15)Fig.9 C p distributions for different stations located on wing chordwise(without and with propellers,Ma＝0.088,α＝8°,T C＝0.15)

圖9給出了低速起飛迎角8°時無動力及帶動力機翼的弦向截面站位壓力系數對比圖,可以看到經過機翼的導流作用后,滑流對襟翼的洗流不像機翼那么明顯,但由于滑流區襟翼上的動壓還是得到明顯增加,同時高速滑流也起到一個吹除在大偏度襟翼上堆積的附面層作用,延遲了襟翼上氣流的分離,這兩個因素明顯提高了襟翼效率。另外左側(即關鍵發動機側)的襟翼升力比右側偏大一些,這是因為左側襟翼的當地來流迎角比右側的略大一點。

圖10(a~c)給出了無動力起飛外形隨迎角變化的物面極限流線圖。圖中顯示在小迎角2°條件下起飛狀態的機翼沒有分離;迎角到8°時靠近襟翼外側的機翼上表面出現了小范圍的展向分離流動,它是由襟翼的翼尖渦誘導而引起的;迎角增大到14°后,機翼上表面的出現了大范圍的分離,靠近襟翼外側的機翼展向分離已經接近到了吊艙的外側機翼位置,而吊艙內側的機翼翼根大面積分離是由吊艙的干擾誘發的,機翼明顯產生了失速現象。圖10(d~f)是帶動力起飛外形不同迎角的物面極限流線圖,與無動力外形相比,迎角14°時,左機翼上靠近襟翼外側的機翼展向分離面積擴大,這是由于滑流旋轉與機翼翼尖渦方向同一側而更容易誘發展向分離,左邊發動機內側的機翼翼根分離被抑制是受到滑流對左機翼當地迎角的上洗影響,同理右機翼上靠近襟翼外側的機翼展向分離被抑制是由于滑流旋轉與右機翼翼尖渦方向反向,關鍵發動機內側的機翼翼根分離被擴大是因為滑流對右機翼當地迎的下洗影響。

圖10 無動力和帶動力全機物面極限流線隨迎角變化情況(Ma＝0.088)Fig.10 Comparison of the pressure coefficients and streamlines on airplane’s surface(without and with propellers,Ma＝0.088)

圖11 帶動力起飛狀態空間渦量等值面圖Fig.11 Iso-surfaces of the Q-criterion colored by C p

3.3 滑流的空間流場特性

圖11給出的是帶動力在低速起飛狀態下的空間渦量等值面圖。這時可以看到螺旋槳后拖出的渦量,經過機翼及襟翼進一步被切割打散后拖出的尾渦圈數明顯減少,相互融合成帶狀分布,并在向下游發展的過程中逐漸淡化變得不再明顯。

為了研究滑流的周期時均影響效果,這里采用有無滑流不同情況下對比當地動壓增量來定義滑流的加速效應邊界,對比當地氣流角增量來定義滑流的洗流效應邊界,期望通過對滑流起主要影響作用的參數分解來判定滑流的影響區域。為了達到以上目的,需要將螺旋槳停止旋轉時定常狀態的流場和螺旋槳正常旋轉的非定常滑流的周期時均流場進行對比,得到兩者間的加速作用和洗流的增量變化。圖12~圖14分別給出了在滑流影響作用下,低速起飛狀態空間流場不同截面的動壓增量Δq、當地迎角增量Δα和當地側滑角增量Δβ的變化情況。可以發現空間中除槳葉卷出的滑流流管以外的地方動壓增量Δq近似為0,流管中機翼上下翼面附近處的動壓增量能達到0.006左右,相比流場遠處來流的動壓q增加了近1.5倍以上(遠方來流的動壓q等于0.00387),相比巡航狀態,起飛狀態下滑流流管的范圍經過襟翼后變得更大,且動壓增量Δq的值更高,說明起飛狀態下襟翼打開后相比巡航狀態對滑流的阻塞效應進一步增強,會增大螺旋槳拉力,使滑流的加速效應和洗流效應更明顯。圖13中能明確看到在襟翼打開后空間中下洗范圍比巡航狀態擴大且被進一步增強,右側平尾的前方有約－3°左右的下洗增量,左側平尾的前方卻有約1°左右的上洗增量,左右不對稱效應被放大,導致飛機更低的縱向靜安定度。圖14中給出了有無滑流情況下當地氣流側滑角增量Δβ的空間剖面圖,能夠看到滑流對空間流場側洗的影響,在立尾翼根的前緣處來流相比巡航狀態有更明顯的右偏斜,而機身后段右側的來流相比巡航狀態不再存在明顯的右偏斜,這樣使得機身受到的側力相比巡航狀態發生了方向性的變化。

圖12 帶動力起飛狀態動壓增量變化的空間剖面圖(α＝8°)Fig.12 Incremental of dynamic pressure for different sections in flow field(with propellers,α＝8°)

圖13 帶動力起飛狀態當地迎角增量變化的空間剖面圖(α＝8°)Fig.13 Incremental of local angle of attack for different sections in flow field(with propellers,α＝8°)

圖14 帶動力起飛狀態當地側滑角增量變化的空間剖面圖(α＝8°)Fig.14 Incremental of local angle of sideslip for different sections in flow field(with propellers,α＝8°)

4 結 論

本章對有無滑流影響作用的某雙發同向旋轉渦槳飛機分別進行了巡航和起飛狀態的數值模擬,與試驗結果進行了對比并對滑流的影響作用開展了詳細的分析,得到了以下一些結論:

(1)螺旋槳后拖出的螺旋狀滑流改變了機翼的流態,其影響大致可以分為兩個方面:首先螺旋槳的轉動給氣流注入了能量,使得滑流區域機翼、襟翼上的動壓和靜壓都明顯增加,產生了一個升力增量和誘導阻力增量;其次改變了機翼上的局部迎角,使滑流區內螺旋槳上行槳葉側上洗增強,增加了局部有效迎角,從而在機翼上誘導出一個正升力增量,同理,滑流區內螺旋槳下行槳葉側下洗增強,減小了局部有效迎角,在機翼上誘導出一個負升力增量。

(2)當飛機襟翼打開后,滑流區內襟翼上的動壓得到明顯增強,同時高速滑流也起到一個吹除在大偏度襟翼上堆積的附面層作用,延遲了襟翼上氣流的分離,這兩個因素明顯提高了襟翼效率。

(3)采用有無滑流不同情況下對比當地動壓增量和當地氣流角增量來定義滑流的加速效應和洗流效應邊界,采用該方法能較好地捕捉和解釋由于滑流對飛機尾翼干擾而使得飛機方向安定性呈現的很大非線性現象,增進了滑流對掃過飛機部件的氣動干擾作用機理的認識。