陳希， 招啟軍

南京航空航天大學(xué) 直升機旋翼動力學(xué)國家級重點實驗室， 南京 210016

考慮遮蔽區(qū)影響的旋翼三維水滴撞擊特性計算新方法

陳希， 招啟軍*

南京航空航天大學(xué) 直升機旋翼動力學(xué)國家級重點實驗室， 南京 210016

針對直升機旋翼三維黏性流場特有的復(fù)雜環(huán)境，建立了一種基于歐拉法的旋翼三維水滴撞擊特性計算的新方法。首先，在旋翼槳葉嵌套網(wǎng)格的基礎(chǔ)上，發(fā)展了一套用于預(yù)測旋翼繞流流場的計算流體力學(xué)(CFD)模擬方法。然后，為克服傳統(tǒng)直升機旋翼二維水滴撞擊特性計算方法的不足，充分考慮旋翼流場的三維效應(yīng)，在嵌套網(wǎng)格中基于歐拉法求解旋翼三維水滴撞擊流場。其中，為解決尾流等區(qū)域的密度脈沖現(xiàn)象所引起的穩(wěn)定性和收斂性問題，提出并建立了遮蔽區(qū)擴散模型。該模型通過判斷遮蔽區(qū)變量，在計算域中動態(tài)生成遮蔽區(qū)域，并隨迭代步數(shù)逐漸擴散。最后，通過與NACA0012翼型及國外UH-1H槳葉的試驗和計算結(jié)果的對比，驗證了旋翼三維水滴撞擊特性計算新方法的可靠性，并進行了溫度和水滴當量直徑(MVD)對旋翼三維水滴撞擊特性的影響分析。結(jié)果表明：遮蔽區(qū)擴散模型的加入，使二維情況的計算時間減少了22%，并增加了三維情況的計算穩(wěn)定性，顯著提高了旋翼三維水滴撞擊特性的計算效率；沿著旋翼槳葉展向位置增大的方向，旋翼槳葉剖面水滴撞擊范圍有所增大，最大水滴局部收集系數(shù)呈先增加后減少再增加的變化趨勢，其變化幅度接近50%；旋翼槳葉表面的水滴撞擊區(qū)域和水滴局部收集系數(shù)隨水滴當量直徑的增加而增加。

旋翼； 結(jié)冰； 水滴撞擊特性； 歐拉法； 遮蔽區(qū)； 計算流體力學(xué)

直升機在結(jié)冰環(huán)境下飛行時，空氣中的過冷水滴會撞擊到槳葉表面，造成旋翼結(jié)冰現(xiàn)象。旋翼結(jié)冰是一個危害直升機飛行安全的重要因素，它會改變旋翼繞流流場，破壞旋翼氣動性能，進而影響直升機的操縱性，甚至導(dǎo)致空難的發(fā)生[1-3]。由于過冷水滴的運動軌跡決定了旋翼表面水滴局部收集系數(shù)，直接影響結(jié)冰情況，所以水滴撞擊特性計算是直升機旋翼結(jié)冰數(shù)值模擬中的一個重要環(huán)節(jié)。

目前，水滴撞擊特性的計算方法主要有拉格朗日法和歐拉法[4-8]。針對旋翼等復(fù)雜流動環(huán)境，拉格朗日法的實現(xiàn)過程較繁瑣，計算效率和精度也需要提高。同時，歐拉法可以應(yīng)用于非定常流場，固壁邊界可以移動或轉(zhuǎn)動，更適合于旋翼槳葉的水滴撞擊特性計算[9-10]。

在進行旋翼水滴撞擊特性計算時，研究者往往選擇忽略旋翼的三維流動效應(yīng)[11-14]。Narducci和Kreeger在旋翼結(jié)冰數(shù)值計算中，通過提取特征剖面的迎角和來流速度，對各剖面進行了二維的水滴撞擊特性計算[11]。這種簡化的方法給旋翼水滴撞擊特性計算和結(jié)冰預(yù)測提供了便利，但是由于忽略了旋翼三維流場效應(yīng)，計算得到的水滴撞擊特性將會與真實值產(chǎn)生一定的偏差，正如文獻[11]中的結(jié)果表明，槳葉中前段的結(jié)冰預(yù)測結(jié)果與試驗數(shù)據(jù)的誤差較大。因此，為獲得更準確的旋翼水滴撞擊特性，三維流場效應(yīng)需要被充分考慮。

在歐拉法的水滴撞擊特性求解過程中，由于控制方程的高度非線性，Bourgault等建立了基于體積分數(shù)的數(shù)值擴散項，以解決易出現(xiàn)的體積分數(shù)的數(shù)值振蕩問題[15]。Tong和Luke將其歸因于局部密度脈沖引起的高表觀密度(Apparent Density)區(qū)域，并提出了對應(yīng)的數(shù)值擴散系數(shù)[16]。楊勝華等進一步指出高表觀密度區(qū)域易出現(xiàn)于流經(jīng)機翼上下表面的兩股氣流交匯處，并提出自適應(yīng)的數(shù)值擴散項，以解決計算不穩(wěn)定、局部奇異解等數(shù)值問題[17]。在對固定翼的水滴撞擊特性計算中，研究者往往將背風面的壁面密度設(shè)定為極小值[8,18]，且由于相對來流速度較快，局部密度脈沖現(xiàn)象引起的高密度區(qū)域可以很快被抹平。但是，旋翼三維水滴撞擊流場具有槳尖渦干擾、相對槳葉的來流速度沿展向變化大、非定常流動等特征，氣流交匯區(qū)域?qū)⒏鼮轭l繁，局部密度脈沖現(xiàn)象所引起的高密度區(qū)域也將難以在迭代過程中被順利抹平，給旋翼三維水滴撞擊特性的計算增加了難度。

