非定常旋渦-翼型相互作用問題數值模擬研究

趙文梅，趙 軍，胡 偶

(中航工業直升機設計研究所，江西 景德鎮 333001)

?

非定常旋渦-翼型相互作用問題數值模擬研究

趙文梅，趙 軍，胡 偶

(中航工業直升機設計研究所，江西 景德鎮 333001)

針對直升機旋渦-槳葉相互作用問題，以其二維簡化模型旋渦-翼型相互作用問題為研究對象，發展了基于雷諾平均Navier-Stokes方程和自適應混合笛卡爾網格的流場數值模擬方法，通過給定合適的自適應判據，動態網格自適應技術能夠準確捕捉旋渦位置和保持旋渦強度。通過兩類旋渦-翼型相互作用問題的數值模擬研究，揭示了旋渦與翼型相互作用過程中產生的復雜流動現象和非定常氣動力變化。數值模擬結果與相關文獻中的試驗數據和數值計算結果吻合，且對三維槳葉的氣動載荷研究和振動分析具有指導意義。

混合笛卡爾網格；自適應；非定常；氣動力；旋渦-翼型相互作用

0 引言

直升機旋翼在轉動過程中槳葉的槳尖處會產生非常強烈的尾跡渦結構，且隨著尾跡渦的發展會與后續槳葉發生相互作用，即直升機空氣動力學中一個重要的研究課題：旋渦-槳葉相互作用問題(Blade-Vortex Interaction, BVI)。BVI問題對直升機的旋翼振動、氣動噪聲等都有顯著的影響。目前，針對BVI問題，為簡化問題模型和發展相關數值方法，一般將其簡化為二維的旋渦-翼型相互作用問題(Airfoil-Vortex Interaction, AVI)。通過對AVI問題的研究，發展相應的數值模擬方法，揭示旋渦的發展過程以及旋渦與翼型相互作用過程中引起的復雜流動現象和氣動力變化。

AVI問題涉及非常復雜的流動現象，例如旋渦的傳播、激波的形成、二次渦的產生、旋渦碰撞、以及聲波傳播等等。諸多學者采用計算流體動力學(CFD)技術開展了AVI問題的數值模擬研究，主要有：Lee和Berchader[1]的高精度數值方法，Falissard和Lerat[2]提出的渦量保持格式(vorticity-preserving schemes)，Morvant[3]的渦量限制格式(vorticity-confinement schemes)，Lei Tang和Bader[4]基于三階精度的TVD格式并結合網格重分布方法，Woo等人[5-6]的非結構網格自適應技術等。

傳統的基于靜態網格的CFD流場數值模擬技術，在處理AVI問題時，由于網格數目相對固定，因而不能準確地捕捉由于旋渦的運動引起的非定常流動現象，包括旋渦強度的保持、位置的捕捉以及碰撞過程中產生的復雜流動現象等。為此，本文基于雷諾平均Navier-Stokes方程，針對旋渦運動問題，發展了一套自適應混合笛卡爾網格(Adaptive Hybrid Cartesian Grid, AHCG)流場數值計算方法。

本文首先給出了AVI問題的幾何模型；其次，建立AVI問題的數值模擬方法；最后，將發展的數值方法用于AVI問題的數值模擬研究。重點分析了旋渦的發展歷程，旋渦與翼型相互作用過程中復雜流場現象，以及作用在翼型上的氣動力變化。

1 AVI問題幾何與數學模型

AVI問題的試驗結果來源于Lee和Bershader[1]等人的研究，通過激波與NACA0018翼型衍射作用形成旋渦結構，然后與下游的NACA0012翼型碰撞。AVI問題的幾何模型如圖1所示。

圖1 AVI問題幾何模型

其中，(ΔX,ΔY)為旋渦的放置位置，Γ為旋渦的強度，V∞為來流速度，c為翼型弦長。

根據試驗測量結果，Scully旋渦模型被用于定義數值模擬中的旋渦結構。在極坐標系(r,θ)下，Scully旋渦模型的切向速度Xθ的表達式為：

其中，r為網格點到渦核的距離，rc表示渦核半徑。另外，還需要給出壓強和密度關系，由于試驗中旋渦是由激波衍射產生的，因此根據等焓假設和徑向動量方程[3]，可得壓力關于半徑的解析表達式：

其中

