潘衛(wèi)軍，羅昊天，羅玉明，王靖開，姜沿強

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

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

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

1 控制方程

數(shù)值模擬使用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為由于平均速度梯度產(chǎn)生的湍動能；Gω為ω的梯度；Гk和Гω分別為k和ω的有效擴散率；Yk和Yω分別為k和ω由于湍流的耗散；Sk和Sω為用戶定義源項；Gb和Gωb為浮力項。

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

前面描述的BSL模型結(jié)合了Wilcox模型和k-ω模型的優(yōu)點，但仍無法正確預測從光滑表面開始的邊界層分離和數(shù)量。其主要原因是這兩個模型都沒有考慮到湍流剪切應(yīng)力的傳輸。這導致了渦流黏度的過度預測，適當?shù)膫鬏斝袨榭梢酝ㄟ^對渦流黏度表述的限制器得到。

(6)

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

(7)

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

2 前期處理工作

2.1 尾渦模型

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

(8)

式(8)中：Γ0為飛機尾渦初始環(huán)量；MLW為飛機最大著陸重量；g為當?shù)丶铀俣龋籅為翼展；V為飛機速度。

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

(9)

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

(10)

rc0=0.052b0

(11)

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

(12)

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

無量綱時間t*可表示為

(13)

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

(14)

2.2 網(wǎng)格和初始化

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

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

表 1 A330-300機型性能和尺寸參數(shù)

表 2 A330-300尾渦參數(shù)

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

圖1 網(wǎng)格無關(guān)性收斂曲線

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

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

圖2 計算域和初始渦位置

圖3 初始切向速度剖面圖

2.3 邊界條件和計算參數(shù)

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

(15)

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

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

(16)

(17)

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

圖4為尾渦對側(cè)風切變效應(yīng)示意圖。

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

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

表3 不同類型側(cè)風的函數(shù)關(guān)系

表4 流體域環(huán)境變量

3 仿真結(jié)果

3.1 數(shù)值模擬結(jié)果

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

圖5 側(cè)風下的尾渦運動示意圖

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

(18)

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

圖6 均勻側(cè)風1 m/s

圖7 均勻側(cè)風3 m/s

圖8 均勻側(cè)風7 m/s

圖9 均勻側(cè)風10 m/s

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

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

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

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

(19)

(20)

(21)

(22)

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

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

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

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

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

圖15 不同側(cè)風類型下的渦流跡示意圖

3.2 實地探測數(shù)據(jù)對比

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

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

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

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

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

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

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

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

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

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

圖19 CDL探測和數(shù)值模擬結(jié)果的尾渦無量綱環(huán)量與無量綱時間的關(guān)系

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

θ0為最大傾斜角

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

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

4 結(jié)論

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

(2)沿高度y垂直方向的線性側(cè)風(S*=0)對尾渦對傳輸?shù)挠绊戄^大，側(cè)風導致的誘導速度的變化引起尾渦對的衰減和運動。

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

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