不同B-V頻率下的飛機尾渦數值模擬研究

潘衛軍,姜沿強,張鈺沁

(中國民用航空飛行學院 空中交通管理學院, 四川 廣漢 618307)

0 引言

飛機在飛行中會形成與大氣環境和自身機型特征相關的、強度逐漸減弱的尾渦流場,這會對后機的飛行安全產生威脅,制約機場的起降效率。因此,研究不同氣象條件下的尾流的演化規律對保障飛機安全運行以及提升機場容量具有重要意義[1]。

尾渦的發展可以分成近區階段、擴展近區階段、中、遠區階段和衰減區域5個階段[2],近區階段主要是指在上下翼面壓強差的作用下機翼、尾翼處渦旋的卷起形成。隨后在翼尖渦的主導下,不斷發展形成1對穩定的反向渦對,并在氣象因素的作用下不斷衰減[3]。目前尾渦演化的研究主要依靠激光雷達的探測、數值模擬和流體實驗等手段進行[4]。大量的數值模擬研究填補了尾渦發展規律的諸多空白[5-6],并促進了尾渦快速預測模型的搭建與完善[7]。尾渦演化的模擬是分析尾流遭遇安全和尾流間隔的重要基礎。Gerz等[8]在對商用飛機尾渦的研究中,給出了尾流遭遇模型和尾流安全評價模型。Ruith等[9]通過直接數值模擬研究證明了尾渦耗散過程是一個自身流場穩定性決定的軸對稱模型。Ghalandari等[10]通過對典型機翼的研究發現對于整機尺度而言,機翼顫振行為引起的運動較小。環境因素是尾渦演化的重要影響因素,諸多學者對不同氣象條件下的尾渦演化進行了研究[11-12]。Misaka等[13]通過大渦模擬研究了飛機尾渦從卷起到衰減的演變過程,得到了在不同環境湍流的情況下主翼脫落的渦面詳細卷起的過程以及渦流場的特征參數之間的相關性。周金鑫等[14]基于歐拉多相流模型,研究得出降雨會加速尾渦強度衰減,降低尾流遭遇的危險性。潘衛軍等[15-16]通過對機翼后緣網格加密提高網格質量,使用添加旋轉修正的SST-RC模型模擬全機尾渦場結構等方法對側風影響下的飛機尾渦演變規律進行了深入的探究。在現代科學中,浮力(buoyancy or brunt-v?is?l?,B-V)頻率多應用于層狀自然介質的動力學分析中,其在溫躍層研究之中具有良好的表現[17],也稱浮力振蕩頻率,是流體對垂直位移(例如對流引起的位移)的穩定性的量度。近地大氣情況復雜,分層現象也較為多樣。對于近地飛行而言,大氣分層、風速風向以及湍流耗散率的強弱都對尾渦演變有影響[18]。Greene[19]最先將B-V頻率、湍流和雷諾數組合在同一個尾渦耗散預測模型中,并且預測模型與飛行測試的結果良好吻合。

部分尾渦耗散預測模型考慮了B-V頻率對飛機尾流耗散的影響。但由于大氣分層情況復雜,難以提出一個適用于所有耗散環境的模型。本文中通過數值模擬對其進行研究,使用UDF(user-defined functions)將B-V頻率引入飛機尾流的環境場中。以空客A333為研究對象,使用ANSYS選擇SSTk-ω湍流模型進行飛機尾流數值模擬。由于溫度分層的存在,模擬也考慮了初始溫度對尾渦耗的影響散。以得到不同初始溫度和B-V頻率組合下的運動與衰減特征,探究大氣分層對尾渦演化的影響規律。

1 數值模擬方法

1.1 控制方程

本文中研究集中在尾渦耗散階段,為節約計算資源采用雷諾平均(reynolds-averaged navier-stokes,RANS)方法進行瞬態計算。RANS方程如式(1)和式(2)所示。

(1)

(2)