鑒于此，為獲得更準確的旋翼水滴撞擊特性，本文建立了適用于旋翼槳葉的三維水滴撞擊特性計算新方法，并提出一種遮蔽區(qū)擴散模型以解決局部密度脈沖對計算穩(wěn)定性的影響。為了驗證本文方法的有效性，采用了美國直升機結(jié)冰飛行試驗(The Helicopter Icing Flight Test)的旋翼結(jié)冰數(shù)據(jù)[19]。通過對比，遮蔽區(qū)擴散模型在能夠保證計算精度的同時，可有效避免局部密度脈沖對計算穩(wěn)定性的影響，并能夠節(jié)省計算資源的消耗，提高了水滴撞擊特性計算環(huán)節(jié)的效率。在此基礎(chǔ)上，針對不同的環(huán)境溫度和水滴當量直徑的參數(shù)條件，分別對三維旋翼水滴撞擊特性進行了計算和分析，獲得了一些有意義的結(jié)果。

1 旋翼三維水滴流場的求解

在水滴撞擊特性計算前，本文采用計算流體力學(xué)(CFD)方法對直升機旋翼三維流場進行了計算，其中，旋翼嵌套網(wǎng)格技術(shù)和旋翼CFD方法在文獻[20]基礎(chǔ)上進行了相應(yīng)的改進。

1.1 水滴流場控制方程

在建立歐拉法求解水滴撞擊流場前，先作如下假設(shè)[21]：① 水滴簡化為直徑為MVD(Median Volumetric Diameter)參數(shù)的球形，在運動過程中不會變形或被破壞；② 水滴之間沒有碰撞或合并現(xiàn)象，水滴撞擊在壁面上無飛濺現(xiàn)象；③ 水滴在空氣流場中沒有發(fā)生傳熱傳質(zhì)現(xiàn)象；④ 水滴上的作用力，只考慮重力和阻力，忽略其余非穩(wěn)態(tài)力。

引入水滴體積分數(shù)(Water Volume Fraction)α的概念，它表示單位體積內(nèi)水滴相所占有的體積，并由此構(gòu)造水滴的表觀密度ρd為

ρd=αρw

(1)

式中：ρw為水的密度。

通過水滴表觀密度，建立水滴撞擊流場的連續(xù)方程和動量方程。由于水滴與旋翼槳葉存在相對運動，將坐標系固連在旋轉(zhuǎn)槳葉上。水滴撞擊流場的積分形式控制方程為

(2)

源項由旋翼旋轉(zhuǎn)引起的旋轉(zhuǎn)通量和水滴上的作用力組成，它們的表達式為

(3)

式中：[uavawa]=qa為空氣流場的絕對速度；Ki為慣性因子;ω為旋翼轉(zhuǎn)速；g為重力加速度。

空氣流場對水滴的作用力大小受慣性因子Ki的影響，慣性因子的表達式為

(4)

式中：μa為空氣的動力黏度；dd為球形水滴的直徑；ρa為當?shù)乜諝饷芏?；CDRe/24為Schiller-Naumann模型中定義的阻力函數(shù)，根據(jù)文獻[17]的經(jīng)驗公式：

(5)

Re是相對雷諾數(shù)，其表達式為

(6)

在時間推進上，采用了五步顯式Runge-Kutta迭代格式。為解決源項在控制方程中占優(yōu)的問題，對時間步長進行了隱式處理。處理后的時間步長表達式為

(7)

式中：Δti為第i個單元的當?shù)貢r間步長；m為五步Runge-Kutta法的迭代次數(shù)。

1.2 邊界條件

水滴流場的邊界條件設(shè)置如下：

1) 來流邊界：qd=qa，α=LWC/ρw，LWC(Liquid Water Content)為液態(tài)水含量，表示單位體積內(nèi)所含的液態(tài)水的質(zhì)量。

2) 出流邊界：所有變量法向一階導(dǎo)數(shù)為0。

3) 壁面邊界條件：與空氣流場的固壁面不同，只允許水滴流出計算域，不允許水滴流入計算域。在不考慮水滴的碰撞飛濺現(xiàn)象后，水滴撞擊在壁面時，會凍結(jié)成冰或以水膜形式吸附在壁面上，相當于水滴流出了計算域。

上述這種壁面邊界條件可通過壁面法矢n與相對于槳葉的水滴速度矢量qd+qω的判定關(guān)系進行變換。若n·(qd+qω)≤0,則水滴撞擊壁面，流出計算域；若n·(qd+qω)>0,則水滴未撞擊壁面，以無滑移壁面邊界條件處理。

1.3 局部收集系數(shù)的計算

水滴局部收集系數(shù)β為旋翼表面局部區(qū)域?qū)嶋H所收集的水量與該區(qū)域可能收集水量最大值的比值。該參數(shù)可以反映水滴的撞擊范圍和撞擊區(qū)域內(nèi)的收集水量分布，是影響最終結(jié)冰量的重要參數(shù)，其計算式為

