劉薇，宋國萍，褚雙磊，韓博

(1.中國民航大學(xué)空中交通管理學(xué)院，天津300300;2.天津市空管運行規(guī)劃與安全技術(shù)重點實驗室，天津300300)

0 引言

從翼尖脫落的渦會卷起并形成一對反向旋轉(zhuǎn)的漩渦，即翼尖渦，在向下游流動的過程中翼尖渦逐漸形成尾流［1］。不斷發(fā)展的尾渦系統(tǒng)，根據(jù)大氣狀況會持續(xù)很長時間，具有很高的動能，因而對尾隨的航空器造成潛在危險［2］。隨著空中交通流量的持續(xù)增長，在繁忙的終端區(qū)尾流的影響日益明顯，相對保守的尾流間隔已經(jīng)成為妨礙航空運輸系統(tǒng)容量增加的主要因素之一。直升機旋翼槳葉脫落的翼尖渦與葉片之間相互作用會引起噪音和顫振;船舶和潛水艇的螺旋槳推進器產(chǎn)生的翼尖渦也是噪音的主要來源，當(dāng)執(zhí)行反偵查任務(wù)時，希望減小或避免這種噪音。

為了解決上述這些問題，深入分析翼尖渦的特性就顯得尤為重要。目前對于翼尖渦的分析主要有理論、實驗和數(shù)值模擬等手段，其中，數(shù)值模擬是一種精度高、使用靈活、尺度范圍可控、成本低的方法。國內(nèi)外已有許多研究者采用數(shù)值模擬方法對翼尖渦展開研究。Srinivasan等［3］采用Baldwin－Lomax湍流模型對NACA0015和ONERA翼型在亞聲速和跨聲速范圍內(nèi)的翼尖渦的形成進行了數(shù)值研究。結(jié)果發(fā)現(xiàn)，機翼表面壓力與實驗結(jié)果吻合很好，但是并沒有發(fā)現(xiàn)實驗中的負壓峰值，同時指出機翼翼尖的形狀會顯著地影響計算的精度。Dacles－Mariani等［4－6］采用RANS方法的不同模型對NACA0012三維翼型進行了一系列翼尖渦的數(shù)值模擬，并將結(jié)果與Chow等［7－9］在風(fēng)洞實驗中的結(jié)果進行了對比，得出數(shù)值計算中空間離散至少不小于5階，才能得到渦核的核心不至于太離散;但數(shù)值計算結(jié)果較實驗還有較大差別，其中帶旋流修正的S－A模型的方法模擬相對其他模型更準確。Sung－Eun等［10］基于商用Fluent軟件與Chow實驗進行了對比，并對比了不同網(wǎng)格分辨率和不同方法，包括渦粘模型(S－A，RKE，SST)和雷諾應(yīng)力模型(RSTM)。結(jié)果發(fā)現(xiàn)，綜合考慮精度和計算開銷，局部加密的非結(jié)構(gòu)化網(wǎng)格是最為理想的一種網(wǎng)格分布;RSTM模型模擬結(jié)果與實驗最吻合，計算開銷也最大。

本文的研究目的是采用雷諾平均方法(Reynolds Averaged Navier Stokes，RANS)，對雷諾應(yīng)力項分別采用渦粘模式中的一方程模型——S－A模型和二方程模型——Realizable k-ε模型封閉。對NACA0012翼型的近場翼尖渦的形成、發(fā)展和演化進行數(shù)值模擬并與Chow實驗進行對比，同時為計算遠場尾渦流場提供初值、網(wǎng)格及方法的參考，并為其他翼型和整機的模擬計算提供參考。

1 物理模型與計算方法

1.1 幾何模型與網(wǎng)格劃分

基于Ansys軟件SCDM模塊，NACA0012翼型和外流場尺寸按照文獻［7－9］中Chow在NASA Ames研究中心的低速風(fēng)洞的實驗段繪制，如圖1所示。

圖1 幾何模型及外流場計算域Fig.1 Geometry model and computational domain of outer flow field

機翼弦長c=1 m，半翼展0.75c，以機翼前緣為原點，入口(inlet)在上游x=－2.5c處，出口(outlet)在下游x=5c位置，展向z=1c，法向y=－1/3c～1/3c，機翼以迎角α=10°放置在流場中。