對于(2)式，當r→∞時，p(r)→p∞。密度則可由等焓假設[3]計算得到。

2 AVI問題數值模擬方法

精確模擬AVI問題，需要重點關注的兩大問題是：旋渦位置的捕捉和強度的保持，以及旋渦在與翼型相互作用過程中產生的復雜流動現象的準確模擬。因此，本文基于雷諾平均Navier-Stokes 方程(RANS)，發展了一套能夠根據旋渦特征跟蹤旋渦結構和動態自適應分布計算網格的自適應混合笛卡爾網格方法(AHCG)。該方法采用混合笛卡爾網格系統創建初始計算網格，在流場計算過程中，采用基于流場特征的網格自適應技術動態分布網格密度，從而達到流場精確數值模擬的目的。

2.1 混合笛卡爾網格方法

如圖2所示，混合笛卡爾網格是由物面附近的貼體結構網格和填充流場其余區域的笛卡爾網格構成，兩套網格之間信息傳遞是通過兩套網格之間的網格交接面完成的。

與傳統的混合網格方法類似，混合笛卡爾網格在笛卡爾網格與貼體網格之間需要進行數據傳遞。網格交接面的信息的傳遞直接影響計算的守恒性、數值模擬的精度以及計算效率。如圖3所示，圖中的實線表示笛卡爾網格的交接面，虛線表示貼體網格的交接面，當需要計算通過笛卡爾網格交接面上的網格邊IJ的數值通量時，需要知道網格邊IJ的左右狀態的信息，而其左側信息可以由笛卡爾網格單元I提供；右側信息則從貼體網格中尋找一個“貢獻單元J(donor cell)”提供。“貢獻單元”滿足需如下的幾何條件：

|dIR|=min(|dIj|),且dIj與nI之間的夾角小于60°

(4)

其中，dIj表示網格單元j的中心相對于網格邊IJ中心的位移向量，nI為網格邊IJ的單位法向量。

圖2 混合笛卡爾網格示意圖

圖3 交接面“貢獻單元”的確定

2.2 基于流場特征的網格自適應技術

在流場數值模擬過程中可以根據流場特征，對流場計算網格進行動態的加密或粗化，達到流場的精確模擬，即基于流場特征的網格自適應技術。本文以速度散度和旋度為解自適應判據，其表達式如下：

然后就可以采用如下的判斷標準，確定哪些網格單元需要加密，哪些網格單元需要粗化：

1)如果一個網格單元滿足τcI>σc或者τdI>σd，那么這個網格就需要加密；

2)如果一個網格單元滿足τcI<0.1·σc并且τdI<0.1·σd，那么這個網格就需要粗化。

采用如上的判據就可以實現計算網格隨著流場特征的變化而動態的自適應，從而達到更為準確地捕捉流場信息的目的。

2.3 控制方程與數值方法

考慮可壓縮雷諾平均Navier-Stokes方程(RANS)的積分守恒形式如下：

其中，W為守恒變量，Fc和Fv分別為對流通量和粘性通量[7]。在任意的有限控制體ΩI上對方程(7)進行空間離散可以得到如下的半離散格式：

上式中，m表示網格單元I和J之間的網格面，NF表示包圍控制體ΩI的網格面的總數目，nm和ΔSm為網格面m的單位外法向矢量和面積。

本文采用有限體積法求解方程(7)，對流項使用具有低數值耗散特性的AUSM+格式[8]求解，并通過解的線性重構[7]獲得二階精度，同時采用Venkkatakrishnan限制器[9]抑制解的非物理振蕩；粘性項采用二階中心型格式離散；時間離散采用LU-SGS隱式迭代方法[10]，同時為提高非定常流動問題的時間模擬精度，采用了雙時間步方法(dual-time stepping)[7]。湍流模型采用SSTk-ω湍流模型[11]。

3 AHCG方法數值驗證

以二維非定常繞圓柱層流流動問題為例，分析采用AHCG方法和傳統靜態網格方法，計算由于非定常脫落渦引起的氣動力變化方面的差距。基于來流速度和圓柱直徑的雷諾數為200，來流馬赫數為0.3。計算域為[-5.0,15.0]×[-5.0,5.0]，圓心坐標為(0.0,0.0)。

