張馨元

摘 要：采用Roe通量微分分裂格式離散Navier-Stokes（N-S）方程組的對流項(xiàng)，用統(tǒng)一時(shí)間步長的偽時(shí)間子迭代隱式LU-SGS方法進(jìn)行時(shí)間推進(jìn)，利用變形網(wǎng)格方法求解任意拉格朗日-歐拉（ALE）坐標(biāo)下的N-S方程組，對NACA-0012翼型小幅振蕩和動(dòng)態(tài)失速以及三段翼型振蕩襟翼等非定常流動(dòng)進(jìn)行了數(shù)值分析；通過與可獲得的實(shí)驗(yàn)結(jié)果對比，吻合良好。與此同時(shí)，還對比分析了代數(shù)B-L、兩方程SST及k-g等三種湍流模式的表現(xiàn)，總體上SST和k-g模式結(jié)果優(yōu)于B-L模式；在動(dòng)態(tài)失速流動(dòng)中，SST模式較k-g模式更靠近測量值；采用SST模式對三段翼型振蕩襟翼流場進(jìn)行詳細(xì)研究，分析流動(dòng)特點(diǎn)及其與升力特性之間的關(guān)系。

關(guān)鍵詞：動(dòng)網(wǎng)格 深失速 振蕩翼型/襟翼

中圖分類號：V21 文獻(xiàn)標(biāo)識碼：A 文章編號：1674-098X（2017）10（b）-0014-06

Abstract：The Roes Flux Difference Splitting is employed as the spatial discretization scheme for the convective terms of the Navier-Stokes equations. The fully implicit lower-upper symmetric- Gauss-Seidel with pseudo time sub-iteration is applied as the temporal marching methods for both N-S equations and turbulence model equations. Flows around pitching airfoils or flap are investigated by solving N-S equations based on a moving grid system which described by Arbitrary Largangian-Eulerian method. Three turbulence models， such as algebraic Baldwin- Lomax， two-equation shear stress transport （SST） and k-g models， are employed to study the flows around pitching and dynamic stall NACA-0012， and the pitching flap of the three-element airfoil. The SST and k-g models perform better than B-L model. In the dynamic stall case， the SST model is better than the k-g model. Flowfields and flow characteristics of the pitching flap are investigated in detail by SST with respect to lift coefficients.

Key Words：Moving grid； Dynamic stall； Pitching airfoil/flap

隨著計(jì)算機(jī)技術(shù)和計(jì)算方法的飛速發(fā)展，計(jì)算流體力學(xué)（CFD）在飛機(jī)和導(dǎo)彈設(shè)計(jì)中得到了越來越廣泛的應(yīng)用。具有運(yùn)動(dòng)邊界的非定常復(fù)雜流場，如機(jī)翼顫振、操縱面嗡鳴、繞流片擺動(dòng)、襟翼偏轉(zhuǎn)、外掛物投放等，一直是計(jì)算流體力學(xué)研究的熱點(diǎn)和難點(diǎn)問題。若上述流動(dòng)發(fā)生在跨音速階段，則還需考慮激波/邊界層的相互干擾；若飛機(jī)經(jīng)歷動(dòng)態(tài)失速，需要準(zhǔn)確模擬大范圍的流動(dòng)分離；當(dāng)飛機(jī)起飛和著陸時(shí)，襟翼偏轉(zhuǎn)或振動(dòng)對飛機(jī)的升力、阻力以及俯仰力矩特性產(chǎn)生重要影響，若嚴(yán)重時(shí)可影響飛行安全和誤導(dǎo)駕駛員操縱，有必要深入研究俯仰翼型/襟翼流動(dòng)。

當(dāng)數(shù)值求解Navier-Stokes方程組時(shí)，需要引入湍流模式使方程組進(jìn)行封閉。非定常動(dòng)邊界問題時(shí)，尤其是激波/邊界層干擾及大范圍流動(dòng)分離問題，湍流模式的作用更加明顯，其表現(xiàn)不容忽視。在飛機(jī)氣動(dòng)性能分析中，代數(shù)B-L、一方程S-A模式及兩方程SST模式的應(yīng)用極為廣泛；然而，在分析非定常動(dòng)邊界繞流時(shí)，很多研究仍只考慮B-L模式。