(8)

式中：α∞為遠場的水滴體積分數(shù)。

2 遮蔽區(qū)擴散模型

在旋翼三維水滴撞擊流場求解過程中，某些局部單元會出現(xiàn)水滴表觀密度的數(shù)值振蕩，這可能會導(dǎo)致求解的不穩(wěn)定。文獻[16]將這種問題歸源于歐拉法中對水滴速度值的單一假設(shè)。圖1給出了沿槳葉表面運動的水滴路徑圖。

在現(xiàn)實情況中，水滴在空間中分布稀疏，較難發(fā)生相互碰撞的現(xiàn)象。在尾流處交匯的時候，它們往往會互不干擾，彼此保持原有的狀態(tài)繼續(xù)運動。這正如拉格朗日法中，追蹤的不同水滴經(jīng)常會經(jīng)過同一單元格，彼此間并無影響。在歐拉法中，來自不同路徑的水滴交匯于同一單元格時，由于該單元只能以一個速度矢量來表示，這就使原本流向不同方向的水滴聚在一起運動，從而造成了一個突然增加的密度脈沖。這種局部密度脈沖現(xiàn)象在翼型或固定翼的流場中出現(xiàn)較少，它們引發(fā)的高密度區(qū)域往往會在迭代過程中被及時抹平。而在三維旋翼水滴撞擊流場中，尾流、槳尖渦、槳轂附近等氣流交匯頻繁的區(qū)域極易發(fā)生密度脈沖現(xiàn)象，從而影響計算的穩(wěn)定性，產(chǎn)生數(shù)值振蕩。

為避免上述情況的發(fā)生，本文提出了遮蔽區(qū)擴散模型。在迭代過程中，通過使用擴散模型，計算域內(nèi)會逐漸在背風面生成遮蔽區(qū)，并向下游擴散。同時，這些區(qū)域會被劃出計算域，減少了計算資源的消耗，并有效避免了密度脈沖帶來的數(shù)值問題。

圖1 沿槳葉表面運動的水滴路徑圖Fig.1 Paths of droplet moving along blade surface

在遮蔽區(qū)擴散模型中，引入遮蔽區(qū)權(quán)重變量shadow，變量shadow有如下定義：

1) shadow的初始化。在水滴撞擊流場求解前，需要對其進行初值處理，賦值為0。

2) shadow在槳葉網(wǎng)格壁面層的更新。在每一次計算迭代完后，都需要對壁面層的shadow進行重新賦值，并與前文壁面邊界條件的判斷同時進行。當n·(qd+qω)>0，表示沒有水滴撞擊到物面區(qū)域，即遮蔽區(qū)的底層，賦值為2；當n·(qd+qω)≤0，表示有水滴撞擊到物面區(qū)域，即槳葉表面的撞擊區(qū)，賦值為0。由此，隨著水滴撞擊流場速度的更新，槳葉表面上的撞擊區(qū)域和遮蔽區(qū)也在逐漸變化，并最終趨于穩(wěn)定。

3) shadow的傳遞。當?shù)赑個單元格的shadowP為0時，且鄰近單元存在shadowNi等于2時，賦值為1。這表示該單元緊鄰遮蔽區(qū)，有成為遮蔽區(qū)的可能性。

4) shadow的增值與減值處理。當單元格shadowP為1時，則根據(jù)單元格的水滴表觀密度的增/減分別對該單元格進行減/增值處理，其計算式為

(9)

式中：η為收斂因子，控制shadow值的變化速度，影響遮蔽區(qū)生成的速度大小，本文取值為0.1；αP為單元P的水滴體積分數(shù)；上標n和n-1分別表示t和t-Δt時刻的值。由此，在單元格內(nèi)，水滴表觀密度的減少量越大，單元格成為遮蔽區(qū)的可能性就越大，在擴散模型中體現(xiàn)為shadow值的增加。

5) shadow值的范圍限定。在計算迭代過程中，為減少不必要的計算量，shadow值范圍為[0，2]，2表示單元格成為遮蔽區(qū)的可能性最大，0表示可能性最小。

6) 遮蔽區(qū)的擴散判定。當單元格shadow值為2時，且滿足式(10)的條件，則該單元標記為遮蔽區(qū)。

(10)

式中：ε為比例因子，即定義遮蔽區(qū)內(nèi)單元與來流值的表觀密度之比，可以控制遮蔽區(qū)生成的范圍大小，本文取值為0.01；下標inf表示表觀密度的初始值。shadowcri為遮蔽區(qū)標準值，其范圍為(3,4]，本文取值為3.5，該數(shù)值的意義為：鄰近單元格中存在一個可能成為遮蔽區(qū)的單元時(shadowNi>1.5)，當前單元格才能確定為遮蔽區(qū)。

圖2給出了旋翼槳尖附近的遮蔽區(qū)在迭代計算中的擴散過程。為顯示清晰，遮蔽區(qū)的單元格均隔行隔列顯示。圖中可見：遮蔽區(qū)首先依附于槳葉表面，生成遮蔽區(qū)的大致范圍；隨著迭代步數(shù)的增加，遮蔽區(qū)底層向槳葉前緣推進，厚度也有所增加；在遮蔽區(qū)底層停止推進后，遮蔽區(qū)開始迅速增厚，并向后緣擴散。

