不同側風類型影響下的飛機尾渦數值模擬研究

潘衛軍，羅昊天，羅玉明，王靖開，姜沿強

(中國民航飛行學院空中交通管理學院，廣漢 618307)

航空器在進入前機尾流區易出現俯仰、滾轉等影響安全的現象。尾流作為影響民航安全的重要因素，盡管對著陸過程有著不良影響，但相關研究表明，側風可以迅速將尾渦吹離跑道反而能減少了飛機進近階段所需的著陸間隔距離，提高了機場的空域容量。側風對尾渦的擾動十分復雜，包括輸送以及對渦強度和衰減率的明顯干擾，由于進近階段各種阻力和湍流的影響，以及溫度變化的影響，復雜和非線性的風切變梯度側風在大氣邊界層的最低高度很常見，因此對側風場中尾流的演變與探測進行更詳細的研究是必不可少的。中外學者借助計算流體力學(computational fluid dynamics, CFD)方法對飛機尾流做了較為深入的研究。Holz?pfel等[1]在維也納機場進行了在飛機尾渦在不同風速、熱分層等條件下的尾渦耗散研究評估當前歐洲尾流間隔重新分類標準(European wake turbulence categories and separation minima on approach and departure, RECAT-EU)以縮減飛機間隔；Zholtovski等[2]通過可壓縮雷諾平均(Reynolds-averaged Navier-Stokes, RANS)解析了矩形 NACA0012幾何機翼渦流周圍的近場及其邊界層問題；Landa等[3]應用剪切應力(shear stress transport, SST)渦流黏度模型和(Speziale-Sarkar-Gatski/Launder-Reece-Rodi, SSG/LRR-ω)雷諾應力模型的數值模擬結果驗證湍流模型捕捉渦流發展的能力；Crouch等[4]采用二階矩RANS CFD方法對尾渦初始耗散階段和演化階段進行了模擬并通過尾渦強度評估機場終端區飛機之間安全間隔；Knopp等[5]提出了一個在顯著的正壓力梯度和高雷諾數下的湍流邊界層流動實驗，使用SSTk-ω(k為湍動能，ω為耗散率)、SSG/LRR-ω和JHh-v2雷諾應力模型進行RANS模擬；文獻[6-8]將SST雙方程k-ω湍流模型和雙方程間歇性-過渡起始動量厚度雷諾數(γ-Reθ)過渡模型在CFD求解器中的傳輸方程一般框架下進行整合進行了關于多級低壓渦輪機中的非穩態湍流和過渡效應的模擬；Imiela等[9]采用SSTk-ω湍流模型和C型網格的組合實現了對空氣動力系數的預測，并提出了一種通過CFD和計算氣動聲學(computational aeoro acoustics，CAA)生成360°翼型極曲線和氣動聲學特性的方法；Maruyama等[10]引入了使用機器學習方法對封閉系數進行統計推斷的框架，以提高CFD在飛行包線限制范圍內可靠湍流建模的預測能力；Sedlacek等[11]通過數值模擬對三角翼和雙三角翼尾渦主導流的網格靈敏度和建模誤差依賴性分析了不同湍流閉合的影響；錢宇等[12]通過建立飛機著陸狀態仿真模型，采用結構化網格對模型進行網格劃分，利用轉捩修正的SSTk-ω湍流模型，通過求解不可壓縮的N-S(Narier-Stokes)方程對生成的網格進行數值計算，得到了著陸狀態下機翼及近翼流場翼尖渦的連續演化過程；谷潤平等[13]為研究尾流特性，降低飛機運行風險，基于數值模擬的研究情況，采用大渦模擬的方法，借助ANSYS軟件對尾流進行仿真模擬；張宇軒等[14]以NACA0012機翼為對象，采用大渦模擬方法，研究了3組不同的馬赫數和雷諾數下翼尖渦中主渦和次級渦的演化特性以及其對機翼氣動力的影響;張鈞鐸等[15]應用升力面模型和自適應網格大渦模擬方法，模擬了國產ARJ21飛機尾渦在3種側風條件下的演化與衰減過程，并對尾渦流場結構特點進行了分析；周金鑫等[16]引入了多相流模型，利用數值模擬方法研究了不同降雨強度條件下，尾渦演化過程中的特征量演變規律；Xu等[17]通過引入地面吹吸區增強機場跑道上的飛機尾渦衰減；Li等[18]通過數值模擬研究均勻側風對流和線性垂直切變側風對流對尾渦強度的影響；Zhou等[19]采用基于解的動態自適應網格方法計算尾渦演化并采用基于自適應網格的大渦模擬(large eddy simulation, LES)對環境湍流強度不同的3種情況進行了研究；潘衛軍等[20]使用添加旋轉修正的SST-RC模型對A330-200飛機進行全機數值模擬研究側風下的翼尖渦耗散。在飛機尾渦探測研究中，K?pp等[21]在塔布斯機場進行的現場試驗，驗證了多普勒相干雷達(coherent Doppler lidar,CDL)實現了尾流從產生的時刻到尾渦衰減階段在長時間內的精確測量；Smalikho等[22]開發了一種通過1.5 μm脈沖相干多普勒激光雷達“流線”進行測量的策略，并提出了一種根據激光雷達數據估計飛機尾渦環量與渦核位置的方法；王筱曄等[23]提出了一種基于CDL頻譜寬度和徑向風速的尾渦快速識別方法,并基于此方法分析了典型機型的尾渦演化過程；Wu等[24]在近地效應(near-ground effect，NGE)下使用脈沖相干多普勒激光雷達(pulsed coherent Doppler lidar, PCDL)評估尾渦特性。為了實時可視化尾渦，開發了尾渦可視化演示器(vortex visualization demonstrator, V2D)。結合徑向速度分布和快速傅里葉變換(fast Fourier transform, FFT)譜表征來識別尾渦的渦核位置。同時基于速度包線和Burnham-Hallock模型修正用于反演NGE下的尾渦環量[24]。因此研究側風下的尾渦演化規律，對CDL激光雷達實地尾渦探測與數值模擬對比驗證有重要指導意義。

