劉建明，蔣向華，王 東，馬永峰

（1.中航工業(yè)沈陽發(fā)動機設計研究所，沈陽 110015；2.北京航空航天大學能源與動力工程學院北京 100191）

0 引言

鳥類與飛機等飛行器之間的碰撞稱為鳥撞。文獻[1]顯示，發(fā)動機是飛機較易受鳥撞擊的部位。隨著航空發(fā)動機技術的發(fā)展，在未來的新發(fā)動機研制中，由于葉尖切向速度增加，葉片厚度日趨變薄，葉片采用新結(jié)構(gòu)（整體葉盤、整體葉環(huán)等）和先進材料（金屬基、樹脂基復合材料等），使發(fā)動機在推力增加、質(zhì)量減輕的同時，葉片的抗鳥撞能力也可能受到影響，從而影響發(fā)動機的可靠性和安全性[2]?；谥袊娇瞻l(fā)動機研制的實際情況，前期采用大量的數(shù)值模擬，最后進行少量的試驗驗證是比較合理的發(fā)展模式[3]。目前，用于鳥撞數(shù)值模擬的大型通用有限元程序有Pam-Crash[1,4]、Abaqus[5]和 MSC.Dytran[6-7]等，還有專門針對鳥撞發(fā)展起來的程序，如 NOSAPM[8]、PW/WHAM[9]、MAGNA[10]等。目前，在采用流固耦合算法研究葉片鳥撞的問題上，普遍采用殼元劃分葉片，而殼元在模擬非薄壁結(jié)構(gòu)問題上存在困難。如果采用實體元劃分葉片，在失效單元的處理及耦合面的推進上，流固耦合程序還缺乏完善的理論。

本文針對目前通用流固耦合算法在模擬實體元結(jié)構(gòu)破壞上存在的不足，基于鳥撞鋁板的試驗結(jié)果，采用MSC.Patran軟件建立了鳥撞實體元平板的有限元模型，同時研究和驗證了實體元平板鳥撞數(shù)值模擬方法，進一步模擬了考慮破壞的實體元空心葉片的鳥撞過程。

1 流固耦合方法研究及試驗驗證

為了進行實體元空心葉片的鳥撞數(shù)值模擬，本文采取以實體元自由表面的節(jié)點作為殼元節(jié)點的方法，建立了覆蓋在實體元表面的1層4節(jié)點殼元。殼元起輔助作用，用來形成封閉的耦合面，并向?qū)嶓w元傳遞力，而自身對原結(jié)構(gòu)的影響較小。

為了驗證方法的可行性，針對文獻[11]所做的鳥撞鋁板試驗，建立了鳥撞鋁板模型。其中鋁板尺寸為410mm×500mm×10mm，4邊固支。鋁板用8節(jié)點6面體實體元劃分，沿厚度方向等厚劃分2層單元；在厚度方向的4個面上建立1圈啞元；以鋁板實體元的上下自由表面的節(jié)點為公共節(jié)點，分別建立1層4節(jié)點四邊形殼元；啞元和殼元形成包圍實體元的、封閉的耦合面。所建立的鳥撞鋁板模型如圖1所示，圖中歐拉域2由包圍鋁板的歐拉網(wǎng)格組成，歐拉域1由歐拉域2之外不發(fā)生鳥撞的歐拉網(wǎng)格組成。

圖1 鳥撞鋁板模型

鋁板用分段線性塑性材料模擬，材料參數(shù)見表1，表中：ρ為密度；E為彈性模量；σy為靜態(tài)屈服強度；μ為泊松比。取殼元與實體元的材料參數(shù)一致進行初始計算。鳥體采用線性流體模型，質(zhì)量M=1.8 kg，密度ρ=928.15 kg/m3，體積模量 K=2200 MPa，用長徑比為2∶1的二端半球中間圓柱實體來模擬，以30°角斜撞向鋁板中央。試驗測得了鋁板中心點的法向位移隨時間的變化情況，以及撞擊開始1.24 ms后沿鋁板寬度方向的測點的位移。

表1 鋁板材料參數(shù)

1.1 殼元厚度h的確定

考慮殼元厚度h對精度的影響，對此只改變h進行計算。分別計算了 h=0.01、0.1、1、0mm（實體元表面未鋪殼元）4種情況下鋁板中心點的法向位移隨時間的變化情況，并與試驗值進行對比，如圖2所示。從圖中可見，當h=0.01mm，即實體元與h之比為102數(shù)量級時，計算值與試驗值的誤差較小。因此，選取h=0.01mm。

圖2 不同殼元厚度時板中心點的位移-時間曲線

