基于動網(wǎng)格法的翼型啟動過程數(shù)值模擬

沈如松，袁書生，王 強，彭紹雄

（海軍航空工程學院飛行器工程系，山東煙臺264001）

工程問題中存在許多翼型啟動引起的非定常流動現(xiàn)象。目前，公開發(fā)表的針對啟動引起的流動研究主要集中于2 個方面，一是極限狀態(tài)下物體的突然啟動[1-4]，即假設只在t=0時刻存在一個加速度，這大大簡化了理論問題的描述，但它只能在理論分析和數(shù)值計算中適用，而工程實踐中的啟動必然存在加速過程；二是翼型大攻角啟動繞流場的流場結構揭示，如Soria J 和New T H 等[5]采用PIV 技術對繞NACA0015 翼型30°攻角加速流動進行了實驗測量，通過對不同加速度下的流動測試發(fā)現(xiàn)，啟動過程剪切層強度、渦結構的發(fā)展與加速度有密切關系。

文獻[6]對以固定加速度加速流過NACA 0015 翼型的空氣流動進行了參數(shù)化研究，采用流動顯示技術提供了復雜渦結構圖像。

吳大轉等人[7]采用變速度入口法對繞30°大攻角NACA0015 翼型的粘性不可壓縮加速流動進行了二維數(shù)值模擬，給出了2 種加速度下從靜止到穩(wěn)定流速的加速過程和完成加速后的詳細瞬態(tài)流動結構和演化過程。

工程實踐中，為避免失速發(fā)生，飛行器的啟動攻角不能過大，因而有必要研究中小攻角下的翼型啟動流場。另外，啟動過程中翼型的非定常非線性氣動力特性也是工程上關心的問題。本文采用動邊界法對15°攻角下NACA4421翼型的啟動過程進行了數(shù)值計算，與對應的實驗結果進行比較，并分析啟動過程中的氣動力特性。

1 數(shù)值方法

本文所研究的二維非定常流動問題的控制方程組可描述為連續(xù)性方程和動量方程。

①連續(xù)性方程：

②動量方程：

式（2）、（3）中：ρ、μ、u和v分別是空氣的密度、粘性系數(shù)和速度分量。

本文求解的是從靜止到低雷諾數(shù)湍流的過渡流動，湍流模型采用基于雷諾平均的Spalart-Allmaras模型。

對控制方程組采用有限體積法進行計算，將其表達為以下通用的離散形式：

式（4）、（5）中：?為通用變量，代表u和v，下標M是標號為P的單元周圍的網(wǎng)格序號。

時間導數(shù)項歸入到和源項中。離散后的控制方程組采用SIMPLE方法進行求解。

2 網(wǎng)格劃分

取弦長為1 m 的NACA4421 翼型，t=0時刻，將計算域取為長40 m、寬10 m的矩形，將翼型置于該計算域的中心位置。在翼型表面以0.01 m 的間距劃分網(wǎng)格，在外邊界上以0.5 m的間距劃分網(wǎng)格，控制體內部以陣面推進法生成自外向內逐漸加密的三角形非結構網(wǎng)格。翼型在加速運動過程中的動態(tài)控制體網(wǎng)格采用彈簧近似法和網(wǎng)格變形重構相結合的方法繪制。其中，為保證翼型表面附近有足夠的網(wǎng)格密度，取彈簧倔強系數(shù)為0.9；在翼型運動引起網(wǎng)格變形過程中，若網(wǎng)格拉伸長度達到0.05 m 或壓縮長度超過0.001 m，則進行網(wǎng)格重繪。翼型附近的初始網(wǎng)格劃分如圖1 所示。翼型表面邊界條件設置為無滑移的壁面，控制體外邊界上的壓強設為1 標準大氣壓的壓力出口。

圖1 初始網(wǎng)格示意圖Fig.1 Sketch of the initial mesh

3 啟動流場結構

考慮翼型在t=0時刻開始以15°攻角的姿態(tài)和1 m/s2的加速度做勻加速前進，在t=3 s時刻加速至3 m/s，之后維持勻速直線運動，對應雷諾數(shù)為2.1×105。