選取shear-stress transport(SST)k-ω湍流模型[20]進行方程封閉,相比于baseline(BSL)k-ω該模型考慮了湍流剪切應力的傳遞,有著更高的計算精度。適用于旋轉流動等情況。使模型封閉RANS方程。湍流動能k和比耗散率ω的求解如下:

Gk-Yk+Sk+Gb

(3)

Gω-Yω+Dw+Sω+Gωb

(4)

式(3)和式(4)中:Gk為由于平均速度梯度而產生的湍流動能;Gω為ω的產生;σk和σω分別為k和ω的湍流prandtl數;ut為湍流粘度;Yk和Yw分別為k和ω由于湍流引起的耗散;Dw為交叉擴散項;Sk與Sω分別為湍動能項與湍流耗散原項;Gb和Gωb為浮力項的解釋。

1.2 尾渦模型

飛機尾流是指飛機在飛行過程中,由于飛機自身升力的存在,機翼上下表面的壓力差通過對機翼后方氣流的擾動形成的1對反向旋轉的漩渦。飛機尾渦參數一般由初始尾渦環量Γ0,初始渦核半徑r0以及初始渦間距b0,3個基本參數描述。由kutta-joukowski圓柱繞流定理將環量與升力結合起來,如式(5),得到Γ0的計算公式:

(5)

式(5)中:m為飛機質量;g為重力加速度ρ為空氣密度;V為飛行速度。Gerz等[8]研究發現尾渦的初始間距與飛機翼展b長度之間的比例關系與機翼載荷系數有關,并且該比例接近于π/4。初始渦核半徑與初始渦核間距的比值在0.01～0.05之間[9],本文中取0.035,如式(6)所示。

b0=πb/4,r0=0.035b0

(6)

Holz?pfezl研究發現使用hallock-burnham(HB)渦模型進行尾渦耗散快速預測的結果更加符合雷達探測的實際值[7]。如式(7),借助HB模型可以良好表示單個尾渦環量和切向速度的分布,將左右雙渦流場疊加可得飛機尾渦雙渦分布模型:

(7)

式(7)中:vθ表示尾渦的切向速度;Γl、Γr為左渦和右渦的初始環量,其大小等于Γ0;rl、rr為流場中一點到左渦和右渦心的距離。根據Γ0與b0可以進一步定義尾渦演化的特征速度v0與特征時間t0。

v0=Γ0/2πb0,t0=b0/v0

(8)

式(8)中,v0是根據biot-savart定律所得尾渦在相互誘導下的初始下降速度;t0為尾渦以v0下降距離b0所需時間。

本文中以A333作為尾渦演化的對象,其翼展為60.3 m,航空器質量m為最大著陸質量(MLW) 187 000 kg,飛行速度70 m/s。尾渦初始位置處的初始溫度分別設置為263、283、303 K,其對應的相關尾渦參數如表1。

表1 不同初始溫度下的尾渦參數

使用飛機翼展,初始環量和特征時間對物理量進行歸一化處理,得到無量綱化時間t′=t/t0,無量綱化各向位移x′=x/b,h′=h/b,無量綱化環量?！?Γ/Γ0。x為翼展方向上一點距離原點的距離,h為距初始渦核位置的高度差。

(9)

1.3 環境場

當尾流在大氣層中下降時,會受到粘性力與浮力的影響。表示流體分層的參數是B-V頻率N,表示流體粒子從平衡位置位移的假設振蕩頻率,即垂直浮力振蕩固有頻率[19]。對于理想氣體,B-V頻率可以表示為

(10)

式(10)中:Ti表示大氣溫度;h為飛行高度;A為尾流橢積;θ為位溫,表示干氣塊干絕熱壓縮或膨脹至1 000 hpa所具有的溫度;Cp為定壓比熱。

N*=Nt0,283

(11)

如式(11),使用T=283 K時的特征時間t0對B-V頻率無量綱化,選取無量綱B-V頻率N*為0,0.5和1.0的3種情況進行研究。

2 尾渦的演化與驗證

2.1 尾渦演化設置

構建分辨率為0.5 m的結構化網格,以初始尾渦左右2渦渦核連線中點為網格的原點。通過UDF將HB尾渦模型,在ANSYS Fluent中編譯添加到搭建的環境場之中,得到A333的初始切向速度分布,如圖1所示。