1.2 殼元彈性模量E s的影響

為了驗證殼元材料的Es對計算精度的影響，只改變殼元的Es，分別計算Es=108000、36000MPa時鋁板中心點的法向位移隨時間的變化情況，并與試驗值進行對比，如圖3所示。

計算表明，3種殼元彈性模量所對應的鋁板中心點的法向位移的差別很小。

根據(jù)以上計算，當殼元厚度取為0.01 mm，材料參數(shù)與實體元一致時，撞擊開始1.24ms后沿鋁板寬度方向的測點的位移的計算值與試驗值的對比，如圖4所示。從圖中可見，二者趨勢比較一致。

1.3 殼元最大塑性應變εs的影響

鑒于文獻[11]所做的鳥撞試驗，鋁板沒有發(fā)生破壞，為了考察實體元（鋪殼）結(jié)構(gòu)發(fā)生破壞時、殼元最大塑性應變εs對計算結(jié)果的影響，建立了鳥體撞擊鈦合金實體元（鋪殼）平板有限元模型。平板尺寸為100 mm×120 mm×1.2 mm，4周固支，材料參數(shù)見表2，εp為實體元最大塑性應變。鳥體半徑r=15 mm，長度L=50mm，以600m/s的速度正撞擊平板中央。

表2 鈦合金材料參數(shù)

當 εs=0.05、0.01、0.1時，對實體元（鋪殼）平板的失效情況進行計算，并與εs=0.05、殼元劃分的平板的基準計算結(jié)果進行對比。某時刻2種劃分時平板的失效如圖5所示。由圖可以清楚地看到，當εs=εp時，實體元(鋪殼)和殼元的失效情況比較一致。因此，計算時取 εs=εp較為合理。

計算結(jié)果表明，對于實體元劃分且厚度不大的葉片，可以采用在實體元的自由表面鋪1層殼元來進行葉片鳥撞數(shù)值模擬，且實體元與殼元的厚度之比為102數(shù)量級，2種單元的材料參數(shù)應一致。殼元的作用是傳遞鳥體與實體元之間的相互作用力。如果實體元表面沒有鋪殼元，目前的Dytran還不具備模擬實體元失效的功能。若實體元在鳥撞過程中失效，則該單元被刪除。

圖5 某時刻實體元(鋪殼)和殼元平板的失效

2 葉片鳥撞有限元建模

采用Patran建立了8節(jié)點6面體空心葉片鳥撞模型。為節(jié)省計算時間，只建立了發(fā)動機中與鳥體有初始碰撞行為的第1級轉(zhuǎn)子葉片模型；鳥體采用歐拉單元模擬。為了進行考慮葉片破壞的鳥撞數(shù)值模擬，參考第1節(jié)中所用的方法，以葉片的葉盆和葉背的實體元自由表面上的節(jié)點為公共節(jié)點，建立了覆蓋在實體元表面的1層4節(jié)點殼元，h=0.01mm（實體元與h比為102數(shù)量級）；葉片的葉尖、葉根以及2葉片之間建立了1層啞元，如圖6所示。鳥撞有限元模型如圖7所示。在圖7中，歐拉域2由自適應網(wǎng)格組成，由包圍相鄰葉片區(qū)域的歐拉網(wǎng)格以及包圍葉片的歐拉網(wǎng)格組成，如圖8所示。歐拉域1包含了葉片轉(zhuǎn)動范圍之外的區(qū)域。歐拉域1和歐拉域2通過啞元加以分隔。

設鳥體質(zhì)量 M=0.35 kg，密度 ρ=930 kg/m3，體積模量K=2200 MPa。文獻[12]的研究表明：圓柱體鳥體和兩端半球中間圓柱鳥體對計算結(jié)果的影響較小，參照航空發(fā)動機適航規(guī)定，用長徑比為2∶1的圓柱體代替鳥體，其速度V=102.8m/s。葉片受離心載荷作用，轉(zhuǎn)速為12000 r/min，初始狀態(tài)為穩(wěn)定旋轉(zhuǎn)，由Nastran的預應力分析計算獲得。

圖8 自適應網(wǎng)格

葉片材料模型見參考文獻[3]，實體元和殼元材料一致，材料參數(shù)見表3。D、P為表征動態(tài)硬化的材料系數(shù)，D=100、P=10，EH為硬化模量。

表3 葉片材料參數(shù)

3 不考慮葉片失效的鳥撞數(shù)值模擬

