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

趙文梅，趙 軍，胡 偶

(中航工業(yè)直升機(jī)設(shè)計(jì)研究所，江西 景德鎮(zhèn) 333001)

?

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

趙文梅，趙 軍，胡 偶

(中航工業(yè)直升機(jī)設(shè)計(jì)研究所，江西 景德鎮(zhèn) 333001)

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

混合笛卡爾網(wǎng)格；自適應(yīng)；非定常；氣動(dòng)力；旋渦-翼型相互作用

0 引言

直升機(jī)旋翼在轉(zhuǎn)動(dòng)過程中槳葉的槳尖處會(huì)產(chǎn)生非常強(qiáng)烈的尾跡渦結(jié)構(gòu)，且隨著尾跡渦的發(fā)展會(huì)與后續(xù)槳葉發(fā)生相互作用，即直升機(jī)空氣動(dòng)力學(xué)中一個(gè)重要的研究課題：旋渦-槳葉相互作用問題(Blade-Vortex Interaction, BVI)。BVI問題對(duì)直升機(jī)的旋翼振動(dòng)、氣動(dòng)噪聲等都有顯著的影響。目前，針對(duì)BVI問題，為簡(jiǎn)化問題模型和發(fā)展相關(guān)數(shù)值方法，一般將其簡(jiǎn)化為二維的旋渦-翼型相互作用問題(Airfoil-Vortex Interaction, AVI)。通過對(duì)AVI問題的研究，發(fā)展相應(yīng)的數(shù)值模擬方法，揭示旋渦的發(fā)展過程以及旋渦與翼型相互作用過程中引起的復(fù)雜流動(dòng)現(xiàn)象和氣動(dòng)力變化。

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

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

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

1 AVI問題幾何與數(shù)學(xué)模型

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

圖1 AVI問題幾何模型

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

根據(jù)試驗(yàn)測(cè)量結(jié)果，Scully旋渦模型被用于定義數(shù)值模擬中的旋渦結(jié)構(gòu)。在極坐標(biāo)系(r,θ)下，Scully旋渦模型的切向速度Xθ的表達(dá)式為：

其中，r為網(wǎng)格點(diǎn)到渦核的距離，rc表示渦核半徑。另外，還需要給出壓強(qiáng)和密度關(guān)系，由于試驗(yàn)中旋渦是由激波衍射產(chǎn)生的，因此根據(jù)等焓假設(shè)和徑向動(dòng)量方程[3]，可得壓力關(guān)于半徑的解析表達(dá)式：

其中

對(duì)于(2)式，當(dāng)r→∞時(shí)，p(r)→p∞。密度則可由等焓假設(shè)[3]計(jì)算得到。

2 AVI問題數(shù)值模擬方法

精確模擬AVI問題，需要重點(diǎn)關(guān)注的兩大問題是：旋渦位置的捕捉和強(qiáng)度的保持，以及旋渦在與翼型相互作用過程中產(chǎn)生的復(fù)雜流動(dòng)現(xiàn)象的準(zhǔn)確模擬。因此，本文基于雷諾平均Navier-Stokes 方程(RANS)，發(fā)展了一套能夠根據(jù)旋渦特征跟蹤旋渦結(jié)構(gòu)和動(dòng)態(tài)自適應(yīng)分布計(jì)算網(wǎng)格的自適應(yīng)混合笛卡爾網(wǎng)格方法(AHCG)。該方法采用混合笛卡爾網(wǎng)格系統(tǒng)創(chuàng)建初始計(jì)算網(wǎng)格，在流場(chǎng)計(jì)算過程中，采用基于流場(chǎng)特征的網(wǎng)格自適應(yīng)技術(shù)動(dòng)態(tài)分布網(wǎng)格密度，從而達(dá)到流場(chǎng)精確數(shù)值模擬的目的。

2.1 混合笛卡爾網(wǎng)格方法

如圖2所示，混合笛卡爾網(wǎng)格是由物面附近的貼體結(jié)構(gòu)網(wǎng)格和填充流場(chǎng)其余區(qū)域的笛卡爾網(wǎng)格構(gòu)成，兩套網(wǎng)格之間信息傳遞是通過兩套網(wǎng)格之間的網(wǎng)格交接面完成的。

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

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

(4)

其中，dIj表示網(wǎng)格單元j的中心相對(duì)于網(wǎng)格邊IJ中心的位移向量，nI為網(wǎng)格邊IJ的單位法向量。

圖2 混合笛卡爾網(wǎng)格示意圖

圖3 交接面“貢獻(xiàn)單元”的確定

2.2 基于流場(chǎng)特征的網(wǎng)格自適應(yīng)技術(shù)

在流場(chǎng)數(shù)值模擬過程中可以根據(jù)流場(chǎng)特征，對(duì)流場(chǎng)計(jì)算網(wǎng)格進(jìn)行動(dòng)態(tài)的加密或粗化，達(dá)到流場(chǎng)的精確模擬，即基于流場(chǎng)特征的網(wǎng)格自適應(yīng)技術(shù)。本文以速度散度和旋度為解自適應(yīng)判據(jù)，其表達(dá)式如下：