網(wǎng)格的尺度和質(zhì)量決定了數(shù)值計算的精度、效率以及得到的湍流脈動信息量的多少，因此網(wǎng)格必須有足夠的空間分辨率和較好的質(zhì)量。本文基于Ansys Icem劃分流場網(wǎng)格，采用結(jié)構(gòu)化六面體網(wǎng)格，包圍機翼的網(wǎng)格采用自適應(yīng)的O－網(wǎng)結(jié)構(gòu)，以提高壁面附近聚集網(wǎng)格點的效率和網(wǎng)格的正交性，如圖2所示。為了顯示清楚，只顯示了部分邊界網(wǎng)格和機翼前后緣的局部O－網(wǎng)網(wǎng)格。按照文獻［2］中提到的綜合前人研究工作，近場區(qū)域網(wǎng)格總數(shù)要達到百萬級，才有可能捕捉到翼尖渦，因此，本文網(wǎng)格分配情況如表1所示。

圖2 流場部分六面體網(wǎng)格和局部O－網(wǎng)放大圖Fig.2 Local hexahedral grid and amplified O－grid

表1 流場網(wǎng)格分配情況Table 1 Grid distribution in flowfield

精確地計算湍流邊界層非常重要。因此在翼面第一層網(wǎng)格達到10－6m，以保證y+小于5，達到粘性底層，如圖3所示，其中縱軸是用壁面摩擦速度u*和運動粘性系數(shù)ν無量綱化的法向位置y+=yu*/ν，橫軸x/c是用弦長c無量綱化的第一層網(wǎng)格點在機翼表面的分布位置。邊界層內(nèi)的網(wǎng)格增長率為1.1。

圖3 機翼表面第一層網(wǎng)格點的法向位置分布圖Fig.3 Normal distribution of the first layer grids along wing surface

1.2 數(shù)值方法與邊界條件

1.2.1 控制方程

笛卡兒坐標(biāo)系下不可壓的雷諾平均 Navier－Stokes方程為［11］:

式中:＜ ＞為系綜平均;xi和 xj(i，j=1，2，3)為笛卡兒坐標(biāo);ρ和ν分別為流體密度和運動粘性;p為等效壓力;fi為質(zhì)量力;為雷諾應(yīng)力項，分別用S－A 模型［12］和帶旋流修正的 k-ε 模型(Realizable k-ε，簡稱 RKE)［13-15］方程封閉。

1.2.2 邊界條件和數(shù)值方法

為了與文獻［7-9］在低速風(fēng)洞中實驗的情況相一致，設(shè)置空氣密度ρ=1.225 kg/m3，基于弦長的雷諾數(shù)Rec=4.6×106，來流速度約為67 m/s。入口和出口的邊界條件為壓力遠場，其他表面和機翼表面為壁面并滿足無滑移條件。

基于密度基隱格式的時間推進求解方法，采用有限體積法進行離散，壓力、動量和能量方程以及湍流擴散項采用二階迎風(fēng)格式，時間項采用一階隱式格式。

2 仿真結(jié)果及分析

2.1 翼面表壓和平均流場仿真結(jié)果

圖4中給出了在展向z=1/3b位置，NACA0012翼型上下翼面表面靜壓系數(shù)分布圖。由圖可以看出，采用兩種計算方法得出的結(jié)果與實驗值均吻合較好。

圖4 在z/b=1/3處翼面靜壓系數(shù)分布圖Fig.4 Surface Cp profile at z/b=1/3

圖5 (a)和(b)中分別給出了不同流向位置(分別為 x/c=0.40，0.50，0.60，0.70，0.80，0.95，2.50，3.00，4.00，4.50)的靜壓系數(shù)Cp的等值線圖。由于翼尖渦的渦核處壓力要小于來流壓力，因此這些圖能夠很好地反映出它的位置、大小和形狀。總體上看，在上翼面存在較大的負壓區(qū)即吸力面，下翼面為正壓區(qū)即壓力面。在這樣的壓差作用下，在翼尖的后緣下游逐漸形成了近似軸對稱形狀的翼尖渦，沒有明顯的耗散，具有相對穩(wěn)定的渦核壓力。進一步分析后發(fā)現(xiàn)，S－A模型比RKE模型體現(xiàn)出更多的渦耗散特性。

