基于黏性渦粒子/沙粒DEM的直升機沙盲建模

譚劍鋒，何龍，于領軍，周國臣

1. 南京工業大學 機械與動力工程學院，南京 211816 2. 中國空氣動力研究與發展中心 低速空氣動力研究所 旋翼空氣動力學重點實驗室，綿陽 621000 3. 陸軍航空兵學院 航空機械工程系，北京 101123 4. 南京航空航天大學 航空學院， 南京 210016

直升機在沙漠、戈壁等環境下垂直起降和近地飛行，旋翼高速旋轉產生的下洗流將掀起細沙顆粒和碎石，形成顯著的“沙盲(Brownout)”現象，從而阻擋飛行員視線，導致直升機起降、貼地飛行困難，甚至因無法識別周圍環境(如障礙物等)導致飛行事故。美國國防部航空安全技術報告表明：“沙盲”導致的弱視環境(Degraded Visual Environments, DVEs)引發了幾乎一半的空軍直升機墜機事故，同時也是美國陸軍直升機損傷的主要因素。因此，研究直升機“沙盲”現象以及建立直升機“沙盲”現象分析方法顯得尤為重要。

為研究直升機“沙盲”現象的形成過程， Whitehouse等基于模型旋翼開展“沙盲”試驗以測試“沙盲”狀態下的可視度，然而試驗并未測量沙粒運動軌跡和流場特性。隨后Tritschler等采用高速PIV技術研究旋翼流場干擾兩相流狀態下的沙粒流動，結果表明旋翼槳尖渦的相互融合顯著增加了靠近地面的上洗速度，并影響沙粒聚集和揚起高度。為填補全尺寸直升機“沙盲”現象試驗數據空白，美國陸軍聯合研究項目實驗室的Wong等在美國陸軍Yuma試驗場開展了EH-60L的3D-LZ“沙盲”飛行試驗，獲得了全尺寸直升機“沙盲”形成過程數據。Leishman等采用單相和雙相流測試技術研究旋翼“沙盲”狀態下的沙粒起跳問題，結果表明受地面影響，旋翼槳尖渦發生拉伸、融合，并形成高度三維非定常特性，誘發地面渦的形成，從而導致沙粒起跳、懸浮。為進一步明確旋翼渦系的影響，2017年，Leishman等采用高速流體顯示技術(Flow Visualization，FV)和高速兩分量PIV研究旋翼地效下的旋翼渦系結構，研究表明受地面影響，旋翼槳葉內段渦系與槳尖渦存在速度差異，并在靠近地面融合，形成復雜的湍流結構，影響地面射流特性，最終影響“沙盲”形成。上述試驗表明直升機“沙盲”現象受旋翼與地面非定常氣動干擾影響顯著。

為研究直升機“沙盲”現象的形成過程，眾多研究者開展直升機“沙盲”現象的數值分析方法研究。Quackenbush等采用基于自由尾跡模型的CHARM軟件研究旋翼靠近地面流場，并研究直升機“沙盲”現象，結果表明直升機貼地飛行狀態下地面渦引發沙粒移動。隨后，Whitehouse等基于自由尾跡Free-wake和面元法建立了旋翼/地面氣動干擾模型，研究旋翼靠近地面的流場特性，并耦合沙粒運動模型和可視化模型模擬直升機“沙盲”現象。結果表明，此方法可實現直升機“沙盲”現象模擬，但并未考慮旋翼槳尖渦的黏性效應和地面黏性效應，也未考慮沙粒的摩擦碰撞特性。為此，Andrea采用基于面元-時間步進自由尾跡模型的ADPANEL軟件研究旋翼靠近地面的非定常流動現象，并耦合粒子輸運模型(Particle Transport Model，PTM)研究EH101和縱列式直升機“沙盲”特性，但方法中仍未考慮黏性效應、地面復雜流動以及沙粒起跳真實特性。此外，Keller等基于渦量輸運模型(Vorticity Transport Model，VTM)和無黏地面模型研究旋翼靠近地面的復雜流動現象，并耦合基于沙粒起跳通量和起跳速度罰值的半經驗模型構建沙粒輸運模型PTM，從而模擬直升機“沙盲”現象，但仍未考慮沙粒間以及與地面的碰撞。為此，Leishman等采用時間精確自由尾跡模型和地面鏡像模型模擬旋翼近地面流場，耦合基于沙粒輸運方程的沙粒動力學模型，建立拉格朗日沙粒跟蹤方法。方法采用Bagnold沙粒揚起速度罰值模型和沙粒摩擦速度罰值模型，并嵌入Shao沙粒夾帶和起跳概率模型，模擬直升機“沙盲”現象。結果表明此方法可實現“沙盲”現象模擬，且旋翼與地面干擾的非定常尾跡結構對沙云特性影響顯著。然而，方法仍未考慮旋翼尾跡和地面的無黏效應，且采用基于半經驗的沙粒夾帶模型、沙粒起跳模型、以及速度罰值的沙粒通量模型計算沙粒揚起過程，而未考慮沙粒碰撞存在的能量交換。為此，2020年，徐國華等采用CCM++軟件和離散動力學EDEM軟件開展了全尺寸直升機“沙盲”現象初步研究，但研究的沙粒數量較少(萬級)，沙粒間距較大，因此在“沙盲”形成過程和計算方法方面仍需進一步研究。

