林長亮，王益鋒，王浩文，陳仁良，尚曉冬

(1.南京航空航天大學 直升機旋翼動力學國家級重點實驗室，南京 210016;2.中國商用飛機有限責任公司 上海飛機設計研究院，上海 200232;3.清華大學 航天航空學院，北京 100084;4.哈爾濱飛機工業(yè)集團有限責任公司 飛機設計研究所，哈爾濱 150066)

鳥撞輕則使機體損傷，重則會造成災難性后果，直接威脅人員的生命安全。Dolbeer等［1］在第8屆鳥撞會議上指出“直升機鳥撞事故主要發(fā)生在距地面高度600 m以下，風擋和旋翼是鳥撞事故發(fā)生的主要部位”。低空飛行是直升機的顯著使用特點，因此直升機發(fā)生鳥撞的可能性很大，并且隨著飛行速度的提高，鳥撞事故的危害性也在逐步加大。旋翼是直升機的升力面、推力裝置和操縱面，是區(qū)別于固定翼飛機的主要特征。在直升機飛行過程中旋翼一直處于高速旋轉的狀態(tài)，一旦發(fā)生鳥類撞擊，就會對飛行安全會造成嚴重的威脅，是抗鳥撞研究中需要著重解決的關鍵技術。

對于鳥撞問題的研究，在飛機鳥撞方面，國內(nèi)外學者作了大量的研究工作。Barber等［2-3］采用10%孔隙率的明膠代替真鳥進行試驗研究，發(fā)現(xiàn)鳥撞過程可被描述成一個非恒定的流體動力學過程。Meguid等［4］通過計算研究了鳥體形狀對計算結果的影響，發(fā)現(xiàn)鳥體和目標的最初接觸面積對接觸力的峰值有巨大影響。劉軍等［5-6］進行了鳥撞平板試驗研究，使用神經(jīng)網(wǎng)絡方法對試驗中的鳥體參數(shù)進行反演計算。陳偉和關玉璞等［7-9］在發(fā)動機葉片鳥撞的載荷模型、瞬態(tài)響應計算、試驗方法等方面開展了大量工作。Zhu等［10］對鳥撞飛機全尺寸風擋進行了試驗研究和數(shù)值模擬。萬小朋、李京菁等［11-12］對飛機機翼前緣的抗鳥撞性能進行了研究。

在直升機鳥撞研究方面，僅有美國西科斯基公司在S-92直升機抗鳥撞設計中，對垂尾前緣、尾槳鳥撞進行了動力學仿真，并開展了鳥撞試驗驗證研究［13-14］。目前，我國對直升機旋翼鳥撞問題進行的研究很少。王益鋒等［15］采用Hertz接觸理論處理局部彈性變形，運用Hamilton原理建立了直升機旋翼槳葉的彈性碰撞的動力學模型。溫海濤、林長亮等［16-17］使用有限元軟件對直升機主槳葉鳥撞的過程進行了數(shù)值模擬。

直升機旋翼在旋轉過程中，槳葉載荷工況非常復雜，除槳葉自身重力外，不僅受到氣動升力作用，還承受離心載荷。同時，由于采用槳轂結構形式的不同，槳葉的揮舞、擺振、扭轉運動還存在不同的非線性耦合。采用有限元軟件很難模擬出實際的旋翼工作特點。利用旋翼綜合氣彈分析程序LORA［19-20］，考慮旋翼槳轂結構特點和旋翼旋轉過程中槳葉的揮舞、擺振、扭轉運動的非線性耦合，求解出直升機旋翼槳葉在飛行過程中的穩(wěn)態(tài)響應，以此作為鳥體撞擊的初始狀態(tài)，然后采用非線性流-固耦合算法，考慮載荷與響應之間耦合對槳葉響應的影響，利用Newmark數(shù)值積分的方法求解槳葉的動態(tài)響應，并分析了鳥體速度、質量、撞擊位置、槳葉根部鉸約束、離心力等系統(tǒng)參數(shù)對槳葉響應的影響，從而為直升機旋翼的抗鳥撞設計提供一定的理論依據(jù)。

1 槳葉鳥撞動力學建模