圖2 遮蔽區(qū)在三維旋翼槳葉尖部的擴散過程Fig.2 Dispersion process of shadow zone over 3D rotor blade tip

3 計算結(jié)果與分析

3.1 NACA0012二維翼型的水滴撞擊特性計算

本文對NACA0012二維翼型進行了水滴撞擊特性的計算。翼型弦長c=0.533 4 m，來流馬赫數(shù)Ma=0.302，迎角AOA(Angle of Attack)分別選取了0°、4°、8°，水滴當量直徑MVD=20 μm，液態(tài)水含量為1 g/m3。

圖3給出了不同迎角下的水滴表觀密度ρd的分布云圖。圖中可見，在翼型背風面會形成較大范圍的遮蔽區(qū)域，該區(qū)域內(nèi)水滴表觀密度很低。隨著迎角的增加，上翼面的遮蔽區(qū)范圍逐漸增加，下翼面的遮蔽區(qū)范圍逐漸減小。

圖3 不同迎角下的水滴表觀密度等值圖Fig.3 Contours of droplet apparent density at different AOAs

圖4給出了不同迎角下翼型前緣的水滴局部收集系數(shù)的分布，橫坐標中的S為該單元中心點沿翼型表面到駐點處的距離。由于在二維翼型的水滴撞擊特性計算和結(jié)冰預(yù)測研究中，LEWICE軟件發(fā)展成熟，本文將計算結(jié)果與文獻中LEWICE軟件的計算數(shù)據(jù)[8]進行了對比。圖中可見，水滴局部收集系數(shù)的計算結(jié)果與參考結(jié)果吻合良好，驗證了本文方法在二維翼型水滴撞擊特性計算中的有效性。通過不同迎角的水滴局部收集系數(shù)的結(jié)果對比，可以再次發(fā)現(xiàn)，隨著迎角的增加，水滴撞擊范圍逐漸增加。同時，最大水滴局部收集系數(shù)不受迎角變化的影響，均保持在0.7附近，當迎角增加后，最大水滴收集處將偏離駐點，向上翼面移動。

圖4 不同迎角下水滴局部收集系數(shù)圖Fig.4 Local collection coefficients of droplets at different AOAs

3.2 UH-1H三維旋翼槳葉結(jié)冰環(huán)境下水滴撞擊特性計算

本文對UH-1H型貝爾直升機旋翼在懸停狀態(tài)下的水滴撞擊特性及結(jié)冰情況進行了預(yù)測。其中，槳葉參數(shù)、氣象條件、試驗數(shù)據(jù)和文獻預(yù)測數(shù)據(jù)均來自于文獻[19]。槳葉參數(shù)如下：翼型為NACA0012，翼型弦長為0.533 4 m，槳葉長度為7.315 2 m，槳葉轉(zhuǎn)速為33.9 rad/s，槳葉扭度為10.9°。氣象條件如下：壓強為101 300 Pa，液態(tài)水含量為0.7 g/m3，水滴平均直徑為30 μm，來流溫度為-19 ℃。

在旋翼流場和旋翼水滴撞擊流場的計算過程中，本文采用了同一套旋翼嵌套網(wǎng)格系統(tǒng)，如圖5所示。背景網(wǎng)格采用半圓柱型網(wǎng)格，K=1和K=108為對稱面。槳葉網(wǎng)格采用C-O型網(wǎng)格，并在槳葉前緣進行了加密處理。其中：I方向是指網(wǎng)格單元從槳葉剖面的尾跡開始沿翼型表面以順時針方向遍歷一周；J方向是指網(wǎng)格單元從翼型最底層向最外層遍歷；K方向是指網(wǎng)格單元從槳根向槳尖遍歷。

圖5 懸停狀態(tài)下的旋翼嵌套網(wǎng)格系統(tǒng)Fig.5 Embedded grid system for rotor in hovering flight

圖6 旋翼槳葉三維水滴撞擊特性圖Fig.6 Sketch of 3-D droplet impingement on rotor blade

圖6給出了通過歐拉法計算得到的水滴撞擊特性的示意圖。在槳葉前緣形成了水滴撞擊區(qū)域，其中，撞擊區(qū)域主要存在于槳葉的下表面；遮蔽區(qū)范圍從撞擊區(qū)域邊緣開始，貼附于槳葉表面，向后緣延伸，并逐漸增厚，如圖中水滴體積分數(shù)α的分布所示。

圖7給出了旋翼槳葉表面的水滴局部收集系數(shù)分布圖和3個特征剖面的水滴表觀密度的分布云圖。從整體上看，水滴局部收集系數(shù)在槳葉前半段沿展向呈增加趨勢。從槳根到r=0.35R剖面處，水滴撞擊面積和局部收集系數(shù)迅速增加，變化顯著；在槳葉中段附近，水滴局部收集系數(shù)增加趨勢減緩；在r=0.7R附近，水滴局部收集系數(shù)開始減少，又在r=0.85R(鄰近槳尖)附近時再次增加。

圖7 槳葉表面的水滴局部收集系數(shù)分布圖 Fig.7 Distribution of local collection coefficient of droplets on blade surface

圖8 旋翼三維水滴局部收集系數(shù)計算結(jié)果Fig.8 Calculated results of local collection coefficients of 3-D droplets on rotor blade