基于拉格朗日體系的黏性渦粒子方法(Viscous Vortex Particle Method, VVPM)能較好捕捉旋翼/地面干擾效應和復雜流場特性，本文將采用非定常面元/VVPM/黏性地面氣動模型模擬旋翼時變尾跡結構和旋翼/地面非定常氣動干擾效應，并基于離散動力學，建立沙粒DEM模型和沙粒-流場耦合模型，構建基于黏性渦粒子/沙粒DEM(Discrete Element Method)的直升機“沙盲”現象分析方法。通過與美國陸軍Yuma試驗場的EH-60L“沙盲”飛行試驗測量結果和基于沙粒夾帶模型、沙粒起跳模型、速度罰值沙粒通量模型的拉格朗日沙粒跟蹤方法計算結果對比，驗證本文方法計算“沙盲”的準確性。隨后，研究直升機懸停和前飛狀態下的“沙盲”形成過程和特性。

1 計算方法

1.1 旋翼/地面非定常氣動干擾模型

直升機“沙盲”狀態的旋翼/地面非定常干擾流場受旋翼尾跡和旋翼氣動載荷影響顯著。非定常面元法能較好模擬旋翼懸停和前飛狀態的槳葉氣動載荷，且相比于CFD具有較高的計算效率，已應用于單旋翼、雙旋翼氣動干擾分析。為此，本文采用非定常面元法模擬直升機旋翼槳葉氣動特性，直升機旋翼流場勢函數可表示為

(1)

式中：為槳葉表面匯；為槳葉表面偶極子；為槳葉表面；為尾跡渦面；為表面法矢量；為空間點位置矢量。

直升機旋翼槳葉邊界條件表示為

(2)

(3)

根據Neumann邊界條件和后緣Kutta條件將式(2)轉為代數方程，計算得到旋翼槳葉匯和偶極子分布。根據非定常伯努利方程計算得到旋翼槳葉壓力系數分布：

(4)

式中：為參考壓強；為空氣密度；為當地速度；為當地壓強；為當地參考速度。

根據三維Navier-Stokes方程，采用黏性渦粒子法模擬旋翼尾跡渦量分布：

(5)

式中：為運動黏性系數；為旋翼尾跡渦矢量分布；為旋翼速度場。

為體現地面對旋翼尾跡結構和流場的干擾，采用基于第二類Fredholm方程的黏性地面氣動模型：

(6)

式中：為地面切矢量；為直升機旋翼槳葉和尾跡誘導的滑移速度，即

(7)

采用四邊形單元離散地面，根據式(6)黏性地面氣動模型和離散非均勻Neumann問題的Green積分函數，耦合渦擴散方程計算得到地面渦量擴散到流體的渦量為

(8)

式中：為三維熱擴散核函數；為積分面積。

耦合式(1)、式(3)、式(5)～式(7)實現旋翼槳葉氣動模型、尾跡模型與黏性地面氣動模型的相互誘導耦合，耦合式(5)和式(8)實現旋翼尾跡和黏性地面氣動模型的渦耦合，形成旋翼/地面非定常氣動干擾模型，實現對直升機“沙盲”狀態下的旋翼/地面非定常干擾流場求解。

1.2 沙粒DEM模型

“沙盲”的形成過程是顆粒相為稀相兩相流動過程，在拉格朗日體系下，追蹤顆粒相的運動，從而得到沙云的形成特征。對于沙粒間的相互作用，離散單元法(Discrete Element Model)在離散度大、流動顯著的顆粒流系統分析中得到廣泛應用。基于離散動力學理論，沙粒運動方程由牛頓第二定律得出,并考慮沙粒受到的流體力、碰撞力、重力以及沙粒旋轉：

