基于CFD耦合DEM并行算法的旋翼飛行器槳葉沙盲現(xiàn)象模擬分析

文韜，程黨根

（1.長沙航空職業(yè)技術(shù)學(xué)院航空機電設(shè)備維修學(xué)院，長沙410124）

（2.空軍航空維修技術(shù)學(xué)院 湖南省飛機維修工程技術(shù)研究中心，長沙410124）

0 引 言

旋翼飛行器在沙漠地區(qū)、雪地等具有離散相（Dispersed Phase）顆粒物的特殊環(huán)境中地效懸停IGE（In Ground Effect）或低速前飛時，會由旋翼下洗流帶動并卷起大量沙塵或雪花，造成嚴(yán)重遮擋飛行員視線的現(xiàn)象被稱為“沙盲（Brownout）”［1］。沙盲會極大地威脅飛行安全，還會對旋翼槳葉、發(fā)動機或機身本體造成嚴(yán)重的機械磨損［2-5］。

目前，研究旋翼飛行器沙盲現(xiàn)象采用的實驗方法主要有粒子圖像測速法（Particle Image Velocimetry，簡稱PIV）、半經(jīng)驗公式法、渦量輸運模型法、黏性渦粒子法等。T.E.Lee等［6］基于PIV方法，在室內(nèi)環(huán)境條件下，對縮比旋翼進(jìn)行實驗觀察，發(fā)現(xiàn)了氣動影響下近地附近沙塵擴散的過程、湍流的形成；B.T.Hance［7］對增加了機身的全機縮比模型進(jìn)行沙盲實驗分析，研究在沙盲環(huán)境下機身對旋翼槳尖渦所產(chǎn)生的影響；胡健平等［8］通過實驗觀察縮比旋翼模型在各種懸停高度時渦的狀態(tài)，測量結(jié)果顯示懸停高度在中間值時，槳尖渦的伸展現(xiàn)象明顯超過耗散現(xiàn)象，但在增大懸停高度時，這種趨勢則發(fā)生了逆轉(zhuǎn)，同時采用雷諾平均方程和剪切應(yīng)力輸運湍流數(shù)值模型，通過碰撞接觸模型的離散元模型，基于“多球法”構(gòu)建了用于直升機地效飛行沙盲現(xiàn)象模擬的CFD-DEM耦合計算方法，對直升機在多種高度懸停下產(chǎn)生沙盲現(xiàn)象時沙塵顆粒的受力狀態(tài)進(jìn)行了分析，并對不同時間的沙塵運動形態(tài)和規(guī)律進(jìn)行了計算，其研究結(jié)果表明，流場中大部分沙塵顆粒是不能形成沙塵云的，且外層空間的沙塵濃度比內(nèi)層高，槳盤平面的上、下層區(qū)域的沙塵顆粒運動方向、切向速度各異，并與可得到的試驗結(jié)果進(jìn)行對比分析，驗證了數(shù)值計算方法的有效性；朱明勇等［9］基于非結(jié)構(gòu)網(wǎng)格技術(shù)和動量源模型的旋翼計算流體力學(xué)（CFD）方法，模擬貼地飛行時直升機旋翼非定常氣動特性，且考慮了實際飛行環(huán)境下旋翼的運動操縱和旋翼配平特性，提出并建立了基于遺傳算法和擬牛頓法的高效混合迭代算法，計算了特殊環(huán)境中地效懸停作用下的槳葉拉力增益、功率變化，有效解決了渦流理論方法較難模擬的“小速度前飛旋翼需用功率突增”問題。雖然上述方法能夠較為精確地給出旋翼飛行器沙盲時沙塵的運動細(xì)節(jié)，但卻很難采用全尺寸機身、槳葉進(jìn)行實驗，并且無法對多參數(shù)影響進(jìn)行分析，相比之下，數(shù)值模擬則容易的多。然而，大多數(shù)數(shù)值計算方法都存在一些非常難以解決的問題，主要集中在以下三點：

（1）槳尖渦需要經(jīng)過很長的時間才能傳至地面，而在這個過程中渦極易耗散掉，而要減小耗散會極大地增加計算開銷，尤其當(dāng)旋翼距離地面較高時。

（2）要解決地面邊界層以及邊界處的湍流問題是極其困難的，尤其是采用非CFD方法。

（3）要很好地模擬流場和顆粒的相互作用也是具有極高難度的，顆粒的受力模型非常復(fù)雜，這也是目前的研究熱點之一。