在槳葉鳥撞過程中，槳葉是可變形柔性靶體，其變形會影響撞擊載荷的大小與分布。為準確預估槳葉鳥撞擊響應，因此采用流-固耦合算法［8］計算槳葉鳥撞的動態(tài)響應。

1.1 流－固耦合鳥撞載荷模型

在鳥撞擊載荷作用期間，在ti時刻，假設鳥體以初始速度，入射角θi(θi=90°時為正撞擊)撞擊槳葉，ti+1=ti+Δt時刻，槳葉發(fā)生變形，入射角變?yōu)?θi+1，剖面入射角示意圖如圖1所示。

圖1 剖面撞擊示意圖Fig.1 Profile sketch of bird impact on the blade

入射角θi+1可以表示為:

以Δt時間內(nèi)消耗的鳥體質量Δmi+1為研究對象，應用動量定理可得:

式中:ρ和Ai+1分別為鳥的密度和橫截面積。撞擊載荷的作用時間由撞擊時鳥體在撞擊過程中消耗的長度Si決定:

當Si等于鳥體的長度L時，撞擊過程結束。

圖2為面積隨鳥體消耗長度的示意圖，鳥體橫截面積呈三角形形式變化，在S/L=0.2處面積為最大，a、b分別為橢圓形撞擊區(qū)域的短軸與長軸，當正撞擊時，a=b=1.5R，R為鳥體的初始半徑。

圖2 面積A隨長度的變化Fig.2 Area“A”variation with length

1.2 載荷空間分布

載荷類型采用分布載荷。當斜撞擊時，撞擊區(qū)域為一個橢圓，方程為(Ysinθ)2/R2+Z2/R2=1。式中Y、Z分別沿橢圓的長、短軸方向，碰撞載荷在撞擊區(qū)域上的平均分布力為:

鳥撞擊載荷初始沖擊壓力雖然很高，但持續(xù)時間極短，而恒定流動壓力雖然不高，但持續(xù)時間較長，鳥撞載荷的沖量主要是在此階段傳遞給槳葉的［18］。因此，一般在建立鳥撞擊載荷模型時根據(jù)恒定流動壓力的分布情況來選擇載荷的空間分布形式。本文中載荷空間分布形式為壓力系數(shù)沿長軸線性分布，沿短軸均勻分布，載荷空間分布如圖3所示。具體分布形式可以表示為:

沿橢圓短軸a，

沿橢圓長軸b，

圖3 載荷空間分布示意圖Fig.3 Load spatial distribution map

1.3 槳葉鳥撞動力學方程

根據(jù)Hamilton原理，引入Chopra有限元法對槳葉單元進行離散，組集動能、應變能、空氣動力及撞擊載荷產(chǎn)生的虛功，得到基于廣義力形式的槳葉非線性動力學隱式微分方程:

式中:T、E、A、P表示動能、應變能、氣動力以及撞擊載荷引起的廣義力，在計算過程中，由于撞擊過程時間持續(xù)較短，假設撞擊過程中入流和槳葉迎角沒有發(fā)生變化。

1.4 槳葉動態(tài)響應求解步驟

本文采用有限轉動梁理論處理槳葉的變形［19-20］，在翼型氣動力計算時采用Leishman-Beddoes二維非定常和動態(tài)失速模型計算出槳葉剖面升力、阻力和俯仰力矩，旋翼誘導速度計算采用Glauert入流改進模型，采用Chopra的15自由度梁元對槳葉進行離散，根據(jù)Hamilton原理得出非線性隱式旋翼動力學方程，利用隱式Newmark數(shù)值算法對動力學方程進行了求解得出槳葉的穩(wěn)態(tài)解，在每個時間步NTIME上，先預估得到葉片的響應(位移、速度和加速度)值，接著利用ASSMK模塊求解槳葉的響應，當槳葉單元為NEIMP撞擊單元時，進入ADDIT模塊，計算撞擊載荷值和作用時間，并將載荷和力矩施加到撞擊單元上，這樣通過ADDIT模塊完成在槳葉動力學模型上施加鳥體分布載荷的過程，然后利用NMARK1模塊計算出槳葉響應值，并且在每個時間步上進行迭代，直到響應收斂，得出這個時間步上的槳葉響應值，再將計算結果返回程序進行下一個時間步的計算。如此反復，直到撞擊仿真時間結束為止。流程圖如圖4所示。