圖4.a比較了采用本文AHCG方法(with AMR)和靜態網格(without AMR)方法計算得到的阻力系數CD。從結果可以看出，AHCG方法得到的平均阻力系數為1.331，而靜態網格方法得到的平均阻力系數為1.186。圖4.b則比較了采用網格自適應技術(with AMR)和不采用網格自適應技術(without AMR)兩種方法計算得到的升力系數CL，從圖中可以看出，升力系數圍繞零升力線振蕩，采用網格自適應技術計算得到的Strouhal數為0.195，而不采用網格自適應技術計算得到的Strouhal數為0.166。表1將數值計算得到的阻力系數和Strouhal數與試驗結果[12]進行了比較，從結果可以看出，動態網格自適應技術能夠有效改善氣動力參數的數值模擬結果，與試驗結果更為吻合。

表1 圓柱繞流計算結果比較

圖4 升阻力系數隨時間的變化比較

4 AVI問題數值模擬結果與分析

4.1 低速Ma∞=0.5的AVI問題

基于AVI問題的幾何與數學模型，考慮來流馬赫數為Ma∞=0.5的流動狀態，無量綱的渦核半徑和旋渦強度分別為rc/c=0.018，Γ/(V∞c)=-0.283(順時針方向)，V∞為自由來流速度，旋渦的位置參數為ΔX=-5.0,ΔY=0.0。計算網格采用本文第2節介紹的自適應混合笛卡爾網格。

圖5顯示了旋渦接近翼型以及與翼型發生碰撞的過程中，不同時刻的自適應計算網格和密度云圖，并與試驗觀測結果[1]進行了對比。可以看出，當順時針旋渦接近翼型前緣時，翼型前緣出現下洗速度，導致駐點上移。隨后旋渦與翼型發生碰撞，旋渦沿翼型下翼面移動并在前緣附近誘導引起“袋狀”超音速流動區，之后在翼型的下翼面形成一道弧形激波。當旋渦通過激波之后，會產生一個逆時針旋轉的二次渦結構。隨著這一對旋渦向下游移動，駐點位置開始下移，并伴有壓縮波產生并隨之向四周擴散。

圖6給出了翼型上下表面不同位置(x=0.02,x=0.05,x=0.1)的壓力系數隨時間的變化過程，并與試驗測量結果[1]和相關文獻[6]中的數值模擬結果進行了對比。數值結果表明，當旋渦與翼型發生碰撞時，翼型表面的壓力會發生階躍式劇烈變化，且下翼面有壓力脈動產生。這種翼型表面壓力的劇烈變化反應到三維槳葉上，槳葉表面會產生壓力分布的突變，而這種突變對槳葉的振動和強度都會有影響，因此在設計過程中必須加以考慮。

4.2 高速Ma∞=0.8的AVI問題

考慮跨聲速流動狀態，來流馬赫數為0.8，基于翼型弦長的雷諾數為3.6×106，翼型攻角為0°，無量綱的渦核半徑和旋渦強度分別為rc/c=0.05，Γ/(V∞c)=-0.2(順時針方向)，其中c表示翼型的弦長，V∞為自由來流的速度，旋渦的位置參數為ΔX=-5.0,ΔY=-0.26。該算例沒有試驗數據，因此，將本文數值模擬結果與相關文獻[5,6,13,14]的研究結果進行了對比。

圖5 旋渦與翼型碰撞過程中不同時刻瞬態網格，以及密度云圖與試驗結果的比較

圖7給出了翼型升力系數與旋渦渦核位置的之間的關系曲線，并與相關文獻的研究結果進行了比較。數值結果顯示，本文方法能夠準確地模擬隨著旋渦的運動而引起升力變化的過程。順時針的旋渦在接近翼型的過程中產生下洗作用，從而導致初始升力值為負；而隨著旋渦的運動，當旋渦靠近翼型前緣位置時，升力系數達到最小值，但隨后升力系數迅速恢復；在旋渦沿著翼型下表面傳播時，升力系數轉變為正值。這種升力的變化，對三維槳葉的受力也具有顯著的影響，因此，在槳葉設計過程中需要關注槳渦干擾問題。

圖6 AVI問題翼型表面不同位置的壓力系數隨時間的變化

5 結論

本文針對直升機空氣動力學中的旋渦-槳葉相互作用問題(BVI)，以其簡化模型旋渦-翼型相互作用問題(AVI)為研究對象，發展了基于雷諾平均Navier-Stokes方程和自適應混合笛卡爾網格的流場數值模擬方法。通過給定合適的自適應判據，動態網格自適應技術能夠實現網格的動態分布，進而準確捕捉旋渦位置并保持旋渦強度。