圖1 初始切向速度分布

求解器選擇基于壓力基(分離式)的求解器。算法選擇PISO(pressure-implicit with splitting of operators),該算法將SIMPLE(semi-implicit method for pressure linked equations)所需的重復計算轉化為求解壓力修正方程,使修正后的速度滿足連續性方程和動量方程,從而大大減少收斂所需迭代次數,適用于瞬態計算[21]。梯度離散方案選擇基于單元體的最小二乘法插值。為保障計算精度,選擇二階隱式瞬態公式,壓力項、能量項以及其他項的離散格式均選擇二階迎風格式。

飛機尾渦場流體速度遠小于1.0馬赫數,壓縮性效應可以忽略不計,所以氣體選擇為不可壓縮流體。無量綱大氣湍流耗散率ε′=(εb0)1/3/v0設置為0.23(ε為大氣湍流耗散率)[18]。

2.2 數值模擬與雷達探測對比

提取激光雷達探測到的尾渦信息,與T=283 K,N*=0時的數值模擬計算結果相比較,進一步驗證數值模擬方法的可靠性。

圖2為雷達探測得到的單渦渦核附近的速度分布曲線和數值模擬得到的雙渦速度分布曲線。對于單渦而言,雷達探測結果和數值模擬結果都體現出渦核附近的速度最大,遠離渦核的速度逐漸減小的特征。對于無側風影響的數值模擬,左右渦的變化基本對稱,所以在后續結果分析中環量變化只給出了單渦的結果。受雷達掃描周期和角度的限制,對于同一尾渦,通常只能得到3～4個數據。由于雷達探測受到環境場的影響,所以尾渦徑向速度分布與模擬的數值結果不完全相同。

圖2 尾渦速度分布

3 計算結果分析

3.1 渦結構與渦強度的變化

圖3為不同溫度場下保存步長為20 s的渦量切面圖。由圖3(a)(b)(c)可得,不同溫度條件下,相同B-V頻率下的大氣層結穩定性一致,尾渦結構變化趨勢基本相同。從圖3(b)(d)(e)可得,同一溫度,隨著B-V頻率的增大尾渦下沉速率明顯減小。在尾渦演化至t′=1.5附近時,尾渦渦形發生較大變化,N*=1時,尾渦主渦上方出現了反向的小型次級渦,使尾渦渦心間距出現了小幅的減小。這是由于不穩定的大氣層結,可以加速尾渦的耗散,促進渦形的改變。同時B-V頻率的增大會加強尾渦在下沉時所受的慣性振蕩,抑制尾渦下沉,從而加快尾渦耗散,減小尾渦垂向風險。

圖4中無量綱量均以T=283 K時對應的初始量進行歸一化處理。如圖4所示,演化后期各溫度下的環量值趨近一致。不同的大氣分層對初期環量下降的作用較小,差值在0.5%以內。但由于不同溫度下的初始環量不同,難以判斷環量初始下降速率的不同的變化是由溫度還是由初始尾渦強度引起的。為控制變量,選擇尾渦初始環量為443.41 m2/s進行各初始溫度下的補充計算。

圖3 尾渦結構變化

圖4 不同溫度場下的環量變化

由圖5可得,各初始溫度下的尾渦環量變化趨勢差別不大,初始溫度對尾渦強度減弱影響較小。尾渦耗散前期,具有一定的結構穩定性,大氣分層的影響較小,環量的下降主要受兩渦互誘導效應[18]。隨著渦量的擴散以及尾渦強度的下降,大氣分層對耗散的影響逐漸加強,B-V頻率的增大會抑制尾渦強度的減小。t′=3時,初始溫度283 K,無量綱浮力頻率N*=0、0.5、1條件下的尾渦環量分別減小了71.97%、73.87%、75.58%。

圖5 相同初始尾渦強度下的環量變化

3.2 尾渦位移的變化