圖4 槳葉響應計算流程圖Fig.4 The calculation flowchart of Blade response

2 分析方法驗證

由于目前國內(nèi)外對于直升機旋翼鳥撞問題的研究比較少，還沒有公開的詳細試驗數(shù)據(jù)作為驗證依據(jù)。并且在理論計算方面，直升機旋翼鳥撞問題屬于軟體瞬態(tài)沖擊旋轉梁范疇，而現(xiàn)有的國內(nèi)外文獻主要集中于剛性質量與梁的碰撞研究，對軟體與梁的碰撞問題的研究非常少。因此，為了能夠驗證本文計算方法的可靠性，主要從以下兩個方面進行了驗證。

2.1 鳥體載荷模型驗證

為驗證本文采用鳥體載荷模型的正確性，針對文獻［21］中平板葉片鳥撞試驗進行了驗證。文獻［21］給出了質量為0.161 5 kg的鳥體垂直撞擊平板葉片根部的高速攝影結果及葉尖Z向位移測量結果。在計算中，鳥體直徑為50 mm，長度為80 mm，密度為1 028 kg/m3，速度為156.1 m/s垂直撞擊平板葉片，撞擊中心距葉片根部為60 mm。葉片材料參數(shù)見表1所示。葉片采用8節(jié)點殼單元進行劃分，葉片有限元網(wǎng)格劃分如圖5所示。

表1 葉片材料參數(shù)Tab.1 Blade material parameters

圖5 葉片有限元模型Fig.5 The Finite element model of blade

圖6 葉尖位移時間歷程曲線Fig.6 Tip displacement time history curve

圖6給出了本文與試驗的葉尖位移對比結果?？梢钥闯?，計算結果中含有高階諧波成分，沒有試驗值曲線光滑，主要是由于本文沒有對計算結果進行濾波處理。試驗中，葉尖的最大位移為70.5 mm，本文計算結果中的最大位移為74.7 mm，誤差為5.9%，總體上講，計算結果與試驗是比較吻合的，從而驗證了本文所采用鳥撞載荷模型的有效性。

2.2 旋翼動力學模型驗證

為驗證旋翼動力學模型在計算瞬態(tài)響應時的正確性，本文選取 Keller實驗作為數(shù)值算例［22］。Keller所采用的試件為與美國H-46運輸直升機槳葉1/8Froude數(shù)相似的模型槳葉，槳葉長 1.006 m，重0.998 kg。實驗裝置如圖7所示。整個實驗過程為:槳葉初始處于靜止狀態(tài)，根部固支，揮舞限動角為0°，初始揮舞角為2°至9.7°，在自重作用下槳葉有初始位移。突然釋放根部約束，槳葉下墜，當揮舞角到達0°時，根部與揮舞限動塊發(fā)生碰撞。

圖8～9分別給出了初始揮舞角為4°時，槳尖位移和距槳葉根部距離20%處槳葉上應變的變化曲線。與實驗相比，本文槳尖最大位移誤差為5.1%，20%處最大應變誤差為6.56%。通過與試驗值的比較驗證了本文建模方法的正確性。

圖7 Keller槳葉揚起下墜實驗裝置Fig.7 Configuration of blade droop stop impact test

圖8 槳尖位移隨時間的變化圖Fig.8 Blade tip displacement changing with time

3 槳葉鳥撞動態(tài)響應計算

旋翼翼型為NACA0012，直徑為3 m，弦長0.3 m，槳葉揮舞剛度3.78×105N·m2，擺振剛度9.45×106N·m2，扭轉剛度 3.78 ×106N·m2，拉伸剛度 1.26 ×109N·m2，槳葉線密度48.6 kg/m。在鳥撞擊瞬態(tài)響應研究中，把模擬鳥看做長徑比為2的圓柱體［7，10］，鳥體直徑0.05 m，質量 0.17 kg，鳥體以 120 m/s的相對速度與槳葉發(fā)生撞擊，鳥體入射角為60°，撞擊位置取槳葉彈性軸上徑向位置x/L=0.4處，分析時間為0.04 s。