(9)

(10)

式中：為沙粒顆粒質量；為沙粒速度矢量；為沙粒慣量；為沙粒角速度矢量；為沙粒力矩矢量。通過采用顯式積分技術，計算沙粒的運動軌跡，得到沙粒在空間的分布形態。

(11)

(12)

(13)

式中：為滾動摩擦系數；為接觸點到沙粒質量中心的距離；為沙粒角速度單位矢量。

圖1 沙粒接觸碰撞模型Fig.1 Contact model of sand particle

1.3 沙粒-流場耦合模型

基于沙粒-流場耦合模型，將沙粒DEM動力學模型嵌入直升機旋翼/地面非定常干擾流場，實現直升機“沙盲”現象分析。由于直升機“沙盲”的沙粒直徑為100 μm量級，因此忽略沙粒對旋翼/地面非定常干擾流場的影響，但需體現旋翼/地面非定常干擾流場對沙粒動力學特性的影響。為此，將流場作用于沙粒的氣動力耦合至沙粒動力學方程，從而實現沙粒-流場耦合。

沙粒與流場的氣-固兩相間作用主要通過相間的作用力來耦合，沙粒受到流體力為

(14)

式中：為空氣密度；為流場速度；為沙粒移動速度；為沙粒等效面積。

根據低速流體的Stokes定律，計算得到沙粒在空氣中的阻力系數為

(15)

式中：沙粒雷諾數為

(16)

其中：為沙粒等效直徑。

旋翼/地面干擾流場與沙粒DEM耦合計算流程如圖2所示。圖中：為旋翼/地面干擾流場計算時間；Δ為計算步長；為計算終止時間；Δ為沙粒DEM計算時間步長：

基于上述沙粒-流場耦合模型，將直升機旋翼/地面非定常氣動干擾分析方法與沙粒DEM模型耦合，建立直升機“沙盲”現象數值分析方法。

圖2 直升機“沙盲”計算流程圖Fig.2 Computer flow diagram of helicopter brownout

2 直升機“沙盲”特性

2.1 旋翼/地面干擾流場驗證

為了驗證本文直升機“沙盲”現象數值分析方法計算旋翼/地面干擾流場的準確性，采用美國馬里蘭大學旋翼/地面干擾試驗，模型旋翼的槳葉數目、半徑、扭轉角、弦長、槳尖馬赫數、翼型分別為2、0.086 m、0°、0.019 6 m、0.08、圓弧翼型。旋翼總距和離地高度分別為12°、1.5。

圖3給出了旋翼不同徑向位置處的徑向速度隨高度變化，以及試驗結果和CFD的OVERTUNS軟件計算結果(其中為徑向速度，為垂向速度)。各徑向位置處，計算得到的速度隨高度變化趨勢與試驗測量結果和CFD計算結果基本吻合。隨著高度增加，徑向速度先快速增加后減小，形成典型的地面射流現象。此外，隨著徑向位置增加，徑向速度峰值對應的高度逐漸減小，速度峰值先增加后略有減小，且速度峰值和其對應的高度與試驗結果均吻合較好。因此，本文直升機“沙盲”現象數值分析方法較好模擬旋翼/地面干擾下的流場特性。

圖3 旋翼/地面干擾下的徑向速度變化特性Fig.3 Varations of radial velocities during rotor-ground interaction

2.2 前飛狀態EH-60L“沙盲”驗證

為驗證本文直升機“沙盲”現象數值分析方法的準確性，采用Wong和Tanner在美國陸軍 Yuma試驗場YPG的Oasis著陸區EH-60L直升機“沙盲”飛行試驗。飛行測試試驗區如圖4所示，著陸區包括7條寬度不一的跑道，其中飛行試驗在跑道1、3、5、7中進行，其他跑道設置了包括電線桿、電纜、小型建筑、汽車等障礙物，以便測量EH-60L機載傳感器靈敏度。試驗通過由多臺相機構成的立體拍攝系統記錄EH-60L飛行時“沙盲”照片，并獲得“沙云”位置和形態。

圖4 “沙盲”試驗的Oasis著陸區示意圖Fig.4 Schematic of Oasis landing zone used for brownout experiments

