張馨元

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

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

中圖分類號：V21 文獻標識碼：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

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

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

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

1 數(shù)值方法

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

1.1 控制方程組

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

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

1.2 湍流模式

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

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

2.1 NACA0012翼型小幅振蕩

振蕩NACA0012翼型具有較為詳細的實驗數(shù)據(jù)，被許多研究者所引用，也被作為本文的驗證算例。翼型繞1/4弦點的振蕩方式為：α（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ù)文獻[6]，振蕩翼型的流動參數(shù)為：Ma=0.60，Re=4.8×106，αm=4.86°，α0=2.44°，k=0.0810。無量綱時間步長0.05，一個振蕩周期約800步，內(nèi)迭代20次，湍流模式為SST。

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

2.1.2 跨音速算例

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

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

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

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

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

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

2.3 三段翼型振蕩襟翼

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

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

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

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

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

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

（4）40.0°到33.2°（圖中4→5），由于襟翼偏角減小，三段翼型的各段升力均降低，襟翼后緣分離區(qū)達到最大，襟翼提供的升力達到最??；但是由于襟翼的俯仰運動所造成的擾動傳播尚未返回，因而主翼和前緣縫翼的升力系數(shù)還處于下探過程中。

（5）33.2°到20.0°（圖中5→8），襟翼升力隨偏角減小反而升高，然而受到襟翼俯仰運動的影響，總升力減??；襟翼后緣分離區(qū)逐漸減小，直至重新達到附著狀態(tài)。endprint

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

3 結(jié)語

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

參考文獻

[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] 肖志祥.復雜流動Navier-Stokes方程數(shù)值模擬及湍流模型應(yīng)用研究[D].西安：西北工業(yè)大學，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] 錢煒祺，符松，蔡金獅.翼型動態(tài)失速的數(shù)值研究[J].空氣動力學學報，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