圖6反映了渦核運動情況,在重力、浮力和空氣黏性力的作用下尾渦呈現出不同程度的下沉。由圖6(a)可得,B-V頻率對尾渦運動的影響較為明顯。圖6(b)表示尾渦橫向位移隨時間的變化,各演化場景下的橫向位移的范圍為-0.24～0.26個翼展?？梢钥闯?B-V頻率對尾渦橫向運動影響較大,而初始溫度對尾渦橫向移動的作用較小。N*=0時兩渦向內側運動,N*=0.5和N*=1兩渦向外側運動。

圖6(c)則體現出初始溫度與B-V頻率共同作用于尾渦的下沉運動。初始溫度越高所對應的尾渦下沉速度越大,但影響程度較小在6.2%以內。受B-V頻率的影響渦核下沉速率明顯減小,t′=3時N*=0所對應的下沉距離為2.7個翼展,為N*=1時下沉距離的2.71倍。

3.3 溫度場的變化

由于N*=1,N*=0 2種情況下溫度分層存在明顯差異,本文中選取初始溫度為288 K,不同B-V頻率進行溫度變化的研究。提取兩渦渦心連線中點左右各6.6倍翼展的溫度變化,結果如圖7所示。渦旋形成初期會擾亂原本的溫度分層,尾渦內部溫度低于環境溫度。隨著尾渦自身的向內向下的渦旋以及尾渦的下沉,會導致N*=1條件下渦核間的溫度顯著大于左右兩渦,整體尾渦溫度大于環境溫度;N*=0條件下渦核間的溫度顯著小于左右兩渦,整體尾渦溫度小于環境溫度。由PROCTOR的三階段耗散模型[22]可知,隨著尾渦長時間的發展,與外界不斷地熱交換,尾渦內部溫度會逐漸趨近環境溫度,直至尾渦消散。

結合圖6和圖7可得,尾渦的橫向位移由于無重力影響,主要受到大氣分層中浮力效應影響。N*=1時,外部溫度小于尾渦內部,渦核間距不斷減小;N*=0時,外部溫度小于尾渦內部,兩渦向外側移動,間距不斷增大。

圖6 尾渦位移變化

圖7 尾渦橫向溫度場變化

3.4 安全性分析

近地面氣象條件復雜,B-V頻率的值較為多變,這是建立尾流間隔預測系統時需要重點考慮的。由前面的環量分析可得,雖然較大的B-V頻率會加快尾渦強度的減小,提高了后機遭遇前機尾渦時的安全性,但是減弱的程度十分有限。圖8為受B-V頻率的影響尾渦下沉速率不同所引起的不同前后機尾流安全間隔的示意圖。淡藍色與淡綠色的區域表示N*=0和N*=1大氣條件下,隨尾渦下沉前機尾渦危險區向下位移的情況?？傻?N*=1條件下小于N*=0條件下的危險區下沉速率,相應的尾流間隔更大。這是由于較大的B-V頻率會導致尾渦下沉大幅度的減慢,促使前機尾流在后機的飛行航路上停留更久,造成后機需要與前機保持更大的安全距離以避免尾流遭遇對后機飛行安全的影響。

圖8 不同B-V頻率下的尾流間隔

4 結論

1) 不同初始溫度主要會引起尾渦初始強度的不同,但對尾渦特征變化的影響較小。B-V頻率對尾渦下沉運動的影響較大。隨B-V頻率的增大,尾渦下沉速率明顯減小,導致前后機的尾流安全間隔增大。

2) 大氣分層主要作用于t′=1.5之后,會對尾渦的結構變化產生較大影響。無量綱B-V頻率N*=1時,誘使尾渦上方出現小尺度的反向渦旋,促進尾渦形變潰散,加快尾渦強度衰減。

3) 尾渦發展前期,由于渦旋作用兩渦間的溫度會出現一定時間的明顯高于或低于尾渦溫度。隨著尾渦的下沉運動,尾渦溫度與外界溫度差值逐漸增大,這會導致兩渦渦核的運動呈現出不同方向。B-V頻率越大,尾渦演化后期的渦核間距越小。