實體元空心葉片鳥撞流固耦合研究及數值模擬

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

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

0 引言

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

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

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

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

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

圖1 鳥撞鋁板模型

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

表1 鋁板材料參數

1.1 殼元厚度h的確定

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

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

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

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

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

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

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

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

表2 鈦合金材料參數

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

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

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

2 葉片鳥撞有限元建模

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

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

圖8 自適應網格

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

表3 葉片材料參數

3 不考慮葉片失效的鳥撞數值模擬

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

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

圖9 鳥撞瞬間

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

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

3.2 鳥體密度ρ的影響

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

圖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 考慮葉片失效的鳥撞數值模擬

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

5 結論

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

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

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

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

[1]劉軍,李玉龍,劉元鏞.基于SPH方法的葉片鳥撞數值模擬研究[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]陳偉,高德平,尹晶.航空發動機葉片的鳥撞擊損傷研究[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]蔣向華,王延榮.采用流固耦合方法的整級葉片鳥撞擊數值模擬[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]劉建明,蔣向華.材料參數對葉片鳥撞動響應影響數值模擬[J].航空發動機,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]劉建明,蔣向華,秦銀雷等.采用流固耦合的實體元空心葉片鳥撞數值模擬 [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]白金澤.基于神經網絡方法的鳥撞飛機風擋反問題研究[D].西安:西北工業大學,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.