首先采用流-固耦合和流-固非耦合兩種算法分別計算了槳葉的動態(tài)響應。圖10給出了槳葉鳥撞載荷曲線。當考慮載荷與響應之間的耦合效應時，碰撞載荷峰值為1.91×105N，載荷作用時間為 1.06 ms，當不考慮耦合響應時，碰撞載荷峰值為2.0×105N，載荷作用時間為1.0 ms;圖11給出了槳葉碰撞位置位移隨時間的變化曲線，可以看出，考慮耦合效應的局部彈性變形要稍小于不考慮耦合效應的局部變形，隨著時間的推移，兩種算法幾乎同時到達向下位移峰值。通過以上比較可以看出，采用流固耦合算法時，雖然載荷作用時間延長，但是由于槳葉變形起到了卸載作用，因此使得撞擊載荷、槳葉的整體變形要小于不考慮耦合效應時的結果;同時可以看出，算法的選擇對槳葉到達最大變形的時間以及槳葉的振動形式影響不大。

圖9 x/L=0.2處應變隨時間的變化圖Fig.9 The strain changes with time of the position x/L=0.2

圖10 碰撞載荷曲線Fig.10 The impact load curve

圖11 撞擊位置位移曲線Fig.11 The displacement curve of impact position

4 系統(tǒng)參數(shù)對槳葉響應的影響

4.1 速度影響分析

鳥體的質量和速度，決定了鳥撞載荷沖量的大小，為研究速度對槳葉響應的影響，保持其他參數(shù)不變，單獨改變鳥體的初始速度，對槳葉的響應進行了分析.圖12給出了不同撞擊速度的載荷曲線，圖13為不同速度下的碰撞位置位移曲線?？梢钥闯鲭S著鳥體速度的增加，碰撞載荷變大，載荷作用時間減小，槳葉撞擊位置局部彈性變形和整體彈性變形增大。

4.2 質量影響分析

分別選取質量0.15 kg，0.17 kg，0.19 kg，采用流 -固耦合算法，對槳葉鳥撞響應進行了分析，計算結果如圖14、15所示。

可以看出鳥體質量對槳葉的鳥撞響應影響較明顯，響應的差異基本與質量成比例關系，單獨增加鳥體質量時，碰撞載荷增加，撞擊位置位移增大。

4.3 撞擊位置影響分析

為比較槳葉對不同加載位置的響應，分別計算了三組算例。撞擊位置分別取在x/L=0.1，x/L=0.4處，x/L=0.7處。計算結果如圖16、17所示。

圖12 不同速度下的碰撞載荷曲線Fig.12 The impact load curve with different speed

圖13 不同速度下的碰撞位置位移曲線Fig.13 The displacement curve of impact location under different speed

圖14 不同質量的碰撞載荷曲線Fig.14 The impact load curve with different quality

圖15 不同質量的碰撞位置位移曲線Fig.15 The displacement curve of impact location under different quality

圖16 不同撞擊位置的碰撞位置位移曲線Fig.16 The displacement curve of impact location with different impact position

圖17 不同撞擊位置的槳尖位移曲線Fig.17 The displacement curve of blade tip with different impact position

圖18 不同約束的碰撞位置位移曲線Fig.18 The displacement curve of impact location with different constraint condition

圖19 不同約束的槳根揮舞彎矩曲線Fig.19 The flap moment curve of blade root with different constraint condition

圖20 不同轉速的碰撞載荷曲線Fig.20 The impact load curve with different rotation speed

從圖16、17碰撞位置位移、槳尖位移時間歷程中可以看出，在撞擊載荷作用期間，槳葉局部彈性變形相差不大，但是隨著時間的推移，當撞擊點離槳根越遠時，碰撞位置、槳尖位移越大。對于x/L=0.1處的撞擊，碰撞位置處的位移振動后迅速返回平衡位置，但槳尖位移存在明顯的滯后，這是由于從撞擊力在撞擊處產(chǎn)生加速度到槳尖呈現(xiàn)位移，中間經(jīng)過了槳葉的彈性變形及加速度兩次對時間的積分過程，因此槳葉尖位移出現(xiàn)了滯后。

4.4 根部約束影響分析

為了研究根部鉸約束對槳葉響應的影響，保持其他參數(shù)不變，單獨改變根部約束形式，采用鉸接式同時根部施加角彈簧，角彈簧剛度為6.0×102N/rad。從圖18碰撞位置位移時間歷程中可以看出，由于根部鉸約束剛度的降低，槳葉剛性運動影響增強，整體彈性位移增大。從圖19槳根揮舞彎矩隨時間的變化曲線可以看出，雖然槳葉整體位移增大，但是由于槳葉剛性運動的影響，槳根揮舞彎矩減小。