目前中外學者已對線性分布的側風和尾渦之間的相互作用有了一定的研究，但非線性分布或均勻的側風究竟是如何影響渦流跡的并沒有得到較好的模擬與驗證，因此對于側風下的激光雷達尾渦探測與數值模擬結果的對比驗證非常重要。為此，建立了空客A330-300的尾渦模型，使用ICEM繪制計算域網格，在ANSYS Fluent中采用SSTk-ω湍流模型進行RANS計算，通過Tecplot和MATLAB對數值模擬后的流場結構進行后處理。在上述研究的基礎上，還研究了均勻或非線性變化的側風中的尾渦的運動。在成都雙流機場進行了CDL尾渦探測實驗，并結合數值模擬結果與實地探測結果進行對比。更加全面地驗證了不同側風影響下的尾渦演化規律。

1 控制方程

數值模擬使用ANSYS Fluent的雙精度求解器進行，控制方程中的非定常項使用二階隱式格式離散化；對流項和擴散項分別使用二階迎風格式和二階中心差分格式進行離散。SSTk-ω模型[25]用于求解不可壓縮流動。

在雷諾平均中，瞬時(精確)Navier-Stokes方程中的解變量被分解為平均值(整體平均或時間平均)和波動分量。對于速度分量ui有

(1)

同樣，對于壓力和其他標量，可表示為

(2)

(3)

湍動能k和比耗散率ω可由式(4)、式(5)獲得。

Gk-Yk+Sk+Gb

(4)

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

(5)

式中：Gk為由于平均速度梯度產生的湍動能；Gω為ω的梯度；Гk和Гω分別為k和ω的有效擴散率；Yk和Yω分別為k和ω由于湍流的耗散；Sk和Sω為用戶定義源項；Gb和Gωb為浮力項。