鳥體的材料特性是鳥撞仿真分析的重點和難點，真實的鳥體有骨有血有肉，鳥體的本構(gòu)方程很難描述[12]。本文分析了實體元表面殼元厚度，鳥體的密度、體積模量，以及葉片材料的彈性模量、屈服應力、硬化模量的變化對葉片鳥撞響應的影響。葉片和鳥體材料計算以第2節(jié)中的數(shù)據(jù)為基準，只改變單個參數(shù)，其他不變。

3.1 殼元厚度及單元材料的驗證

圖9 鳥撞瞬間

為了驗證第1節(jié)中所用方法的準確性，計算了殼元厚度h=0.01、0.1 mm的2種情況下葉片應力峰值最大的單元應力（簡稱葉片應力）隨時間的變化，并與未鋪殼實體元葉片（h=0）的結(jié)果對比。某時刻的鳥撞瞬間如圖9所示，不同h的葉片應力曲線如圖10所示。從圖9中可見，由于鳥撞產(chǎn)生的應力值超過了葉片的屈服應力，葉片產(chǎn)生局部塑性變形，巨大的相互作用力使得鳥體發(fā)生流變。從圖10中可見，當h=0.01mm時，結(jié)果相對h=0的誤差較小，為1.3%，從而驗證了實體元與h之比為102數(shù)量級、材料一致的準確性。以下取h=0.01mm進行計算。文獻[13]將鳥撞分為4個階段：即初始撞擊、壓力衰減、恒定流動、流動終止。葉片應力響應趨勢與鳥體正撞擊剛性靶的應力趨勢比較一致（圖10）。

圖10 不同殼元厚度的葉片應力-時間曲線

3.2 鳥體密度ρ的影響

選取鳥體密度 ρ=1300、930、500 kg/m3進行計算。3種ρ鳥體撞擊下葉片應力曲線如圖11所示。從圖中可見，ρ對葉片響應的影響較大。葉片應力峰值隨ρ的增大而增大，因為ρ越大，與葉片相撞的鳥體的質(zhì)量越大，使得撞擊能量越大，應力峰值也越大。

圖11 不同鳥體密度的葉片應力-時間曲線

3.3 鳥體體積模量的影響

選取鳥體體積模量K=2200、6000、10000 MPa進行計算。3種K鳥體撞擊下葉片應力曲線如圖12所示。從圖中可見，鳥體K對葉片響應的影響較小。

圖12 不同鳥體體積模量的葉片應力-時間曲線

3.4 葉片彈性模量的影響

選取彈性模量 E=200000、112500、56250 MPa進行計算。3種E的葉片應力曲線如圖13所示。從圖中可見，E對葉片初始撞擊的應力峰值影響較小，但對恒定流動的應力峰值的影響較大，E越大，恒定流動的應力峰值也越大。因此，在其他條件相同的情況下，E大的葉片能吸收更多的鳥體能量。

圖13 不同葉片彈性模量的葉片應力-時間曲線

3.5 葉片屈服應力σs的影響

選取屈服應力σs=400、900、1300 MPa進行計算。3種σs的葉片應力曲線如圖14所示。從圖中可見，σs對葉片應力響應的影響較大，σs越大，葉片初始撞擊和恒定流動的應力峰值也越大。

圖14 不同葉片屈服應力的葉片應力-時間曲線

3.6 葉片硬化模量E H的影響

選取硬化模量 EH=14286、7000、3000 MPa進行計算。3種EH的葉片應力曲線如圖15所示。從圖中可見，葉片EH對葉片應力響應的影響較大，EH越大，初始撞擊的應力峰值也越大，而恒定流動的應力峰值越小。

4 考慮葉片失效的鳥撞數(shù)值模擬

由于葉片的轉(zhuǎn)速很高，撞擊過程時間較短，因此，鳥體撞擊葉片產(chǎn)生的應力和應變很大，從而易造成葉片的損傷和破壞。為了更加準確真實地模擬鳥體對葉片的撞擊過程，給葉片添加最大塑性應變破壞準則進行計算，最大塑性應變εb=0.05。葉片破壞時的鳥撞瞬間如圖16所示。從圖中可見，葉片前緣由于受到撞擊而造成破壞。

5 結(jié)論

（1）基于鳥撞鋁板的試驗結(jié)果，驗證了在實體元表面鋪一層殼元進行流固耦合計算的可行性。計算結(jié)果表明，當實體元與殼元的厚度之比為102數(shù)量級，且2種單元的材料參數(shù)一致時，鋪1層殼元對計算結(jié)果的影響較小。

（2）鳥體密度對葉片應力響應的影響較大，葉片應力峰值隨鳥體密度的增大而增大；鳥體體積模量對葉片應力響應的影響較小。