本文發(fā)展了自行研究開發(fā)的多塊結(jié)構(gòu)網(wǎng)格程序NSAWET[1]（N-S Analysis based on Window- Embedment Technology），求解隨網(wǎng)格點(diǎn)一起移動(dòng)的運(yùn)動(dòng)坐標(biāo)系即Arbitrary Lagrangian-Eulerian （ALE）的N-S方程組作為控制方程，利用變形網(wǎng)格方法對振蕩翼型/襟翼進(jìn)行數(shù)值研究；采用3種湍流模式對NACA0012小幅振蕩和動(dòng)態(tài)失速，以及三段翼型振蕩襟翼等非定常流動(dòng)進(jìn)行數(shù)值研究。

1 數(shù)值方法

基于動(dòng)網(wǎng)格技術(shù)研究俯仰翼型/襟翼的非定常流動(dòng)特點(diǎn)，因此變形或運(yùn)動(dòng)邊界對控制方程組、湍流模式、時(shí)間推進(jìn)方法以及邊界條件的影響被作為本文工作的重點(diǎn)，現(xiàn)分述如下。

1.1 控制方程組

本文所采用的控制方程組為隨網(wǎng)格點(diǎn)一起移動(dòng)的運(yùn)動(dòng)坐標(biāo)系即ALE的N-S方程組[2，3]，與Euler描述下的N-S方程組相比，ALE的N-S方程組對流項(xiàng)中增加了由網(wǎng)格運(yùn)動(dòng)而引起的通量。

對流項(xiàng)采用通量微分分裂類Roe格式離散，粘性項(xiàng)采用二階中心差分格式離散。

1.2 湍流模式

主要采用B-L、兩方程k-ωSST以及k-g模式使N-S方程組封閉。代數(shù)B-L模式形式簡單，勿需求解額外方程，SST及k-g模式需要求解關(guān)于k、ω或g的湍流模式方程組。endprint

2 計(jì)算結(jié)果及分析

2.1 NACA0012翼型小幅振蕩

振蕩NACA0012翼型具有較為詳細(xì)的實(shí)驗(yàn)數(shù)據(jù)，被許多研究者所引用，也被作為本文的驗(yàn)證算例。翼型繞1/4弦點(diǎn)的振蕩方式為：α（t）=αm+α0sin（2kt），其中αm為平均攻角，α0最大振幅攻角，k=為無量綱的振蕩頻率，為振蕩頻率，c為弦長，為來流速度。采用C型網(wǎng)格，翼型弦長為1，周向和法向網(wǎng)格數(shù)目分別為501和101，法向第一層網(wǎng)格距離物面為1.0×10-5。

2.1.1 亞音速情形

根據(jù)文獻(xiàn)[6]，振蕩翼型的流動(dòng)參數(shù)為：Ma=0.60，Re=4.8×106，αm=4.86°，α0=2.44°，k=0.0810。無量綱時(shí)間步長0.05，一個(gè)振蕩周期約800步，內(nèi)迭代20次，湍流模式為SST。

圖1是升力系數(shù)和俯仰力矩系數(shù)隨攻角變化的遲滯曲線及其與實(shí)驗(yàn)數(shù)據(jù)以及文獻(xiàn)[6]的結(jié)果對比，圖2是幾個(gè)典型時(shí)刻對應(yīng)的壓力系數(shù)Cp曲線。從圖1與圖2可以看出，本文數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)吻合非常好，且升力系數(shù)比文獻(xiàn)[6]中利用k-ε模式的結(jié)果更接近實(shí)驗(yàn)值。通過Cp分布可以看出，振蕩過程中，上揚(yáng)和下俯時(shí)激波的表現(xiàn)差異非常大，即便下俯時(shí)的瞬間攻角較?。é?3.49°），激波仍很強(qiáng)，可是翼型上揚(yáng)時(shí)瞬間攻角較大（α=4.28°）激波反而較弱，從客觀上反映了流動(dòng)的遲滯效應(yīng)。當(dāng)流動(dòng)中出現(xiàn)激波并逐漸增強(qiáng)后，激波與邊界層的干擾非常嚴(yán)重，導(dǎo)致了流動(dòng)的強(qiáng)烈非線性，增加了數(shù)值分析的難度。

2.1.2 跨音速算例

該算例[3，5，6]中流動(dòng)參數(shù)為Ma=0.755，Re=5.5×106，αm=0.016°，α0=2.51°，k=0.0814。

圖3給出了B-L、SST及k-g模式計(jì)算的升力和力矩系數(shù)隨攻角變化的遲滯曲線，及其與實(shí)驗(yàn)及文獻(xiàn)[6]參考值的對比?？梢钥闯?，無論是文獻(xiàn)還是本文方法，所得到的升力系數(shù)均比實(shí)驗(yàn)偏低。從圖3還可以看出，SST和k-g模式對升力和力矩系數(shù)的模擬均較好，而且相互之間的差別非常小。B-L模式所得的升力遲滯環(huán)稍小，而力矩遲滯環(huán)稍大，計(jì)算效果較差。