SSTk-ω模型包括BSLk-ω模型的所有增強，此外在湍流黏度的定義中還考慮了湍流切應力的傳遞。這些特點使SSTk-ω模型對于更廣泛的流量類別(如逆壓梯度流、翼型、跨音速沖擊波)比標準k-ω模型和SSTk-ω模型更準確和可靠。

前面描述的BSL模型結合了Wilcox模型和k-ω模型的優點，但仍無法正確預測從光滑表面開始的邊界層分離和數量。其主要原因是這兩個模型都沒有考慮到湍流剪切應力的傳輸。這導致了渦流黏度的過度預測，適當的傳輸行為可以通過對渦流黏度表述的限制器得到。

(6)

式(6)中：μt為渦流黏度；S′為應變率；常數α1= 0.31；系數α*抑制湍流黏度導致低雷諾數修正；F2可表示為

(7)

式(7)中：μ為流體黏性系數；y為到下一個表面的距離。

2 前期處理工作

2.1 尾渦模型

當飛機垂直方向受力平衡時，飛機尾渦垂直動量等于飛機所受重力，尾渦初始環量可表示為[26]

(8)

式(8)中：Γ0為飛機尾渦初始環量；MLW為飛機最大著陸重量；g為當地加速度；B為翼展；V為飛機速度。

尾渦初始流場采用Burnham-Hallock模型[27]模擬，切向速度Vθ(r)由初始環量Γ0、渦核半徑rc和到渦心距離r決定，其表達式為

(9)

Gerz等[28]對初始渦核半徑rc0通過初始渦核間距b0定義為

(10)

rc0=0.052b0

(11)

時間尺度t0描述了飛機或飛機模型產生下沉的尾渦對向下沉降一個初始渦核間距的時間，可表示為

(12)

式(12)中：s為翼展載荷系數，取值π/4。

無量綱時間t*可表示為

(13)

式(13)中：t為尾渦耗散時間。

(14)

2.2 網格和初始化

選擇A330-300作為研究對象，其機型參數如表1所示。

計算后得到A330-300尾渦參數如表2所示。

表 1 A330-300機型性能和尺寸參數

表 2 A330-300尾渦參數

為提高尾渦數值模擬精度，迎風耗散格式采用二階迎風，湍流模型采用SSTk-ω模型。選取固定位置渦核周圍處渦量，進行的網格無關性驗證如圖1所示，全局網格尺寸0.7 m較1 m網格尺寸特征位置處渦量縮減13.8%；在0.5 m網格尺寸時渦量值較0.7 m網格尺寸時特征位置處渦量縮減3.84%，并在0.3 m網格尺寸時趨于平緩，較于0.5 m的網格尺寸渦量減少1.69%，遠小于5%，符合網格無關性驗證條件，在此基礎上繼續增加網格數量對結果影響較小。

圖1 網格無關性收斂曲線

同時采用基于有限體積法的結構化網格完成計算域網格劃分，相比于非結構化網格，結構化網格更為精細，可以更好地提高參數收斂性與計算精度，設置網格尺寸為0.3 m。以渦核間距中心點為坐標原點;以垂直于地面方向左側為y軸正方向；以速度入口來流方向為x軸正方向。在Fluent UDF中解釋H-B模型，計算域的初始渦核位置和生成的初始尾渦速度場分別如圖2和圖3所示。

（一）疫病防控 將保險所用資金用于疫病防控，除疫苗注射和圈舍定期消毒投資以外，當出現重大動物疫病需要對病死牛羊進行無害化處理時，在條件允許的情況下可以進行補貼。

圖2 計算域和初始渦位置

圖3 初始切向速度剖面圖

2.3 邊界條件和計算參數

風切變效應的非線性函數可采用對數律分布或指數律分布來進行描述，但指數率計算的風速值與實測值偏差較小，且用指數律分布計算風速輪廓線比較簡便。風切變的指數律分布可表示為[29]

(15)