4.5 旋翼轉速影響分析

當槳葉旋轉時，離心力對槳葉的響應有很大的影響。保持其他參數(shù)不變，單獨改變旋轉角速度，對槳葉鳥撞響應進行了計算。計算結果如圖20、21所示。

從圖20碰撞載荷變化的時間歷程可以看出，當旋轉角速度增加時，碰撞載荷增大，這是由于旋轉角速度的增加導致鳥體相對槳葉的撞擊速度增大，但是迎角相對減小。

圖21 不同轉速的碰撞位置位移曲線Fig.21 The displacement curve of impact location with different rotation speed

從圖21碰撞位置位移時間歷程可以看出，當旋轉角速度增加時，由于離心力的作用，槳葉的局部剛度和整體剛度得到增強，槳葉的局部彈性變形變小，碰撞位置位移減小，當載荷作用結束后，碰撞位置位移明顯減小。

5 結論

(1)以旋翼綜合氣彈分析程序LORA為基礎，采用流固算法，建立了直升機槳葉鳥撞動力學方程，并利用數(shù)值積分求解出了鳥撞槳葉的動態(tài)響應，為研究直升機旋翼鳥撞問題提供了一種參考方法。

(2)采用流固耦合算法時，作用時間延長，但碰撞載荷以及槳葉局部彈性變形相對要小一些，同時發(fā)現(xiàn)算法的選擇對槳葉到達最大變形的時間以及振動形式影響不大。

(3)研究了鳥體速度、質量、撞擊位置、槳葉根部約束和離心力等參數(shù)對槳葉動態(tài)響應的影響，從而為直升機槳葉抗鳥撞設計提供一些理論依據(jù)。

［1］Dolbeer R A，Wright S E，Cleary E C.Bird strikes to civil helicopters in the United States［C］.1990-2005.the 8th meeting ofBirdStrikeCommitteeUSA/Canada，2006:45-50.

［2］ Barber J P，Taylor H R，Wilbeck J S.Characterization of bird impact on a rigid plate:partⅠ［R］.AFFDL-TR-75-5，AD A021142，1975:1-44.

［3］ Barber J P，Wilbeck J S，Taylor H R.Bird impact forces and pressures on rigid and compliant targets［R］.AFFDL-TR-77-60，AD A 061313，1978:1-78.

［4］Meguid S，Mao R，Ng T.FE analysis of geometry effects of an artificial bird striking an aeroengine fan blade［J］.International Journal of Impact Engineering，2008，35(6):487-498.

［5］劉 軍，李玉龍.鳥體本構模型參數(shù)反演Ⅰ:鳥撞平板試驗研究［J］.航空學報，2011，32(5):802-811.

LIU Jun， LIYu-long. Parameters inversion on bird constitutive model PartⅠ:study on experiment of bird striking on plate［J］.Acta Aeronautica et Astronautica Sinica，2011，32(5):802-811.

［6］劉 軍，李玉龍，石霄鵬，等.鳥體本構模型參數(shù)反演Ⅱ:模型參數(shù)反演研究［J］.航空學報，2011，32(5):812-821.

LIU Jun，LI Yu-long，SHI Xiao-peng，et al.Parameters inversion on bird constitutive model PartⅡ:study on model parameters inversion［J］.Acta Aeronautica et Astronautica Sinica，2011，32(5):812-821.

［7］陳 偉，關玉璞，高德平.發(fā)動機葉片鳥撞擊瞬態(tài)響應的數(shù)值模擬［J］，航空學報，2003，24(6):531-533.

CHEN Wei， GUANYu-pu， GAODe-ping. Numerical simulation of the transient response of blade due to bird impact［J］.Acta Aeronautica et Astronautica Sinica，2003，24(6):531-533.

［8］陳 偉，高德平.載荷與響應耦合下葉片鳥撞擊響應分析［J］，航空動力學報，1998，13(1):93-96.

CHEN Wei，GAO De-ping.Transient responses of blades to bird impact coupled with deformations of blade［J］.Journal of Aerospace Power，1998，13(1):93-96.