然后就可以采用如下的判斷標(biāo)準(zhǔn)，確定哪些網(wǎng)格單元需要加密，哪些網(wǎng)格單元需要粗化：

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

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

采用如上的判據(jù)就可以實(shí)現(xiàn)計(jì)算網(wǎng)格隨著流場(chǎng)特征的變化而動(dòng)態(tài)的自適應(yīng)，從而達(dá)到更為準(zhǔn)確地捕捉流場(chǎng)信息的目的。

2.3 控制方程與數(shù)值方法

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

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

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

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

3 AHCG方法數(shù)值驗(yàn)證

以二維非定常繞圓柱層流流動(dòng)問題為例，分析采用AHCG方法和傳統(tǒng)靜態(tài)網(wǎng)格方法，計(jì)算由于非定常脫落渦引起的氣動(dòng)力變化方面的差距。基于來流速度和圓柱直徑的雷諾數(shù)為200，來流馬赫數(shù)為0.3。計(jì)算域?yàn)閇-5.0,15.0]×[-5.0,5.0]，圓心坐標(biāo)為(0.0,0.0)。

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

表1 圓柱繞流計(jì)算結(jié)果比較

圖4 升阻力系數(shù)隨時(shí)間的變化比較

4 AVI問題數(shù)值模擬結(jié)果與分析

4.1 低速M(fèi)a∞=0.5的AVI問題

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

圖5顯示了旋渦接近翼型以及與翼型發(fā)生碰撞的過程中，不同時(shí)刻的自適應(yīng)計(jì)算網(wǎng)格和密度云圖，并與試驗(yàn)觀測(cè)結(jié)果[1]進(jìn)行了對(duì)比。可以看出，當(dāng)順時(shí)針旋渦接近翼型前緣時(shí)，翼型前緣出現(xiàn)下洗速度，導(dǎo)致駐點(diǎn)上移。隨后旋渦與翼型發(fā)生碰撞，旋渦沿翼型下翼面移動(dòng)并在前緣附近誘導(dǎo)引起“袋狀”超音速流動(dòng)區(qū)，之后在翼型的下翼面形成一道弧形激波。當(dāng)旋渦通過激波之后，會(huì)產(chǎn)生一個(gè)逆時(shí)針旋轉(zhuǎn)的二次渦結(jié)構(gòu)。隨著這一對(duì)旋渦向下游移動(dòng)，駐點(diǎn)位置開始下移，并伴有壓縮波產(chǎn)生并隨之向四周擴(kuò)散。

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

4.2 高速M(fèi)a∞=0.8的AVI問題

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

圖5 旋渦與翼型碰撞過程中不同時(shí)刻瞬態(tài)網(wǎng)格，以及密度云圖與試驗(yàn)結(jié)果的比較

圖7給出了翼型升力系數(shù)與旋渦渦核位置的之間的關(guān)系曲線，并與相關(guān)文獻(xiàn)的研究結(jié)果進(jìn)行了比較。數(shù)值結(jié)果顯示，本文方法能夠準(zhǔn)確地模擬隨著旋渦的運(yùn)動(dòng)而引起升力變化的過程。順時(shí)針的旋渦在接近翼型的過程中產(chǎn)生下洗作用，從而導(dǎo)致初始升力值為負(fù)；而隨著旋渦的運(yùn)動(dòng)，當(dāng)旋渦靠近翼型前緣位置時(shí)，升力系數(shù)達(dá)到最小值，但隨后升力系數(shù)迅速恢復(fù)；在旋渦沿著翼型下表面?zhèn)鞑r(shí)，升力系數(shù)轉(zhuǎn)變?yōu)檎怠＿@種升力的變化，對(duì)三維槳葉的受力也具有顯著的影響，因此，在槳葉設(shè)計(jì)過程中需要關(guān)注槳渦干擾問題。

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

5 結(jié)論

本文針對(duì)直升機(jī)空氣動(dòng)力學(xué)中的旋渦-槳葉相互作用問題(BVI)，以其簡(jiǎn)化模型旋渦-翼型相互作用問題(AVI)為研究對(duì)象，發(fā)展了基于雷諾平均Navier-Stokes方程和自適應(yīng)混合笛卡爾網(wǎng)格的流場(chǎng)數(shù)值模擬方法。通過給定合適的自適應(yīng)判據(jù)，動(dòng)態(tài)網(wǎng)格自適應(yīng)技術(shù)能夠?qū)崿F(xiàn)網(wǎng)格的動(dòng)態(tài)分布，進(jìn)而準(zhǔn)確捕捉旋渦位置并保持旋渦強(qiáng)度。

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

圖7 翼型升力系數(shù)與渦核位置的關(guān)系

[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-)，女，甘肅敦煌人，碩士，助理工程師，主要研究方向：直升機(jī)旋翼動(dòng)力學(xué)。

1673-1220(2015)01-006-07

V211.52

A