為了進一步顯示旋翼三維水滴撞擊特性，圖8 給出了水滴局部收集系數(shù)的數(shù)值曲線圖。圖8(a)給出了3個典型剖面的水滴局部收集系數(shù)的分布，可以看出最大水滴局部收集系數(shù)均分布在駐點附近，同時下翼面撞擊區(qū)域較大，上翼面撞擊區(qū)域較?。浑S著徑向距離的增加，水滴撞擊范圍在弦向上有所增加。圖8(b)給出了最大水滴局部收集系數(shù)沿展向的分布情況，從整體上看，最大水滴局部收集系數(shù)在展向上有近50%的變化幅度。其中，最大水滴局部收集系數(shù)在槳根附近增加速度較快，從0.6增加到了0.85；在0.35R到0.7R間，最大水滴局部收集系數(shù)增加量較少，僅增加了0.1；隨后，最大水滴局部收集系數(shù)先于0.83R處減少到0.6，又在槳尖處回升至0.85附近。

最大水滴局部收集系數(shù)在r=0.8R附近的下降現(xiàn)象主要源于旋翼槳尖渦的影響，圖9(a)給出了旋翼流場的槳尖渦和水滴撞擊流場的低表觀密度區(qū)域(3.5×10-7kg/m3)。槳尖渦經(jīng)過槳葉的尾流區(qū)，這導(dǎo)致了槳尖渦區(qū)域的水滴表觀密度較低。同時，在槳尖渦的誘導(dǎo)作用下，大部分的水滴不經(jīng)過槳尖渦區(qū)域，如圖中剖面所示。由于旋翼做旋轉(zhuǎn)運動，這個低水滴表觀密度的槳尖渦區(qū)域會對后一片槳葉產(chǎn)生很大影響。圖9(b)所示，在r=0.8R剖面的前緣附近，水滴表觀密度較低，這導(dǎo)致了最大水滴局部收集系數(shù)在0.83R處的波谷。

圖9 旋翼槳尖渦對水滴撞擊流場表觀密度分布的影響Fig.9 Influence of blade tip vortex on distribution of droplet apparent density in rotor droplet flowfield

圖10給出了y=0平面的水滴表觀密度的局部分布云圖，圖中ABCD區(qū)域內(nèi)的黑色箭頭表示水滴運動軌跡在y=0平面上的投影。在槳根附近(圖中A所示)，水滴運動方向變化顯著，具有明顯的三維特征，如果采用二維翼型剖面的簡化計算方法，可能會與真實值產(chǎn)生較大誤差。在r=0.45R附近(z=3.2 m，圖中B所示)，水滴運動軌跡仍具有明顯的方向變化，后文選擇了這里的剖面冰形進行了對比驗證。在靠近槳尖附近(圖中C、D所示)，由于相對來流速度的增加、展向速度的影響減小，水滴運動軌跡的方向主要沿x軸方向，此時，與二維翼型的水滴運動軌跡特點一致。

由于該試驗中并未直接記錄槳葉表面的水滴局部收集系數(shù)的試驗值，這里通過計算冰形與試驗值的對比來間接驗證三維水滴撞擊特性計算新方法的可靠性。在結(jié)冰計算環(huán)節(jié)中，采用了文獻[22]的結(jié)冰預(yù)測方法。圖11給出了展向最大結(jié)冰厚度的對比(兩種測量方法的試驗值來自文獻[19]，結(jié)冰時間為180 s)。圖12給出了槳葉剖面的結(jié)冰預(yù)測結(jié)果(結(jié)冰時間為180 s，計算冰形來自文獻[11]的計算結(jié)果，試驗冰形來自文獻[19]中的Tracings方法)。由圖11可見，在徑向最大結(jié)冰厚度的對比上，計算值沿徑向的變化趨勢與試驗值一致，且在各剖面上，預(yù)測冰形與試驗值均吻合良好。由于旋翼槳葉的氣動加熱，槳葉的表面溫度沿展向逐漸上升。剖面r=0.45R處的溫度較低，所結(jié)冰型接近于霜冰類型，結(jié)冰量可以在一定程度上反映水滴收集系數(shù)的大小。因此，該剖面冰形可以用來驗證水滴局部收集系數(shù)的計算結(jié)果。由圖12(a)可見，在結(jié)冰范圍和結(jié)冰量上，計算結(jié)果均與文獻[19]的試驗值一致，并優(yōu)于文獻[11]的計算結(jié)果。需要指出的是，文獻[11]的計算結(jié)果是基于二維剖面的水滴撞擊特性計算的方法，而r=0.45R剖面處于圖10的區(qū)域B，具有較明顯的三維效應(yīng)。因此，通過與文獻[11]計算結(jié)果在r=0.45R的剖面冰形對比，進一步說明了旋翼三維效應(yīng)在水滴撞擊特性計算中的重要性。

圖10 y=0剖面上的水滴表觀密度分布圖Fig.10 Distribution of droplet apparent density in y=0 plane

圖11 旋翼最大結(jié)冰厚度分布圖Fig.11 Maximum depth of ice accretion on rotor

圖12 不同徑向位置的剖面冰形圖 Fig.12 Ice shape on different sections for different radial position

3.3 遮蔽區(qū)擴散模型在水滴撞擊特性計算過程中的影響

為了討論遮蔽區(qū)擴散模型對水滴撞擊特性計算的影響，分別給出了二維翼型和三維旋翼槳葉情況下的典型算例。