［9］關玉璞，陳 偉.粒子分離器渦流葉片鳥撞擊損傷試驗［J］，航空動力學報，2007，22(12):2094-2095.

GUAN Yu-pu，CHEN Wei. Experimental study on bird damage of vortex vane in the particle separator［J］.Journal of Aerospace Power，2007，22(12):2094-2095.

［10］ Zhu S H，Tong M B，Wang Y Q.Experiment and numerical simulation of a full scale aircraft windshield subjected to bird impact［R］.AIAA，2009:2009-2575.

［11］萬小朋，趙美英.基LS-DYNA的飛機機翼前緣抗鳥撞分析［J］.西北工業(yè)大學學報，2007，25(2):285-289.

WAN Xiao-peng，ZHAO Mei-ying.A method for calculating anti-bird impact capability of an aircraft wing using LS-DYNA Software［J］. Journal of Northwestern Polytechnical University，2007，25(2):285-289.

［12］李京菁，趙美英.機翼前緣抗鳥撞設計仿真方法驗證［J］.機械科學與技術，2011，30(10):1757-1760.

LI Jing-jing，ZHAO Mei-ying.Verification of the simulation methods to bird impact resistance design of wing leading edge［J］.Mechanical Science and Technology for Aerospace Engineering，2011，30(10):1757-1760.

［13］ Dobyns A，F(xiàn)ederici F，Young R.Bird strike analysis and test of a spinning S-92 tail rotor［C］.American Helicopter Society 57th Annual forum.May 9-11，2001.

［14］ Dobyns A.Bird strike analysis of S-92 vertical tail cover using DYTRAN［R］.American Helicopter Society.October 7，1998.

［15］王益鋒，王浩文，高 正.直升機旋翼槳葉的彈性碰撞動力學建模［J］.航空動力學報，2009，24(5):2046-2050.

WANG Yi-feng，WANG Hao-wen，GAO Zheng.Dynamic model of elastic impact of helicopter rotor blades［J］.Journal of Aerospace Power，2009，24(5):2046-2050.

［16］溫海濤，關玉璞，高德平.直升機主槳葉的鳥撞有限元數(shù)值模擬［J］.航空動力學報，2009，24(5):1150-1158.

WEN Hai-tao，GUAN Yu-pu，GAO De-ping.Numerical simulations of bird impact on helicopter rotor blades［J］.Journal of Aerospace Power，2009，24(5):1150-1158.

［17］林長亮，王浩文，陳仁良.采用流固耦合方法的直升機槳葉鳥撞數(shù)值模擬［J］.科學技術與工程，2012，12(1):1-6.

LIN Chang-liang， WANG Hao-wen， CHEN Ren-liang.Numerical simulation of bird impact on the helicopter blade by fluid-structure coupling method［J］.Science Technology and Engineering，2012，12(1):1-6.

［18］尹 晶.鳥撞擊的載荷因素對葉片響應的影響［J］.航空動力學報，1992，7(1):51-54.

YIN Jing.Effect of some load factors of bird impact on blade reponse［J］.Journal of Aerospace Power，1992，7(1):51-54.

［19］王益鋒，王浩文，高 正.一種復合材料柔性梁的響應計算方法［J］.南京航空航天大學學報，2008，40(2):180-184.

WANG Yi-feng，WANG Hao-wen，GAO Zheng.Method for calculating response of composite flexible beams［J］.Journal of Nanjing University of Aeronautics and Astronautics，2008，40(2):180-184.

［20］王浩文，高 正，鄭兆昌.前飛狀態(tài)下直升機旋翼系統(tǒng)氣彈響應及穩(wěn)定性分析［J］.振動工程學報，1999，12(4):521-528.

WANG Hao-wen， GAO Zheng， ZHENG Zhao-chang.Aeroelastic response and stability of helicopter rotor blades in forward flight［J］.Journal of Vibration Engineering，1999，12(4):521-528.

［21］文穎娟.葉片鳥撞擊響應分析?；夹g與驗證研究［D］.南京:南京航空航天大學，2006.

［22］Keller J A，Smith E C.Experiment/theoretical correlation of analysis for helicopter rotor blade/droop stop impacts［R］.AIAA-97-1094，1997.