根據民航局CAAC-25文件，對陸上飛機應制定在干跑道上對起飛和著陸演示是安全的90°側風分量，該分量必須至少為20節或0.2倍基準失速速度(0.2VSRO)，取大者，但不必超過25節。因此最大側風風速取25節。

側風切變的垂直梯度可以優先降低渦旋的下降速度，導致它們傾斜，增加分離，有時甚至向上上升，環境側風高度的垂直二階導數S定義為[31-37]

(16)

(17)

式(17)中：S*為無量綱化的S；v0為尾渦初始下降率。

圖4為尾渦對側風切變效應示意圖。

綠色曲線為側風切變梯度為負或正時的側風垂直剖面；左渦由紅色圓圈表示；右渦由藍色圓圈表示；θ為尾渦對的傾角

為研究尾渦在不同側風條件下的演化規律，分別取1、3、7、10 m/s的4種均勻側風、線性側風和非線性側風，速度隨高度變化的具體函數關系如表3所示。流體域環境變量如表4所示。

表3 不同類型側風的函數關系

表4 流體域環境變量

3 仿真結果

3.1 數值模擬結果

尾渦運動示意圖如圖5所示。

圖5 側風下的尾渦運動示意圖

二維Q準則(Qcriterion)可表示為[38]

(18)

式(18)中：u、v分別為x、y的速度方向。

圖6 均勻側風1 m/s

圖7 均勻側風3 m/s

圖8 均勻側風7 m/s

圖9 均勻側風10 m/s

圖10 線性側風(S*=0)

圖11 非線性側風(S*<0)

圖12 非線性側風(S*>0)

對以下變量進行無量綱化處理。

(19)

(20)

(21)

(22)

由圖13可以看出，左右渦橫向偏移量受側風風場大小影響，側風越大偏移量越大。在S*>0或S*<0的垂直非線性側風條件下，橫向偏移量更大。

圖13 不同側風類型下的尾渦無量綱橫向偏移隨無量綱時間變化示意圖

由圖14可知，在任何環境側風條件下，左渦下沉量較右渦大，左右渦下沉趨勢不變。對比分析左右渦縱向運動趨勢發現尾渦的高度變化受不同側風影響較大，側風風速的改變將對尾渦的縱向偏移產生截然不同的運動軌跡和趨勢。這是因為左渦直接和側風接觸，相互作用并融合，導致左渦的空間高度較右渦在垂直方向上有著更強的波動。在施加的持續側風風場擾動下尾渦的速度場將在短時間內發生劇烈改變，尾渦對在空間上的分布也趨于復雜。

圖14 不同側風類型下的尾渦的無量綱下降高度隨無量綱時間變化示意圖

由圖15可以看到，側風風場對尾渦存在明顯的輸運作用，且隨著風速的增加輸運作用更為強烈。在均勻側風和線性側風條件下尾渦渦流跡在等距離X下沉更多，而非線性側風條件下則在等距離Y下橫向側移更多。

圖15 不同側風類型下的渦流跡示意圖

3.2 實地探測數據對比

實地雷達探測在成都雙流國際機場獲得的A330-300著陸時的各參數反演值與假設的A330-300尾流的數值模擬的結果接近。由于尾渦迅速減弱，可能只在有限的時間內被探測到。

圖16為根據H-B模型計算出尾渦初始的切向速度分布與實地探測結果對比。可以看出，由于左渦(順時針)、右渦(逆時針)旋轉，使得尾渦的切向速度帶有方向性，負值代表尾渦切向速度方向為y軸負方向，對比雷達尾渦切向速度變化曲線可以看出，隨徑向距離增加切向速度呈現出先增后減的變化趨勢。尾渦模型在切向速度分布上的差異與探測結果基本一致。

圖16 H-B尾渦模型切向速度分布與CDL雷達實測分布

現場測量通過雷達回波判斷尾渦是否存在并根據回波數據給出尾渦位置、強度等信息。基于CDL觀測到的非線性側風條件下的A330-300型客機尾渦演化的徑向風速圖如圖17、圖18所示。可以看出，對應時刻的尾渦強度較大，正負風速絕對值顯著大于周圍環境風場。據此可以判斷，尾渦已經產生，經過數據處理后，可以判斷為尾渦初始位置。