2.2 NACA0012翼型動(dòng)態(tài)失速

動(dòng)態(tài)失速是指運(yùn)動(dòng)的壓力面在超過其臨界迎角時(shí)流場發(fā)生非定常分離和失速的現(xiàn)象。根據(jù)粘性作用的大小，動(dòng)態(tài)失速可分為輕失速和深失速。輕失速除表現(xiàn)出一般靜態(tài)失速的特征外，還有非定常分離的強(qiáng)粘性/無粘相互作用性質(zhì)，邊界層厚度為翼型厚度量級。深失速則表現(xiàn)為全粘性現(xiàn)象，存在高度非線性的壓力脈動(dòng)和翼型表面上有大尺度旋渦運(yùn)動(dòng)，邊界層厚度可達(dá)到弦長量級[7]。

采用兩方程SST和k-g模式分別計(jì)算深失速問題。文獻(xiàn)[7]給出了NACA0012翼型正弦振蕩深失速的實(shí)驗(yàn)結(jié)果，振蕩規(guī)律與2.1節(jié)所述相同，只是參數(shù)有所改變，Ma=0.20，Re=1.0×106，αm=15°，α0=10°，k=0.15。為準(zhǔn)確分析流動(dòng)中的非定常效應(yīng)，更細(xì)致的描述流場，無量綱時(shí)間步長為0.005，因此每一時(shí)間步攻角變化相對較小。

升力和力矩系數(shù)隨攻角變化的遲滯曲線見圖4，在翼型上揚(yáng)過程中，計(jì)算結(jié)果與實(shí)驗(yàn)符合較好；當(dāng)翼型下俯運(yùn)動(dòng)時(shí)，流動(dòng)中由于受到強(qiáng)逆壓梯度的影響，翼型上表面分離劇烈，兩種模式均不能很好的模擬流動(dòng)的分離特點(diǎn)，對其非線性的力和力矩變化值模擬能力較差，然而可喜的是，他們均可不同程度地模擬下俯過程中力和力矩的變化趨勢。從圖中可以看出，SST比k-g模式效果更好，這得益于SST模式的弱非線性渦粘系數(shù)的定義和對逆壓梯度流動(dòng)中渦粘系數(shù)進(jìn)行限制，使其具有剪切應(yīng)力輸運(yùn)的特點(diǎn)。B-L模式不適于計(jì)算此類流動(dòng)，因而沒有采用。

2.3 三段翼型振蕩襟翼

選擇文獻(xiàn)[8]中的三段翼型，其前緣縫翼和后緣襟翼偏角均為30°。采用法向外推的單塊O型網(wǎng)格（如圖5所示），網(wǎng)格點(diǎn)數(shù)為582×131，法向第一層網(wǎng)格到物面的距離為1.0×10-5。

首先，對馬赫數(shù)0.2，攻角8.109°，Re=9.0×106的定常繞流采用SST模式進(jìn)行數(shù)值分析，三段翼型表面的壓強(qiáng)系數(shù)Cp對比如圖5右所示，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)[8]吻合很好。

其次，在上述定常流動(dòng)的基礎(chǔ)上，數(shù)值分析振蕩襟翼的流動(dòng)特點(diǎn)。給定襟翼繞其10%弦點(diǎn)振蕩的運(yùn)動(dòng)方式為α（t）=αm+α0sin（2kt），其中k=0.10，αm為襟翼基本偏角（30°），α0為10°。計(jì)算得到的三段翼型的總升力系數(shù)與各部件升力系數(shù)隨襟翼偏角變化遲滯曲線如圖6所示；與之對應(yīng)，圖7給出一個(gè)周期內(nèi)典型時(shí)刻的流線。對照圖6與圖7可以看出，一個(gè)周期內(nèi)流動(dòng)狀態(tài)可以根據(jù)襟翼偏角分為以下幾個(gè)階段：

（1）20.0°到36.3°（圖中8→1→2），總升力系數(shù)Cl與襟翼偏角基本保持線性變化，襟翼的升力系數(shù)呈弱非線性增加，斜率逐步減??；偏角較小時(shí)流動(dòng)保持附著，當(dāng)偏角較大時(shí)，在襟翼后緣逐漸出現(xiàn)小的分離泡；襟翼提供的升力在36.3°到達(dá)最大。