根據YPG試驗場Oasis著陸區沙粒樣本分析獲得沙粒主要分布在0～1 000 μm。由于大直徑沙粒重量大，移動和空氣中懸浮困難，而占主要成分的小直徑沙粒重量輕，移動和懸浮容易，因此，本文分析中建立長度為160 m、寬度為160 m的地面模型，采用球形沙粒，直徑為1～100 μm，且沙粒直徑和位置隨機分布。

EH-60L直升機旋翼由4片后掠20°槳葉構成，槳葉半徑、弦長、翼型、后掠角分別為8.18 m、0.527 4 m、SC1095和SC1095R8、20°，槳葉采用非線性負扭,如表1所示。EH-60L直升機平均前飛速度為37 km/h，飛行重量為7 394 kg，飛行軌跡如圖5和圖6所示。

表1 EH-60L直升機參數Table 1 Parameters of EH-60L helicopter

圖5 各時間點直升機空間位置Fig.5 Positions of helicopter at different time points

圖6 直升機飛行軌跡和Oasis著陸區模型Fig.6 Trajectory of helicopter and model of Oasis landing zone

計算得到此飛行狀態下各時間點直升機沙云位置、飛行測試結果、以及基于拉格朗日沙粒跟蹤方法的計算結果如圖7～圖9所示。

圖7 沙云形態t=0 sFig.7 Shape of dust clouds at t=0 s

本文方法計算得到各時間“沙云”俯視圖和側視圖位置均與飛行測試的“沙云”結果比較吻合，且計算獲得的“沙云”輪廓與飛行測試輪廓一致。由此表明本文方法能較好的模擬直升機近地飛行的“沙云”輪廓和形成過程。

在=0,2 s，本文方法和拉格朗日沙粒跟蹤方法計算得到“沙云”前端位置基本相同，均略超前于飛行測試結果(見圖8(a))。=4 s時，拉格朗日沙粒跟蹤方法計算得到“沙云”前端位置明顯超前于飛行測試結果，而本文方法計算得到的“沙云”前端位置與飛行測試結果一致(見圖9(a))。此外，拉格朗日沙粒跟蹤方法計算得到的“沙云”向后收縮，與飛行測試結果的“沙云”擴張形式存在明顯差異(見圖9(a))，而本文方法計算得到的“沙云”向兩側擴散，與飛行測試結果一致。因此，相比于拉格朗日沙粒跟蹤方法，本文方法能更好捕捉“沙云”前端位置和形狀變化特性。

圖8 沙云形態t=2 sFig.8 Shape of dust clouds at t=2 s

圖9 沙云形態t=4 sFig.9 Shape of dust clouds at t=4 s

在=0,2,4 s，拉格朗日沙粒跟蹤方法計算的“沙云”高度低于飛行測試結果 (見圖7(b)、圖8(b)、圖9(b))，尤其在=10～16，而本文計算方法計算的“沙云”高度與飛行測試結果更接近，“沙云”輪廓與飛行測試結果基本一致。此外，隨著“沙云”向后移動，拉格朗日沙粒跟蹤方法計算得到的“沙云”高度先增加，后快速降低，由此表明直升機下方的沙粒揚起，后快速下沉，而本文方法計算得到的沙粒先揚起，而后懸浮，并與飛行測試的沙云形態特性一致，且“沙云”高度變化特性與飛行測試結果更吻合。

因此，相比基于自由尾跡模型和沙粒起跳罰值等模型的拉格朗日沙粒跟蹤方法，本文直升機“沙盲”數值分析方法計算得到的直升機“沙云”現象精度更高，與美國陸軍YPG試驗場Oasis著陸區EH-60L“沙盲”飛行試驗結果更吻合。

2.3 EH-60L直升機懸停“沙盲”形成特性

基于上述EH-60L直升機參數和地面參數，計算模擬EH-60L懸停狀態(=1.5)“沙盲”形成過程，結果如圖10所示。

=0.5 s時旋翼/地面干擾形成的地面射流對沙粒移動影響時間較短，暫未形成比較明顯的“沙云”現象。受到旋翼槳尖渦與地面干擾的非定常流場作用(見圖11，其中為垂直方向速度)，沙粒向外側移動，在旋翼槳尖渦處形成明顯的空白區，且在旋翼槳尖渦之間形成堆積(見圖10(a))，主要原因在于旋翼槳尖渦正下方速度較大，沙粒移動速度較快，而槳尖渦之間的沙粒受到流場作用起跳，并受沙粒間碰撞作用而停止移動形成堆積。

圖10 懸停狀態EH-60L“沙盲”形成過程Fig.10 Process of brownout for EH-60L in hover