圖17(e)、圖18(e)分別為為同時段尾渦演化的頻譜寬度組圖，其呈現的尾渦基本演化特征與徑向風速圖和數值模擬徑向速度云圖[圖17(c)、圖18(c)]一致。因初始渦核間距較大，如圖17(d)、圖18(d)所示，可明顯區分左右尾渦及渦核位置，頻譜寬度較大。

圖17 徑向風速(2018年9月9日 08：37,機型：A330-300,S*<0)

圖18 徑向風速(2018年9月19日 15：19,機型：A330-300,S*>0)

從圖17(c)可以看出，在受左側非線性側風(S*<0)條件下左右渦核的位置距離地面約37 m和42 m，距離CDL分別約517 m和529 m，渦核間距約為11 m。隨后左右尾渦均從高空下降并向外擴散，其中上風渦低于下風渦，符合數值模擬結果。

從圖18(c)可以看出，在受右側非線性側風(S*>0)條件下左右渦核的初始位置距離地面分別約20 m和21 m，距離CDL分別約254 m和288 m，渦核間距約為34 m。隨后左右尾渦均從高空下降并向外擴散，其中上風渦高于下風渦，與數值模擬結果左側非線性側風結果一致。

從圖19數值模擬結果中可以看出，尾渦對耗散開始之前，兩種不同非線性環境側風下的環量都會增加。左渦和右渦的環量在開始時是一致的。在S*<0的情況下，左渦(上風)環量在t*=0.1不斷增加，在t*=1.1之前，右渦(下風)環量大于左渦環量，隨后發生改變；CDL激光雷達探測反演結果顯示右渦環量在t*=0.8時刻之前右渦環量一直大于左渦環量，與數值模擬結果基本一致。相反，在S*>0的情況下，左渦(上風)環量大于右渦(下風)環量。

圖19 CDL探測和數值模擬結果的尾渦無量綱環量與無量綱時間的關系

通過數值模擬可以看到由于誘導速度與尾渦強度成正比，上面討論的環量差異導致兩個尾渦以不同的速率下降。尾渦對的傾角θ定義如圖4所示。從圖20可以明顯看出，當S*為負時傾斜率逐漸增加，在t*=1.0時刻后幾乎呈線性增長；同時，當S*為正時，傾斜率也逐漸增加；這種情況下，尾渦對傾斜率在t*=0.8內增長幅度相對較小，然后迅速增加并最終也呈線性增長趨勢。盡管在兩個非線性側風模擬中|S*|被設置為相同，但傾斜率的變化卻是不同的。這可能有兩個原因；首先，有其他環境因素會影響下降過程。其次，側風不會對稱地改變尾渦對的傾斜。

θ0為最大傾斜角

CDL激光雷達探測結果中，尾渦對側風和側風切變的垂直梯度(圖21)影響尾渦對的傾斜，導致左渦(上風渦)在右渦(下風渦)之間不對稱傾斜。從圖21(b)中可以看出，左渦下降速度較慢，并且比右渦保持在更高的高度。相反，從圖21(a)中可以看出，右渦高于左渦，并出現更大的橫向偏移。

H為尾渦渦心距地面高度，H0為尾渦渦心探測初始位置距地面高度

4 結論

(1)尾渦的空間分布對側風大小非常敏感，側風誘導效應導致尾渦速度場發生劇烈波動。

(2)沿高度y垂直方向的線性側風(S*=0)對尾渦對傳輸的影響較大，側風導致的誘導速度的變化引起尾渦對的衰減和運動。

(3)非線性垂直切變側風(S*≠0)的存在會造成尾渦對下沉率呈非對稱性，從而引起尾渦對傾斜并改變其橫向間隔。

(4)尾渦對的傾斜可歸因于側風垂直剖面的特征。側風的垂直二階導數影響尾渦對的傾斜。