（2）36.3°到38.2°（圖中2→3），總升力繼續(xù)增大，直至到達(dá)最大升力點(diǎn)，但上升斜率下降很快；主翼和前緣縫翼升力繼續(xù)增加，由于襟翼后緣出現(xiàn)了較明顯的旋渦，襟翼提供的升力已有所下降。

（3）38.2°到40.0°（圖中3→4），總升力隨偏角的增加不升反降；襟翼升力急劇下降，后緣分離區(qū)持續(xù)增大；主翼和前緣縫翼的升力均經(jīng)歷了從升高到低的過程。

（4）40.0°到33.2°（圖中4→5），由于襟翼偏角減小，三段翼型的各段升力均降低，襟翼后緣分離區(qū)達(dá)到最大，襟翼提供的升力達(dá)到最?。坏怯捎诮笠淼母┭鲞\(yùn)動(dòng)所造成的擾動(dòng)傳播尚未返回，因而主翼和前緣縫翼的升力系數(shù)還處于下探過程中。

（5）33.2°到20.0°（圖中5→8），襟翼升力隨偏角減小反而升高，然而受到襟翼俯仰運(yùn)動(dòng)的影響，總升力減?。唤笠砗缶壏蛛x區(qū)逐漸減小，直至重新達(dá)到附著狀態(tài)。endprint

由上述分析可知，各部件升力系數(shù)變化均具有較強(qiáng)的非線性特性，其最大值不發(fā)生在最大偏角，襟翼升力最小值也不在最小偏角。雖然主翼對在升力的貢獻(xiàn)占了主要部分，但非線性特性主要由襟翼運(yùn)動(dòng)引起分離區(qū)變化所致。前緣縫翼升力則一直保持較弱的非線性特征，由于和襟翼相隔較遠(yuǎn)，對襟翼的俯仰變化敏感程度低。

3 結(jié)語

基于變形結(jié)構(gòu)網(wǎng)格，采用三種不同表現(xiàn)的湍流模式對振蕩翼型/襟翼流動(dòng)進(jìn)行了數(shù)值分析。通過與標(biāo)準(zhǔn)算例的驗(yàn)證，表明本文方法的正確性和計(jì)算的可靠性；通過計(jì)算結(jié)果發(fā)現(xiàn)，兩方程SST和k-g模式，尤其SST模式，能捕捉流動(dòng)中的非定常效應(yīng)和模擬較強(qiáng)逆壓梯度，在分析運(yùn)動(dòng)邊界流動(dòng)時(shí)獲得了與實(shí)驗(yàn)吻合一致的效果。通過對單段翼型小幅度、深失速等俯仰振蕩的準(zhǔn)確模擬，并將其擴(kuò)展到三段翼型繞流中的襟翼俯仰運(yùn)動(dòng)，可獲得增升裝置中單獨(dú)翼對整個(gè)機(jī)翼氣動(dòng)特性的認(rèn)識和了解，為設(shè)計(jì)增升裝置打下堅(jiān)實(shí)的基礎(chǔ)。

參考文獻(xiàn)

[1] Chen Haixin， Zhang Yufei， Xiao Zhixiang， Fu Song， NSAWET Code and Its Application [A].5th Asia Workshop on Computational Fluid Dynamics[C].2006.

[2] I.Demirdzic， M. Peric， Space conservation law in finite volume calculations fluid flow [J].International journal for numerical methods in fluids，1988（8）：1037-1050.

[3] C.Michler，H.De Sterck，H.Deconinck，An arbitrary Lagrangian Eulerian formulation for residual distribution schemes on moving grids [J].Computers & Fluids，2003（32）：59-71.

[4] 肖志祥.復(fù)雜流動(dòng)Navier-Stokes方程數(shù)值模擬及湍流模型應(yīng)用研究[D].西安：西北工業(yè)大學(xué)，2003.

[5] A.Jahangirian， M. Hadidoolabi， Unstructured moving grids for implicit calculation of unsteady compressible viscous flows [J].INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS，2005（47）： 1107-1113.

[6] Zhang L. P.，Wang Z. J.，A Block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes[J].Computers & Fluids， 2004（33）：891-916.

[7] 錢煒祺，符松，蔡金獅.翼型動(dòng)態(tài)失速的數(shù)值研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2001，19（4）：427-433.

[8] Daryl L. Bonhaus， W. Kyle Anderson， Dimtri J. Mavriplis， Numerical study to assess sulfur Hexafluoride as a medium for testing multielement airfoils NASA technical[Z]. 1995.endprint