而且，當(dāng)使用全尺寸的氣動計算模型和顆粒形狀各異、粒徑分布廣泛、受力復(fù)雜的沙塵顆粒進(jìn)行耦合計算時，會造成巨大的計算資源消耗。

本文假設(shè)顆粒體積分?jǐn)?shù)遠(yuǎn)小于流場計算域的10%，并且沙塵顆粒全部采用球體模型進(jìn)行簡化以控制計算開銷。采用基于雷諾平均N-S方程及k-ω（SST）湍流模型，結(jié)合計算流體力學(xué)（Computational Fluid Dynamics，簡稱CFD）與離散元（Discrete Element Method，簡稱DEM）耦合并行算法，使用Fluent軟件進(jìn)行流場連續(xù)相計算沙塵的完整運動軌跡和顆粒分布。

1 流場計算方法

在沙塵顆粒所有的受力中，流場對其作用力是最為顯著的，因而流場計算的準(zhǔn)確性決定了沙塵顆粒運動的準(zhǔn)確性。目前，研究兩相流的觀點主要有兩種，一種觀點認(rèn)為流體是連續(xù)介質(zhì)，顆粒群為離散體系；另一種觀點認(rèn)為流體是連續(xù)介質(zhì)，顆粒群是擬連續(xù)介質(zhì)或擬流體［10］。在有兩種聚合物的混合體系中，離散相是被包圍的聚合物，其對混合物機械性能的影響較小。連續(xù)相是對黏附、滲透、擴散、導(dǎo)熱等性能參數(shù)影響較大的另一種聚合物。本文引入以變形后瞬時坐標(biāo)為自變量的歐拉（Eulerian）坐標(biāo)系，流場是建立在歐拉模型下的連續(xù)相。

1.1 物理模型

為了平衡計算精度和計算開銷，采用的槳葉為縮比的單片NACA 0012翼型矩形槳葉，弦長0.1 m，根切16.67%，距離地面高度為0.3 m，半徑為0.6 m，轉(zhuǎn)速為200 rad/s，槳距為12°。

設(shè)置沙塵顆粒于槳盤正上方0.4 m處固定產(chǎn)生，以模擬槳盤吸入沙塵后隨流場到達(dá)地面附近重新被流場卷起的過程。沙塵顆粒的物理特性及受力狀態(tài)對于其在流場中的運動狀態(tài)至關(guān)重要，作用于單個沙塵顆粒上的作用力非常復(fù)雜。從理論上分析，流場中運動的顆粒主要受到曳力、重力、Saffman升力、Magnus力、黏結(jié)力、附加質(zhì)量力、Basset力、浮力等影響。通常，黏結(jié)力、浮力、Basset力和附加質(zhì)量力由于相較其他力都非常小可忽略不計，而對于本文中顆粒體積分?jǐn)?shù)遠(yuǎn)小于10%，從而可以忽略其對流場動量的影響，并且沙塵顆粒粒徑微小，從而只考慮曳力和重力對顆粒移動速度和運動軌跡的影響。曳力采用不考慮孔隙率的Free steam模型［11］，計算公式如下：

式中：A為顆粒的投影面積；Cd為曳力系數(shù)。

式中：Re為雷諾數(shù)。

Saffman升力［12］是在有速度梯度的流場中顆粒兩側(cè)流場速度不一致而由低速側(cè)指向高速側(cè)的力：

式中：B為實驗常數(shù)；ρ1為氣相密度；η為氣體黏性；Rd為顆粒半徑；v1為氣相速度；v2為離散相速度；?為哈密頓算子。

Magnus力是由粒子自身旋轉(zhuǎn)而產(chǎn)生的力，其表達(dá)式為

式中：A為經(jīng)驗常數(shù)。

對于沙塵或者粉塵這類顆粒極其微小的離散相，顆粒與顆粒之間的接觸碰撞模型適合使用無滑移Hertz-Mindlin with JKR模型。為了節(jié)約計算資源而將其簡化為球體模型，沙塵顆粒的物理參數(shù)設(shè)置如表1所示，其中，恢復(fù)系數(shù)用碰撞前法向相對速度與碰撞后法向相對速度的比值表示；滾轉(zhuǎn)系數(shù)用滾動摩擦力矩與支持力的比值表示；靜摩擦系數(shù)用最大靜摩擦力與接觸面間面積的比值表示［6］。