圖13給出了前文AOA為0°時二維翼型算例的計算情況。在加入遮蔽區(qū)模型后，殘值在580步附近即達到10-3量級，比未加入模型情況提前了約150步，節(jié)省了近22%的計算時間。同時，通過水滴局部收集系數(shù)的對比，可以看出，遮蔽區(qū)模型的引入，并不會影響二維水滴撞擊特性的數(shù)值計算結(jié)果。

圖14(a)給出了遮蔽區(qū)單元的動態(tài)生成過程。圖中可見，在計算初期，遮蔽區(qū)單元生成速度較快，隨著迭代步數(shù)的推進，遮蔽區(qū)單元增加速度逐漸減少，在約500步附近，達到最大值。之后，遮蔽區(qū)單元個數(shù)將保持穩(wěn)定，遮蔽區(qū)也不再擴散。圖中也分別給出了第200步、第300步和第500步的翼型尾緣的局部遮蔽區(qū)域分布情況。圖14(b)給出了遮蔽區(qū)域整體的動態(tài)擴散過程。在第200步時，遮蔽區(qū)主要依附于翼型表面，此時，遮蔽區(qū)在整個計算域中所占比例較?。辉诘?00步時，翼型表面的遮蔽區(qū)域逐漸增厚，并有向后方擴散的趨勢；在第500步時，遮蔽區(qū)單元進一步向尾流方向擴散，同時，翼型表面遮蔽區(qū)的厚度增加速度減緩。

圖13 遮蔽區(qū)擴散模型對二維翼型計算結(jié)果的影響Fig.13 Influence of shadow zone model on calculated results of 2-D airfoil

圖14 遮蔽區(qū)單元擴散過程Fig.14 Dispersion process of cells among shadow zone

圖15是三維旋翼槳葉的水滴撞擊流場殘值對數(shù)收斂曲線。在未加入遮蔽區(qū)模型時，殘值在第55迭代步開始發(fā)散，不能繼續(xù)進行計算，這正是因為尾流附近發(fā)生了密度脈沖。在加入遮蔽區(qū)擴散模型后，由于三維槳葉遮蔽區(qū)域的底層生成數(shù)量較多，殘值在經(jīng)過初期發(fā)生的小規(guī)模振蕩后，迅速收斂。

圖15 遮蔽區(qū)擴散模型對旋翼三維算例結(jié)果的影響Fig.15 Influence of shadow zone model on calculated results of 3-D rotor

圖16 未加入遮蔽區(qū)擴散模型的第55迭代步的計算結(jié)果Fig.16 Calculated results at the 55th iterative step without present model

考慮到未加入遮蔽區(qū)擴散模型的算例在第55迭代步后發(fā)散，圖16給出了當前迭代步的水滴流場數(shù)據(jù)。由圖可見，由于槳葉上下來流的交匯，在槳葉尾流后呈現(xiàn)大面積的高密度區(qū)域。在r=0.15R處剖面的水滴表觀密度示意圖中，尾流出現(xiàn)呈帶狀的高密度區(qū)域，其中，某些單元的表觀密度高達來流值的40倍。

加入遮蔽區(qū)擴散模型后，密度脈沖引起的高密度區(qū)域面積減小，如圖17所示。第55步時，與未加入模型的結(jié)果相比，高密度區(qū)域雖然存在，但面積和數(shù)值均有所降低。同時，因遮蔽區(qū)域的生成，從槳葉尾緣上脫出的低密度區(qū)域面積有所增加，可以有效抑制密度脈沖現(xiàn)象的發(fā)生，滿足數(shù)值計算的穩(wěn)定性要求。當?shù)綌?shù)為100步時，由于遮蔽區(qū)域的進一步擴散，高密度區(qū)域消失。此時，在尾流區(qū)將不再出現(xiàn)因密度脈沖現(xiàn)象引起的高密度區(qū)域，增強了計算的穩(wěn)定性。

圖17 加入擴散模型后槳葉剖面表觀密度及遮蔽區(qū)的計算結(jié)果(r=0.15R)Fig.17 Calculated results of apparent density and shadow zone of blade with present model (r=0.15R)

3.4 不同氣象參數(shù)下的旋翼槳葉水滴撞擊特性計算

為了分析氣象參數(shù)對旋翼水滴撞擊特性的影響，本文在不同的水滴平均當量直徑和環(huán)境溫度的情況下(其余氣象參數(shù)與前文一致)，對懸停狀態(tài)下的三維旋翼進行了計算。

表1給出了水滴當量直徑(MVD)對水滴撞擊特性的影響。從表1中看出，MVD為15 μm時，水滴撞擊面積占整個槳葉表面積的2.6%，當MVD增加到30 μm時，撞擊面積比例增加到6.75%，單位時間的水滴收集量也隨之增加。可見，MVD參數(shù)對旋翼水滴撞擊特性的影響很大。

表1 不同MVD下的水滴撞擊特性Table 1 Droplet impingement property at different MVDs

圖18給出了在不同MVD參數(shù)下的三維旋翼表面水滴局部收集系數(shù)的分布云圖。圖中可見，隨著MVD參數(shù)的增加，水滴的慣性(保持原有運動方向的能力)隨之增強，這使得水滴更容易撞擊到旋翼槳葉表面，導(dǎo)致了水滴撞擊區(qū)域擴大、水滴局部收集系數(shù)增加。

