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

陳希， 招啟軍

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

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

陳希， 招啟軍*

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

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

旋翼； 結冰； 水滴撞擊特性； 歐拉法； 遮蔽區； 計算流體力學

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

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

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

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

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

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

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

1.1 水滴流場控制方程

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

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

ρd=αρw

(1)

式中：ρw為水的密度。

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

(2)

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

(3)

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

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

(4)

式中：μa為空氣的動力黏度；dd為球形水滴的直徑；ρa為當地空氣密度；CDRe/24為Schiller-Naumann模型中定義的阻力函數，根據文獻[17]的經驗公式：

(5)

Re是相對雷諾數，其表達式為

(6)

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

(7)

式中：Δti為第i個單元的當地時間步長；m為五步Runge-Kutta法的迭代次數。

1.2 邊界條件

水滴流場的邊界條件設置如下：

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

2) 出流邊界：所有變量法向一階導數為0。

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

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

1.3 局部收集系數的計算

水滴局部收集系數β為旋翼表面局部區域實際所收集的水量與該區域可能收集水量最大值的比值。該參數可以反映水滴的撞擊范圍和撞擊區域內的收集水量分布，是影響最終結冰量的重要參數，其計算式為

(8)

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

2 遮蔽區擴散模型

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

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

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

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

在遮蔽區擴散模型中，引入遮蔽區權重變量shadow，變量shadow有如下定義：

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

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

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

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

(9)

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

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

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

(10)

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

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

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

3 計算結果與分析

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

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

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

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

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

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

3.2 UH-1H三維旋翼槳葉結冰環境下水滴撞擊特性計算

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

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

圖5 懸停狀態下的旋翼嵌套網格系統Fig.5 Embedded grid system for rotor in hovering flight

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表2給出了環境溫度對水滴撞擊特性的影響。在易發生結冰的環境溫度下，隨溫度的下降，水滴收集量和撞擊面積都有小幅度的增加。總體上，水滴撞擊特性受溫度的影響較小。

圖18 不同MVD參數下的水滴局部收集系數分布圖 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 結 論

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

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

3) 隨著MVD參數的增加，旋翼槳葉表面的水滴撞擊區域擴大，水滴局部收集系數有所增加，并可能導致后續結冰量的增加。

[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] 胡劍平, 劉振俠, 張麗芬. 發動機整流支板大尺寸過冷水滴撞擊特性[J]. 航空學報, 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] 易賢, 王開春, 桂業偉, 等. 結冰面水滴收集率歐拉計算方法研究及應用[J]. 空氣動力學學報, 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] 張威, 林永峰, 陳平劍. 基于歐拉方法的旋翼翼型結冰數值模擬及參數影響[J]. 南京航空航天大學學報, 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] 楊勝華, 林貴平, 申曉斌. 三維復雜表面水滴撞擊特性計算[J]. 航空動力學報, 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] 招啟軍, 徐國華. 新型槳尖旋翼懸停氣動性能試驗及數值研究[J]. 航空學報, 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] 陳希, 招啟軍, 趙國慶. 計入離心力影響的直升機旋翼翼型結冰數值模擬[J]. 航空動力學報, 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">網絡出版時間

時間：2016-12-051640

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

國家自然科學基金 (11272150)

*

.E-mailzhaoqijun@nuaa.edu.cn

陳希， 招啟軍． 考慮遮蔽區影響的旋翼三維水滴撞擊特性計算新方法J． 航空學報,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