圖5 靜壓系數(shù)CP在不同流向位置的等值線圖Fig.5 Contours of Cp at various streamwise locations

為此，進一步分析渦核處靜壓系數(shù)沿下游的分布，如圖6所示。渦核核心是根據(jù)文獻［2］提到的通過在半翼展面內(nèi)的最大渦量對應(yīng)的位置來確定的。從圖中可以看出，總體上，采用兩種方法的計算結(jié)果在翼面上和剛離開翼面(x/c＜1.5)的位置上都與實驗值(Chow的實驗值只測到流向x/c約為1.6的位置)較為吻合，在離開翼面后核心處壓力突然減小，這可理解為離開翼面后翼尖渦的卷起并形成而導(dǎo)致的;S－A模型從下游x/c=1.5以后核心處壓力開始逐漸增高，在x/c=2.75左右出現(xiàn)急劇上揚，而RKE模型計算的渦核核心壓力出現(xiàn)相對穩(wěn)定的趨勢。

圖6 靜壓系數(shù)沿渦核核心線的分布圖Fig.6 Cp along vortex core centerline

2.2 渦核參數(shù)與近場翼尖渦的形成過程

圖7 和圖8給出了渦核核心分別在x-y和x-z平面內(nèi)的分布情況。其中在x-y平面內(nèi)，兩種方法都預(yù)測出實驗中測出的從x/c=0.5到后緣的渦核核心下移，離開翼面以后上移的趨勢;在x-z平面內(nèi)，渦核向內(nèi)弦(z負方向)移動;在Chow用激光測量的發(fā)煙流動顯示中可觀察到這一現(xiàn)象，并將其解釋為主渦與兩次渦的融合而產(chǎn)生的“扭結(jié)(kink)”現(xiàn)象。整體上，RKE模型比S－A模型預(yù)測得更好。

圖8 渦核核心線在x-z平面的分布圖Fig.8 Vortex core centerline location in the x-z plane

圖9 給出了基于RKE模型，用來流速度U∞無量綱化的軸向渦量 ωx=?w/?y－ ?v/?z沿機翼表面的分布圖。由圖可以清晰地看到翼尖渦的形成過程，在機翼前緣附近x/c=0.05和x/c=0.10位置有一對反向旋轉(zhuǎn)的渦，正的渦量來自于吸力面(上翼面)，負的渦量來自于壓力面(下翼面)，這是由于該區(qū)域存在從翼根向翼尖的強烈的展向壓力梯度;在x/c=0.15以后的位置上，吸力面的正渦量很快消失了，取而代之的是形成的新的負渦，此時吸力面壓力小于翼尖處壓力，流體從翼尖處回流到吸力面，新的負渦量與吸力面正渦量相互作用，并將其耗散掉。從x/c=0.30到x/c=1.00位置可以看到有大量的流體從壓力面經(jīng)過翼尖流向吸力面，原來在吸力面的負渦與在壓力面的負渦相互作用，并在x/c=1.10位置卷起成一個完整的翼尖渦。

圖9 不同流向位置軸向渦量的等值線分布圖Fig.9 Contours of axial component of vorticity at various x/c

圖10 給出了用來流速度U∞和弦長c無量綱化的渦粘系數(shù)ν*=ν/(U∞c)沿渦核核心的分布圖。由圖可以看出，兩種方法差異明顯。從x/c=1到x/c=3，S－A模型對應(yīng)的渦粘系數(shù)增大了將近10倍，而RKE模型渦粘系數(shù)沿流向分布接近于零。渦粘性(渦粘系數(shù))是表征湍流度的物理量，翼尖渦形成于翼面，在壓差作用下，壓力面邊界層分離，充分發(fā)展的湍流從翼尖流向吸力面，分離的剪切層在吸力面卷起并形成翼尖渦，一旦形成，翼尖渦核心處開始趨于穩(wěn)定，因此，渦核核心處的湍流會沿流向減小。RKE模型很好地預(yù)測到了這一趨勢，而S－A模型表現(xiàn)出渦核核心處較大的湍流度，意味著翼尖渦會更快耗散，這也印證了圖5的結(jié)論。