表2給出了環(huán)境溫度對水滴撞擊特性的影響。在易發(fā)生結(jié)冰的環(huán)境溫度下，隨溫度的下降，水滴收集量和撞擊面積都有小幅度的增加?？傮w上，水滴撞擊特性受溫度的影響較小。

圖18 不同MVD參數(shù)下的水滴局部收集系數(shù)分布圖 Fig.18 Distribution of local collection coefficient at different MVDs

表2不同溫度下的水滴撞擊特性

Table2Dropletimpingementpropertyatdifferenttemperatures

CaseT/℃Amountofimpingingwater/(10-3g·s-1)Proportionofimpactareatosurfacearea/%1-47．1826．692-87．1906．703-147．2466．754-197．2956．835-247．3296．89

4 結(jié) 論

1) 建立了一種基于歐拉法的三維旋翼水滴撞擊特性計算新方法。通過對NACA0012翼型和UH-1H旋翼的驗證計算，表明了該方法的有效性。同時，提出的遮蔽區(qū)擴散模型有效解決了密度脈沖所引起的穩(wěn)定性和收斂性問題，并大幅減少了水滴撞擊特性環(huán)節(jié)的計算時間。

2) 在二維翼型上，最大水滴局部收集系數(shù)受來流迎角的變化影響較小，主要受來流馬赫數(shù)影響。在三維旋翼槳葉上，隨著徑向距離的增加(相對來流馬赫數(shù)變大)，旋翼槳葉剖面撞擊范圍有所增加，最大水滴局部收集系數(shù)會呈現(xiàn)先增加后減少再增加的變化趨勢。同時，在槳尖附近，水滴運動軌跡三維特征較??；在槳葉中段至槳根附近，水滴運動軌跡具有較明顯的三維特征，用二維剖面擬合的方法在該部位可能會偏離真實值。

3) 隨著MVD參數(shù)的增加，旋翼槳葉表面的水滴撞擊區(qū)域擴大，水滴局部收集系數(shù)有所增加，并可能導(dǎo)致后續(xù)結(jié)冰量的增加。

[1] FLEMMING R J. The past twenty years of icing research and development at Sikorsky aircraft: AIAA-2002-0238[R]. Reston: AIAA, 2002.

[2] FLEMMING R J, ALLDRIDGE P, DOEPPNER R. Artificial icing tests of the S-92 helicopter in the McKinley Climatic Laboratory: AIAA-2004-0737[R]. Reston: AIAA, 2004.

[3] BRITTON R K, BOND T H, FLEMMING R J. An overview of a model rotor icing test in the NASA Lewis icing research tunnel: AIAA-1994-0716[R]. Reston: AIAA, 1994.

[4] WRIGHT W B, POTAPCZUK M G, LEVINSON L H. Comparison of LEWICE and GlennICE in the SLD regime: AIAA-2008-0439[R]. Reston: AIAA, 2008.

[5] IULIANO E, BRANDI V, MINGIONE G, et al. Water impingement prediction on multi-element airfoils by means of Eulerian and Lagrangian approach with viscous and inviscid air flow: AIAA-2006-1270[R]. Reston: AIAA, 2006.

[6] 胡劍平, 劉振俠, 張麗芬. 發(fā)動機整流支板大尺寸過冷水滴撞擊特性[J]. 航空學(xué)報, 2011, 32(10): 1778-1785.