t=0.05 s翼型繞流場結構如圖2 a）所示，此時，翼型速度0.05 m/s。圖中清晰地顯示翼型剛剛開始啟動時，前緣附近下翼面的流體質點繞過前緣流至上翼面，后緣附近下翼面的流體質點則繞過后緣流至上翼面，形成2個較為明顯的環(huán)流，且這2個環(huán)流都附著在翼型表面，并未在脫離壁面的流場中形成明顯的渦核。圖2 b）是t=1 s翼型繞流場結構圖。此時，翼型速度1 m/s。從中可見隨著翼型的前進，原來附著在后緣的環(huán)流逐漸脫離開尾緣，在翼型后方形成一個渦核，但前緣的環(huán)流仍附著在翼型表面上。圖2 c）是t=3 s翼型繞流結構圖。此時，翼型速度3 m/s，且自此開始維持此速度做勻速直線運動。從中可見，尾緣脫落的渦繼續(xù)留在原地做渦旋運動。而同前3 幅圖相比，原來繞翼型前緣形成的渦已經(jīng)逐漸發(fā)展為繞整個翼型的環(huán)流。圖2 d）是翼型啟動結束后流場結構的實驗結果。對比可見，本文計算結果與實驗結果吻合良好。

圖2 流場結構Fig.2 Flow structure

4 氣動性能

4.1 翼面壓強分布變化

各時刻上下翼面的壓力分布情況如圖3所示。

圖3 壓力分布Fig.3 Pressure distribution

由圖3 a）可以看出，翼型啟動的瞬間，下翼面最大正壓出現(xiàn)在4.9%弦長處，而92.5%弦長至后緣點之間部分的下翼面則提供負壓，不利于升力產生；上翼面最大負壓出現(xiàn)在76%弦長處，而靠近前緣的部分上翼面提供正壓，也不利于升力產生。由圖3 b）～d）可以看出，隨著啟動過程的進行，下翼面的最大正壓位置無明顯變化，仍位于前緣附近，但壓強分布曲線逐漸轉平，在下翼面的大部分區(qū)域產生較為均勻的正壓，且整個下翼面無負壓出現(xiàn)，利于升力的產生。上翼面的最大負壓位置則逐漸向翼型前緣移動。

圖4是上下翼面最大壓差位置隨時間變化圖。從圖中看出，翼型啟動的瞬間，上下翼面最大壓強差出現(xiàn)在66%弦長位置，隨啟動過程進行，最大壓差位置逐漸前移至前緣點附近。

圖4 最大壓差位置Fig.4 Location of maximum pressure difference

上下翼面壓差關于弦向位置的變化率隨時間的變化如圖5所示。

圖5 上下翼面壓差變化率Fig.5 Pressure difference changing rate

從圖5 中可以看出整個翼型啟動過程中，最大壓差變化率均出現(xiàn)在前緣點附近，而翼型中部的絕大部分區(qū)域壓差變化率分布則較為均勻，且隨著啟動過程的進行，最大壓差變化率迅速增加。因此，工程實踐中，應重點加強翼型前緣附近的結構強度。

4.2 升力變化

整個啟動過程中，翼型升力隨時間的變化情況如圖6所示。

圖6 升力曲線Fig.6 Lift curve

由圖6 可見，翼型在啟動加速階段（前3 s），升力近似以速度的二次冪曲線的規(guī)律上升。本文計算了多個加速度啟動的工況，并對加速啟動階段的升力變化數(shù)據(jù)做關于速度的二次多項式擬合，對多項式中的系數(shù)做關于加速度的三次多項式擬合，最終得到近似公式

式（7）中：L為升力；a為加速度；V為瞬時速度大小。

從式（7）中看出，a2項及a3項在V2項和V項系數(shù)表達式中的權重均分別比a項小1和2個量級。因此，當a較小時，V2項系數(shù)和V項系數(shù)均與加速度呈近似的線性關系；而當a較大時，則不能忽略a2項和a3項的影響。該式適用于NACA4421翼型由靜止開始以15°攻角做加速度在1 m/s2到10 m/s2之間的勻加速啟動過程。

第3 s末，翼型不再做加速運動，而是開始以3 m/s的速度做勻速直線運動。翼型在第3 s前后受到的升力如表1所示。

表1 3 s前后翼型升力Tab.1 Lift on the airfoil around 3 second

由表1 可見，在加速階段剛剛結束的第1 個時間步內，即由第3 s至第3.05 s，升力有小幅下降，從7.7 N降低至7.55 N，之后逐漸回升。計算得到勻速段的平均升力系數(shù)為1.38，與實驗值[8]1.35 比較相差不大，驗證了本文計算結果的準確性。

5 結論