=1.0 s時，旋翼地面射流對沙粒的作用時間增長，靠近旋翼的沙粒堆積仍然存在，而受到最外側旋翼槳尖渦形成的地面射流作用，沙粒開始起跳、揚起。時間增加到2.0 s，沙粒起跳高度逐漸增加，揚起現象更加明顯。時間=4.0 s時，受地面射流和槳尖渦作用，細小顆粒揚起高度進一步增加，并沿著旋翼速度場移動。受重力和流體力的作用，沙粒向旋翼上方和內側移動，并逐步形成沙云。時間增加到7.0 s后，沙粒起跳更加明顯，揚起的沙粒形成非常明顯的沙云，并逐步包圍EH-60L直升機，從而形成典型的“沙盲”現象。

圖11 懸停狀態EH-60L/地面干擾流場Fig.11 Flowfield of EH-60L/ground interaction in hover

2.4 EH-60L直升機前飛狀態“沙盲”形成特性

基于上述EH-60L直升機參數和地面參數，計算模擬EH-60L前飛狀態(=1.0)沙盲形成過程，結果如圖12所示。

直升機旋翼槳尖渦與地面干擾產生的兩個卷起槳尖渦，推動直升機兩側沙粒向外移動，而中間的沙粒向后移動，前部沙粒向前移動(見圖12(a)和圖12(b))。時間=0.5 s，直升機兩側和前方的流場作用時間較短，沙粒暫未揚起。時間增加到1.0 s，受到直升機兩側地面卷起槳尖渦的作用，沙粒向外移動速度更明顯，且受卷起槳尖渦向上誘導作用(見圖13(a)，其中為平面合速度)，沙粒出現起跳、碰撞、輕微揚起。時間增加到2.0 s時，直升機兩側地面卷起槳尖渦所產生的沙粒起跳、揚起現象更為明顯。同時，受直升機前方卷起槳尖渦誘導的向前地面射流作用(見圖13(b)，其中為平面合速度)，機頭前部的沙粒向前移動和起跳。當時間增加到5.0 s，受兩側卷起槳尖渦的顯著影響，沙粒揚起現象顯著，形成兩側沙云現象，直升機前部的沙粒揚起現象也更為明顯。當時間增加7.0～10.0 s，直升機兩側的沙粒揚起速度更快，在重力和流體力作用下，沙粒懸浮并逐步回落，形成典型的回流現象。受直升機前方卷起槳尖渦作用，直升機機頭前部的沙粒揚起較為嚴重，出現懸浮的沙云現象，最終形成典型前飛“沙盲”現象。

圖13 前飛狀態EH-60L/地面干擾流場Fig.13 Flowfield of EH-60L/ground interaction in forward flight

3 結 論

耦合旋翼黏性渦粒子和黏性地面氣動模型，建立基于離散動力學的沙粒DEM模型和沙粒-流場耦合模型，構建直升機“沙盲”現象分析方法，并應用于美國陸軍Yuma試驗場EH-60L“沙盲”計算，得到如下結論：

1) 建立的直升機“沙盲”現象分析方法較準確模擬EH-60L全尺寸直升機近地飛行狀態下的“沙云”輪廓。

2) 本文方法計算得到各時間的“沙云”形狀、揚起位置、高度與飛行測試結果基本一致，并能較準確體現各飛行狀態下的“沙云”形成過程。

3) 相比基于自由尾跡模型和沙粒起跳罰值等模型的拉格朗日沙粒跟蹤方法，本文建立的直升機“沙盲”現象分析方法計算得到的“沙云”輪廓精度更高，與飛行測試結果更吻合。

4) 懸停狀態下的旋翼槳尖渦與地面干擾形成地面射流，推動沙粒移動堆積，并導致沙粒碰撞起跳，在旋翼/地面干擾流場作用下揚起，形成沙云，覆蓋直升機前方視野，形成沙盲現象。

5) 前飛狀態下的旋翼槳尖渦與地面干擾形成兩側擴散卷起槳尖渦和機頭前方卷起渦，推動沙粒移動、碰撞起跳，并在槳尖渦誘導作用下揚起，覆蓋直升機前方視野，形成“沙盲”現象。

后續研究可將旋翼CFD方法和“沙盲”求解方法耦合，更好體現現代旋翼槳葉先進翼型對“沙盲”現象的影響，并可進一步研究旋翼構型與氣動設計對“沙盲”現象的影響規律。