HU J P, LIU Z X, ZHANG L F.Supercooled large droplet impact behaviors on an aero-engine strut[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(10): 1778-1785 (in Chinese).

[7] 易賢, 王開春, 桂業(yè)偉, 等. 結(jié)冰面水滴收集率歐拉計算方法研究及應(yīng)用[J]. 空氣動力學(xué)學(xué)報, 2010, 28(5): 596-601.

YI X, WANG K C, GUI Y W, et al. Study on Eulerian method for icing collection efficiency computation and its application[J]. Acta Aerodynamica Sinica, 2010, 28(5): 596-601 (in Chinese).

[8] KIM J W, DENNIS P G, SANKAR L N, et al. Ice accretion modeling using an Eulerian approach for droplet impingement: AIAA-2013-0246[R]. Reston: AIAA, 2013.

[9] KINZEL M P, SAROFEEN C M, NOACK R W, et al. A finite-volume approach to modeling ice accretion: AIAA-2010-4230[R]. Reston: AIAA, 2010.

[10] BAIN J, SANKAR L N,GARZA D. A methodology for the prediction of rotor blade ice formation and shedding[C]//Proceedings of the SAE 2011 Aircraft and Engine icing and Ground Testing Conference. Warrendale, PA: SAE International, 2011.

[11] NARDUCCI R, KREEGER R E.Analysis of a hovering rotor in icing conditions: NASA/TM-2012-217126[R]. Washington, D.C.: NASA, 2012.

[12] CAO Y H, LI G Z, ZHONG G. Tandem helicopter trim and flight characteristics in the icing condition[J]. Journal of Aircraft, 2010, 47(5): 1559-1569.

[13] RAJMOHAN N, BAIN J, NUCCI M, et al. Icing studies for the UH-60A rotor in forward flight[C]//AHS Aeromechanics Specialists’ Conference. Fairfax, VA: American Helicopter Society International, 2010.

[14] 張威, 林永峰, 陳平劍. 基于歐拉方法的旋翼翼型結(jié)冰數(shù)值模擬及參數(shù)影響[J]. 南京航空航天大學(xué)學(xué)報, 2011, 43(3): 375-380.

ZHANG W, LIN Y F, CHEN P J. Numerical simulation of ice accretion and parameter effects based on Eulerian droplet model[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2011, 43(3): 375-380 (in Chinese).

[15] BOURGAULT Y, HABASHI W G, DOMPIERRE J, et al. A Eulerian approach to supercooled droplets impingement calculations: AIAA-1997-0176[R]. Reston: AIAA, 1997.

[16] TONG X L, LUKE E A. Eulerian simulations of icing collection efficiency using a singularity diffusion model: AIAA-2005-1246[R]. Reston: AIAA, 2005.

[17] 楊勝華, 林貴平, 申曉斌. 三維復(fù)雜表面水滴撞擊特性計算[J]. 航空動力學(xué)報, 2010, 25(2): 284-290.

YANG S H, LIN G P, SHEN X B. Water droplet impingement prediction for three-dimensional complex surfaces[J]. Journal of Aerospace Power, 2010, 25(2): 284-290 (in Chinese).

[18] JUNG K Y, JUNG S K, MYONG R S, et al. A three-dimensional unstructured finite volume method for droplet impingement in aircraft icing: AIAA-2013-2576[R]. Reston: AIAA, 2013.

[19] LEE J D, HARDING R, PALKO R L. Documentation of ice shapes on the main rotor of a UH-1H helicopter in hover: NASA-CR-168332[R]. Washington, D.C.: NASA, 1984.

[20] 招啟軍, 徐國華. 新型槳尖旋翼懸停氣動性能試驗及數(shù)值研究[J]. 航空學(xué)報, 2009, 30(3): 422-429.

ZHAO Q J, XU G H. Aerodynamic performance of rotor with new type blade-tip in hover based upon test and numerical investigations[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(3): 422-429 (in Chinese).

[21] CAO Y H, MA C, ZHANG Q, et al. Numerical simulation of ice accretions on an aircraft wing[J]. Aerospace Science and Technology, 2012, 23(1): 296-304.

[22] 陳希, 招啟軍, 趙國慶. 計入離心力影響的直升機旋翼翼型結(jié)冰數(shù)值模擬[J]. 航空動力學(xué)報, 2014, 29(9): 2158-2165.

CHEN X, ZHAO Q J, ZHAO G Q. Numerical simulation for ice accretion on rotor airfoil of helicopter considering the effects of centrifugal force[J]. Journal of Aerospace Power, 2014, 29(9): 2158-2165 (in Chinese).

(責任編輯： 鮑亞平，徐曉)

AeronauticsandAstronautics,Nanjing210016,China

New method for predicting3-D water droplet impingement on rotorconsidering influence of shadow zone

CHENXi,ZHAOQijun*

NationalKeyLaboratoryofScienceandTechnologyonRotorcraftAeromechanics,NanjingUniversityof

For the complex 3-D viscous flowfield of the helicopter rotor, a new 3-D Eulerian method is established for calculating the droplet impingement on the rotor. The Computional Fluid Dynamics (CFD) simulation method for predicting the rotor flowfield is developed based on the embedded grid method. Considering the 3-D effect of the rotor fully, the 3-D flowfield of droplets on the same embedded grids is solved by the Eulerian method to overcome the defects of the traditional 2-D calculation method. To overcome the problem of numerical stability and convergence by the density impulse in the wake area, a shadow zone dispersion model is proposed. In the model, a shadow variable is used to control the generation and dispersion of the shadow cells. The new Eulerian method is validated by comparing the calculation and experiment results of NACA0012 airfoil and UH-1H rotor in hover. In addition, the effects of the atmospheric temperature and Median Volumetric Diameter (MVD ) of the droplet on the droplet impingement on the rotor are calculated and analyzed. Results show that the computation time of the 2-D droplet impingement property is reduced by 22% and the calculation stability of 3-D droplet impingement property is improved by using the shadow zone dispersion model. The droplet impingement area increases with the spanwise distance on the rotor blade. The maximum droplet local collection efficiency presents the change tendency of ‘increase-decrease-increase’, and the variation amplitude is close to 50%. The droplet impingement area and the collection efficiency increase with the increase of the MVD.

rotor; ice accretion; water droplet impingement property; Eulerian method; shadow zone; CFD

2016-09-01；Revised2016-10-11；Accepted2016-11-25；Publishedonline2016-12-051640

NationalNaturalScienceFoundationofChina(11272150)

2016-09-01；退修日期2016-10-11；錄用日期2016-11-25； < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2016-12-051640

www.cnki.net/kcms/detail/11.1929.V.20161205.1640.004.html

國家自然科學(xué)基金 (11272150)

*

.E-mailzhaoqijun@nuaa.edu.cn

陳希， 招啟軍． 考慮遮蔽區(qū)影響的旋翼三維水滴撞擊特性計算新方法J． 航空學(xué)報,2017,38(6):120745.CHENX,ZHAOQJ.Newmethodforpredicting3-DwaterdropletimpingementonrotorconsideringinfluenceofshadowzoneJ.ActaAeronauticaetAstronauticaSinica,2017,38(6):120745.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0313

V211.3

A

1000-6893(2017)06-120745-12

URL：www.cnki.net/kcms/detail/11.1929.V.20161205.1640.004.html

*Correspondingauthor.E-mailzhaoqijun@nuaa.edu.cn