（3）葉片彈性模量對葉片初始撞擊的應力峰值影響較小，但對恒定流動的應力峰值的影響較大，彈性模量越大，恒定流動的應力峰值也越大。葉片屈服應力、硬化模量對葉片應力響應的影響較大，屈服應力越大，葉片初始撞擊和恒定流動的應力峰值也越大；硬化模量越大，初始撞擊的應力峰值越大，而恒定流動的應力峰值越小。

（4）考慮葉片失效的計算進一步模擬了鳥撞對葉片的真實損傷過程。

[1]劉軍,李玉龍,劉元鏞.基于SPH方法的葉片鳥撞數(shù)值模擬研究[J].振動與沖擊,2008,27（9）：90-93.

LIU Jun,LI Yulong,LIU Yuanyong.Numerical simulation study of bird-impact on a blade using SPH method[J].Journal of Vibration and Shock,2008,27（9）：90-93.(in Chinese)

[2]陳偉,高德平,尹晶.航空發(fā)動機葉片的鳥撞擊損傷研究[J].燃氣渦輪試驗與研究,1998,11（4）：34-39.

CHEN Wei,GAO Deping,YIN Jin.Damage investigation of bird impact on areoengine blades[J].Gas Turbine Experiment and Research,1998,11（4）：34-39.(in Chinese)

[3]蔣向華,王延榮.采用流固耦合方法的整級葉片鳥撞擊數(shù)值模擬[J].航空動力學報,2008,23（2）：299-304.

JIANG Xianghua,WANG Yanrong.Numerical simulation of bird impact on bladed rotor stage by fluid-solid coupling method[J].Journalof Aerospace Power,2008,23（2）：299-304.(in Chinese)

[4]劉軍,李玉龍,徐緋.基于PAM-CRASH的鳥撞飛機風擋動響應分析[J].爆炸與沖擊,2009,29（1）：80-84.

LIU Jun,LI Yulong,XU Fei.Dynamic response analysis of bird-impact aircraft windshields based on PAM-CRASH[J].Explosion and Shock Waves,2009,29（1）：80-84.(in Chinese)

[5]ZHU Shuhua,TONGMingbo,Wang Yuequan.Experiment and numerical simulation of a full-scale aircraft windshield subjected to bird impact[R].AIAA-2009-2575.

[6]劉建明,蔣向華.材料參數(shù)對葉片鳥撞動響應影響數(shù)值模擬[J].航空發(fā)動機,2010,36（5）：36-38.

LIU Jianming,JIANG Xianghua.Numerical simulation of blade material effect on dynamic response of bird impact on flat blade[J].Aeroengine,2010,36（5）：36-38.(in Chinese)

[7]劉建明,蔣向華,秦銀雷等.采用流固耦合的實體元空心葉片鳥撞數(shù)值模擬 [J].航空動力學報,2010,25(10),2211-2216.

LIU Jianming,JIANG Xianghua,QIN Yinlei,et al.Numerical simulation of bird impact on solid-element hollow blades by using fluid-solid coupling method[J].Journal of Aerospace Power,2010,25(10),2211-2216.

[8]Storace A F,Nimmer R P,Ravenhall R.Analytical and experiment investigation of bird impact on fan and compressor blading[J].Journalof Aircraft,1984,21（7）:520-527.

[9]Teichman H C,Tadros R N.Analytical and experimental simulation of fan blade behaviour and damage under bird impact[J].Journal of Engineering for Gas Turbines and Power,1991,113（3）：582-594.

[10]McCarty R E,Hart J L.Validation of the MAGNA computer program for nonlinear finite element analysis of aircraft transparency bird impact[R].AD-A-140701,1983.

[11]白金澤.基于神經(jīng)網(wǎng)絡方法的鳥撞飛機風擋反問題研究[D].西安:西北工業(yè)大學,2003.

BAI Jinze.Inverse issue study of bird-impact to aircraft windshield based on neural network method[D].Xi'an:Northwestern Polytechnical University,2003.(in Chinese)

[12]朱書華,童明波.鳥體形狀對飛機風擋鳥撞動響應的影響[J].南京航空航天大學學報,2008,40（4）:551-555.

ZHU Shuhua,TONG Mingbo.Bird shape sensitivity to dynamic response of bird strike on aircraft windshield[J].Journal of Nanjing University of Aeronautics and Astronautics,2008,40（4）：551-555.(in Chinese)

[13]Barber JP.Taylor H R.Wilbeck JS.Bird impact forces and pressures on rigid and compliant targets[R].AD-A-061313,1978.