圖10 無量綱化的渦粘系數(shù)沿渦核核心分布圖Fig.10 Normalized eddy viscosity along vortex core centerline

3 結(jié)論

本文采用數(shù)值模擬的方法研究了有限翼展的NACA0012翼型產(chǎn)生的近場翼尖渦流場。主要結(jié)論如下:

(1)從表面壓力和平均流場來看，S－A模型和RKE模型在機翼表面的壓力分布差別不明顯，都與實驗值吻合較好，但是在流向x/c＞1的位置，S－A模型出現(xiàn)壓力增大的趨勢，而RKE模型渦核處壓力相對穩(wěn)定;

(2)從渦核位置來看，兩種方法都能不同程度地捕捉到流動顯示實驗中可觀察到的由主渦與兩次渦的融合而產(chǎn)生的“扭結(jié)”現(xiàn)象。RKE模型與實驗結(jié)果更為接近。進一步采用RKE方法模擬了翼尖渦形成和卷起的過程，捕捉到了近場翼尖渦形成并逐漸卷起的過程。

(3)從渦核粘性來看，兩種方法差別明顯，S－A模型體現(xiàn)出了翼尖渦渦核更大的耗散性，從而會低估翼尖渦的強度，而RKE模型模擬的結(jié)果則相對穩(wěn)定，與實驗結(jié)果更符合。

綜上，本文工作表明，采用先進的CFD方法可以較好地預(yù)測翼尖渦流場，這將為進一步研究翼尖渦在下游的發(fā)展演化即遠場尾渦流場提供網(wǎng)格和計算方法的參考依據(jù)，也為其他翼型和整機的數(shù)值模擬提供了參考。

［1］ Rossow V.Lift－generated vortex wakes of subsonic transport aircraft［J］.Progress in Aerospace Sciences，1999，35:507－660.

［2］ Thomas Gerz，F(xiàn)rank Holz?pfel，Denis Darracq.Commercial aircraft wake vortices［J］.Progress in Aerospace Sciences，2002，38:181－208.

［3］ Srinivasan G，McCroskey W，Baeder J，et al.Numerical simulation of tip vortices of wings in subsonic and transonic flows［R］.AIAA－1986－1095，1986.

［4］ Dacles－Mariani J，Rogers S，Kwak D，et al.A computational study of wingtip vortex flowfield ［R］.AIAA－1993－3010，1993.

［5］ Dacles－Mariani J，Zilliac G，Chow J S，et al.Numerical/experimental study of a wingtip in the near－field ［J］.AIAA Journal，1995，33(9):1561－1568.

［6］ Dacles－Mariani J，Kwak D，Zilliac G.On numerical errors and turbulence modeling in tip vortex flow prediction［J］.International Journal for Numerical Methods in Fluids，1999(30):65－82.

［7］ Chow J S，Zilliac G，Bradshaw P.Measurements in the near－field of a turbulent wingtip vortex［R］.AIAA－1993－0551，1993.

［8］ Chow JS.Turbulence measurements in the near－field of a wingtip vortex［D］.Stanford，CA:Ph.D Thesis，Dept.of Mechanical Engineering，Stanford Univ.，1994:20－25.

［9］ Chow J S，Zilliac G，Bradshaw P.Mean and turbulence measurements in the near field of a wingtip vortex ［J］.AIAA Journal，1997，35(10):1561－1567.

［10］ Sung－Eun Kim，Shin Hyung Rhee.Prediction of tip－vortex flow past a finite wing［R］.AIAA－2005－58，2005.

［11］張兆順，崔桂香，許春曉.湍流理論與模擬［M］.北京:清華大學(xué)出版社，2005:24－26.

［12］ Spalart P R，Allmaras S R.A one－equation turbulence model for aerodynamic flows［C］//30th Aerospace Sciences Meeting and Exhibit.Reno，Nevada:AIAA，1992:1－22.

［13］ Jones W P，Launder B E.The prediction of laminarization with a two－equation model of turbulence［J］.International Journal of Heat and Mass Transfer，1972，15(2):301－314.

［14］ Launder B E，Sharma B I.Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disc［J］.Letters in Heat and Mass Transfer，1974，1(2):131－138.

［15］ Wilcox David C.Turbulence modeling for CFD［M］.Second Edition.Anaheim:DCW Industries，1998:174－176.