表1 沙塵顆粒的物理參數(shù)Table 1 Physical parameters of dust particles

沙塵顆粒初始狀態(tài)的分布情況如圖1所示，沙塵顆粒分布域為1.8 m2的正方形區(qū)域，平鋪厚度為一層，顆粒數(shù)量為30 000顆，且為均勻分布、靜止生成的球形顆粒。

圖1 沙塵顆粒初始狀態(tài)的分布情況Fig.1 Distribution of initial state of dust particles

采用ANSYS-Fluent軟件計算連續(xù)相的懸停流場，采用EDEM軟件計算離散相的灰塵顆粒，動量信息采用C++語言編譯的API（Application Programming Interface）耦合接口進(jìn)行數(shù)據(jù)傳遞，使用能從機理上反映連續(xù)相和離散相之間相互作用，并能真實反映出耦合特征的單向耦合模型Lagrange-Euler進(jìn)行并行耦合計算。

1.2 網(wǎng)格系統(tǒng)

本文采用SMM（Sliding Mesh Method）技術(shù)，將計算域劃分為槳葉附近的旋轉(zhuǎn)區(qū)域和地效流場的靜止區(qū)域進(jìn)行瞬態(tài)求解，計算域之間通過Interface進(jìn)行數(shù)據(jù)交換，旋轉(zhuǎn)區(qū)域網(wǎng)格對于捕捉流場細(xì)節(jié)至關(guān)重要，弦長雷諾數(shù)為0.5×106，y+≈1，計算出第一層網(wǎng)格厚度約為4×10-6m。為了在保證計算精度，減小截斷誤差的同時控制計算開銷，網(wǎng)格附面層增長率設(shè)置為1.2，并且采用槳葉表面網(wǎng)格以四邊形為主和三角形混合的技術(shù)，在附面層使用各向異性的六面體網(wǎng)格具有比三棱柱網(wǎng)格更好的精度，而且還能在曲率較大的曲面很好地控制網(wǎng)格數(shù)量。在外層，六面體網(wǎng)格和四面體網(wǎng)格通過金字塔形網(wǎng)格進(jìn)行連接。槳葉表面以及旋轉(zhuǎn)區(qū)域網(wǎng)格如圖2所示，計算域網(wǎng)格的剖面形態(tài)如圖3所示。

圖2 槳葉網(wǎng)格模型Fig.2 Blade mesh and rotating area volume mesh

圖3 計算域剖面網(wǎng)格［13］Fig.3 Calculation domain section grid［13］

靜止區(qū)域的地面邊界第一層厚度同樣設(shè)置為4×10-6m，設(shè)置增長率為1.2，并且對槳盤直徑長度附近的垂直地面方向的背景網(wǎng)格進(jìn)行加密，以平衡計算精度和網(wǎng)格數(shù)量。

1.3 數(shù)值方法

DEM方法能夠?qū)μ幱诹鲌鲋械拿恳涣Ｉ硥m顆粒的受力建立數(shù)學(xué)模型，準(zhǔn)確地給出每一粒沙塵在每個時間步的運動狀態(tài)以及受力情況，因而每一粒沙塵受到流場影響的空間擴散方式、運動軌跡等細(xì)節(jié)都能夠詳細(xì)呈現(xiàn)。

在CFD耦合DEM數(shù)值計算中，采用基于雷諾平均N-S方程的壓力基求解器作為計算流體力學(xué)（CFD）軟件的求解器［14］，使用Phase coupled SIMPLE算法求解，湍流模型選用能更好捕捉壁面速度梯度變化的k-ω（SST）兩方程模型［15］，位于槳盤上部2R位置的頂部邊界設(shè)置為壓力進(jìn)口，圓柱形遠(yuǎn)場側(cè)面邊界條件設(shè)置為壓力出口，所有壁面條件設(shè)置為無滑移，對流相采用二階迎風(fēng)格式，擴散相采用中心差分格式，非穩(wěn)態(tài)相采用二階隱式格式求解［3］。

2 懸停流場與沙塵運動規(guī)律分析

2.1 懸停流場分析