通過兩類AVI問題的數值模擬研究，揭示了旋渦與翼型在相互作用過程中產生的復雜流動現象，包括誘導激波與二次渦的產生、非定常氣動力的變化等。數值模擬結果表明，旋渦在與翼型相互作用的過程中，翼型所受的氣動力產生明顯的變化，而這一現象反映到三維槳葉上，即槳葉的氣動力會產生變化。因此，二維AVI問題的數值模擬，對槳葉所受的載荷、振動和強度特性的研究都具有現實的指導意義。

圖7 翼型升力系數與渦核位置的關系

[1] Lee S, Bershader D. Head-on Parallel Blade-Vortex Interaction [J]. AIAA Journal, 1994, 32(1): 16-22.

[2] Falissard F, Lerat A, Sides J. Computation of airfoil-vortex interaction using a vorticity-preserving scheme [J]. AIAA Journal, 2008, 46(7):1614-1623.

[3] Morvant R, Badcock K, Barakos G, et al. Airfoil-vortex interaction using the compressible vorticity-confinement method [J]. AIAA Journal, 2005, 43(1): 63-75.

[4] Tang L, Baeder J D. Adaptive Euler Simulations of Airfoil-Vortex Interaction [J]. International Journal for Numerical Methods in Fluids, 2007, 53(5): 777-792.

[5] Oh W S, Kim J S, Kwon O J. An Unstructured Dynamic Mesh Procedure for 2-D Unsteady Viscous Flow Simulations[R]. AIAA Paper 2002-0121, 2002.

[6] Oh W S, Kim J S, Kwon O J. Numerical Simulation of Two-Dimensional Blade-Vortex Interactions Using Unstructured Adaptive Meshes [J]. AIAA Journal, 2002, 40(3): 474-480.

[7] BLAZEK J. Computational fluid dynamics: principles and applications[M]. Amsterdam: ELSEVIER, 2001.

[8] LIOU M S. A sequel to AUSM: AUSM+ [J]. Journal of Computational Physics, 1996, 129(2):364-382.

[9] VENKATAKRISHNAN V. On the accuracy of limiters and convergence to Steady State Solutions[R]. AIAA Paper 93-0880, 1993.

[10] SHAROV D, NAKAHASHI K. Reordering of 3-D hybrid unstructured grids for vectorized LU-SGS Navier-Stokes calculations[R]. AIAA Paper 97-2102, 1997.

[11] Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[12] Bake W K. Dipole Sound from Cylinders[M]. Mechanics of Flow-Induced Sound and Vibration, 1st Edition, Vol.1 Academic Press, New York, 1996.

[13] Srinivasan G R, McCroskey W J, Beader J D. Aerodynamics of Two-Dimensional Blade-Vortex Interaction [J]. AIAA Journal, 1986, 24(10):1569-1576.

[14] Damodaran M, Caughey D A. Finite-Volume Calculation of Inviscid Transonic Airfoil-Vortex Interaction[J]. AIAA Journal, 1988, 26(11):1346-1353.

Numerical Simulation of Unsteady Airfoil-Vortex Interaction Problem

ZHAO Wenmei, ZHAO Jun, HU Ou

(Aviation Industry Corporation of China Helicopter Research and Development Institute, Jingdezhen 333001, China)

In the background of blade-vortex interaction problem, considering a 2D simplified model, airfoil-vortex interaction (AVI) problem, a novel numerical method for the AVI problem was developed, based on Reynolds average Navier-Stokes equations and adaptive hybrid Cartesian grid method. Given suitable adaptive criterion, the adaptive mesh refinement process could accurately capture the vortex position and keep the vortex intensity. On this basis, both unsteady flow phenomenon and aerodynamic force of two type AVI problems were investigated. Computed results showed good agreement with existing experimental data and computational results. This result has a guiding significance for researching on aerodynamic load and vibration analysis of the three dimensional blade.

hybrid Cartesian grid；adaptive, unsteady；aerodynamic force；AVI

2014-08-29

趙文梅(1986-)，女，甘肅敦煌人，碩士，助理工程師，主要研究方向：直升機旋翼動力學。

1673-1220(2015)01-006-07

V211.52

A