針對15°攻角下NACA4421 翼型的啟動過程，本文采用動邊界法對其進行了研究，獲得了流場結構的生成和發(fā)展過程及啟動過程中翼型的非定常非線性氣動力特性。計算結果表明，在流場結構方面，翼型的啟動過程同時也是尾緣啟動渦的生成、脫落及繞翼型環(huán)流充分發(fā)展的過程。在氣動力方面，啟動瞬間，上翼面最大負壓出現(xiàn)在76%弦長處，上下翼面最大壓強差出現(xiàn)在66%弦長位置，隨啟動過程的進行，二者均逐漸向前緣移動，最終穩(wěn)定在前緣點附近。整個啟動過程中，下翼面最大正壓點和上下翼面壓差隨弦向位置的最大變化率均維持在前緣點附近，因而工程實踐中應重點加強翼型前緣附近的結構強度。加速過程中整個翼型受到的升力近似以速度的二次冪曲線的規(guī)律變化。加速段結束后翼型轉入勻速運動的瞬間則出現(xiàn)升力小幅下降的現(xiàn)象，之后逐漸回升至穩(wěn)定升力。

[1] OHUE H，KAWASHIMA G. Experimental study on the flow induced by a suddenly accelerated concentric rotating drum[J]. Transactions of the Japan Society of Mechanical Engineers：Part B，2003，69（7）：1583-1589.

[2] 袁新，江學忠. 翼型大攻角低速分離流動的數(shù)值模擬[J].工程熱物理學報，1999，20（2）：161-165.

YUAN XIN，JIANG XUEZHONG. Simulation of lowspeed separated flows for the NREL airfoil at high angles of at- tack[J]. Journal of Engineering Thermophysics，1999，20（2）：161-165.（in Chinese）

[3] 王演興，宛輝，胡坤忠.繞大攻角翼型低速非定常流動的SUPG 有限元數(shù)值模擬[J]. 空氣動力學學報，1999，17（3）：309-313.

WANG YANXING，WAN HUI，HU KUNZHONG. Numerical simulation of the unsteady flow past an airfoil at high incidence by SUPG FEM[J].Acta Aerodynamica Sinica，1999，17（3）：309-313.（in Chinese）

[4] SCANDURA P. Two-dimensional perturbations in a suddenly blocked channel flow[J]. European Journal of Mechanics B/Fluids，2003，22（2）：317-329.

[5] SORIA J，NEW T H，LIM T T. Multigrid CCDPIV measur-ements of accelerated flow past an airfoil at an angle of attack of 30[J]. Experimental Thermal and Fluid Science，2003，27（4）：667-676.

[6] FREYMUTH P. Vortex patterns of dynamic separation：a parametric and comparative study[J]. Progress in Aerospace Sciences，1985，22（3）：161-208.

[7] 吳大轉，王樂勤. 繞大攻角翼型加速流動的數(shù)值模擬[J].工程力學，2006，23（12）：68-72.

WU DAZHUAN，WANG LEQIN. Numerical simulation of accelerated flow past an airfoil at high angles of attack[J]. Engineering Mechanics，2006，23（12）：68-72.（in Chinese）

[8] JOHN D ANDERSON JR.Fundamentals of aerodynamics[M]. 3rd ed. New York：The McGraw-Hill Companies，Inc.，2001：332.

[9] 于沖，王旭.y+值對翼型氣動參數(shù)計算精度的影響研究[J]. 空軍工程大學學報：自然科學版，2012，13（3）：25-29.

YU CHONG，WANG XU. The study of effect of y+ on precision of pneumatic parameters of foil[J]. Journal of Air Force Engineering University：Natural Science Edition，2012，13（3）：25-29.（in Chinese）.

[10]潘光，鐘如意.撲翼翼型流體動力特性數(shù)值計算與分析[J].魚雷技術，2012，20（1）：9-13.

PAN GUANG，ZHONG RUYI. Numerical computation and analysis on hydrodynamic characteristics of different profiles of flapping foil[J].Torpedo Trchnology，2012，20（1）：9-13.（in Chinese）.

[11]LICHT S，POLIDORO V，F(xiàn)LORES M，et al.Design and projected performance of a flapping foil AUV[J]. IEEE Journal of Oceanic Engineering，2004，29（3）：786-794.

[12] ANDERSON J D. Aircraft performance and design[M].Boston：The McGraw-Hill Companies，Inc.，1999：28-39.