旋翼飛行器利用地面效應(yīng)懸停IGE狀態(tài)具有和非利用地面效應(yīng)懸停OGE（Out Ground Effect）完全不一樣的流場，如圖4所示。槳尖處的誘導(dǎo)速度在近地面處受到地面的阻礙，從而形成極大的滯止渦。康寧等［16］采用N-S方程對近地飛行的旋翼流場進(jìn)行了數(shù)值模擬，分析了槳盤處誘導(dǎo)速度隨各參數(shù)的變化規(guī)律，結(jié)果表明，在地面渦出現(xiàn)的飛行高度和速度范圍內(nèi)，槳盤前緣處下洗速度顯著增加。在OGE狀態(tài)下，槳根處流場直線向下，而在IGE狀態(tài)下，流場受到地面阻擋形成高壓區(qū)，流場方向與OGE狀態(tài)相反，這與R.Sambaraju等［17］采用PIV的測試結(jié)果完全符合。

圖4 計算域剖面速度場Fig.4 Velocity field of calculation domain section

2.2 沙塵顆粒運動軌跡的分布規(guī)律

流場穩(wěn)定之后的近地面速度場是影響顆粒運動速度與軌跡的最重要因素，影響著顆粒在空間的分布狀態(tài)以及運動軌跡。穩(wěn)定后的流場速度分布如圖5（a）所示，顆粒受流場影響剛到達(dá)地面的空間分布情況分別如圖5（b）和圖5（c）所示，此時的流場并未完全穩(wěn)定，并且顆粒主要受到重力和誘導(dǎo)速度的影響于0.09 s到達(dá)地面附近，但是顆粒已初現(xiàn)受流場影響而呈現(xiàn)的螺旋型分布。

圖5 近地流場速度分布以及沙塵顆粒剛接觸地面后的分布Fig.5 Velocity distribution of near ground flow field and distribution of dust particles just after contacting the ground

隨著流場進(jìn)一步發(fā)展以及沙塵顆粒受到流場持續(xù)的影響，在空間呈現(xiàn)新的分布狀況。0.17 s時沙塵空間分布狀況如圖6所示，可以看出：此時的沙塵已經(jīng)開始受到槳根和槳尖處流場的影響，由于曳力比重力大而開始向空中卷起；槳根與槳尖之間由于并沒有法向與地面的速度而部分開始向槳根中心移動，部分向槳尖處移動，這和流場計算的結(jié)果完全吻合。圖6（c）定量地給出了中心—半徑方向的顆粒數(shù)量分布，大量的顆粒分布集中于距離旋轉(zhuǎn)中心1 400～1 900 mm的半徑上，且槳根處也聚集了較高濃度的顆粒。

圖6 0.17 s時沙塵顆粒分布Fig.6 Dust particle distribution at 0.17 s

當(dāng)流場繼續(xù)發(fā)展至0.3 s時，沙塵顆?？臻g分布狀況如圖7所示，可以看出：槳根附近的沙塵顆粒已被流場卷起到了相當(dāng)高度，如果在試驗時飛行器有機身，則必然在機身貼壁面形成沙盲現(xiàn)象；同時，在大約2.5R之外處會有從地面卷起的揚塵，從顆粒數(shù)量統(tǒng)計可以看到，槳根處和槳尖處會將其間的沙塵向兩端帶動，在流場作用下沙塵顆粒被卷起形成沙盲現(xiàn)象。

3 結(jié) 論

（1）本文所提方法能夠?qū)π盹w行器利用地面效應(yīng)懸停流場進(jìn)行精確模擬，模擬結(jié)果與用PIV測試的結(jié)果完全符合。

（2）在對懸停流場進(jìn)行模擬的同時，采用離散元計算方法能夠?qū)κ芰鲌龀掷m(xù)影響的沙塵顆粒的空間分布狀態(tài)及受力運動軌跡進(jìn)行較好地模擬。

（3）相比較其他實驗方法和半經(jīng)驗方法，本文方法能從形成機理上，多參數(shù)化地了解地效流場發(fā)展規(guī)律、細(xì)微沙塵顆粒受力運動和沙盲形成因果。

本文仍存在不足及需要改進(jìn)的地方，例如流場遠(yuǎn)場計算域給定為4R，而顆粒運動計算域給定為正方形邊長3R，從而導(dǎo)致很多粒子并沒有形成揚塵就離開了計算域。而有些粒子會在壁面處沉積，這不符合物理現(xiàn)象，對于沙盲形成的整個過程的沙塵顆粒捕捉是不足的。其次，地效時的沙塵數(shù)量及其龐大，本文只給出一層沙塵設(shè)置在槳盤上方，這對于從地面卷起的現(xiàn)象模擬并不完全具有代表性，而對地面附面層的處理具有極高的難度。