謝燦軍， 童明波， 劉 富， 郭亞洲， 朱書華

(1.南京航空航天大學(xué) 飛行器先進(jìn)設(shè)計(jì)技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南京 210016;2.中國商飛 上海飛機(jī)設(shè)計(jì)研究院，上海 201210;3.西北工業(yè)大學(xué) 航空學(xué)院，西安 710012)

民用飛機(jī)結(jié)構(gòu)的鳥撞試驗(yàn)是適航取證的重要一環(huán)，合理準(zhǔn)確的鳥撞數(shù)值分析能預(yù)先指導(dǎo)飛機(jī)結(jié)構(gòu)抗鳥撞設(shè)計(jì)，可極大減少適航試驗(yàn)費(fèi)用［1］。SPH方法作為基于拉格朗日技術(shù)的自適應(yīng)無網(wǎng)格粒子法，與有限元方法耦合用于鳥撞分析，能有效求解強(qiáng)非線性、強(qiáng)耦合問題［2－4］。

民用飛機(jī)適航條款規(guī)定飛機(jī)尾翼結(jié)構(gòu)必須能承受3.63kg鳥重撞擊。即在輕量化設(shè)計(jì)目標(biāo)下，對尾翼結(jié)構(gòu)鳥撞數(shù)值分析可靠性提出更高要求。故考慮瞬態(tài)高速?zèng)_擊下鋁合金材料的應(yīng)變率效應(yīng)，建立準(zhǔn)確的材料本構(gòu)模型;尾翼結(jié)構(gòu)中存在大量連接件，如何合理模擬結(jié)構(gòu)間連接亦會(huì)影響計(jì)算精度。對于鋁合金材料的本構(gòu)特性主要通過元件級材料試驗(yàn)獲得［5－6］，中應(yīng)變率以下材料力學(xué)試驗(yàn)一般用液壓伺服拉伸試驗(yàn)機(jī);高、超高應(yīng)變率試驗(yàn)則采用分離式霍普金森桿裝置。連接件的斷裂失效試驗(yàn)、鋁合金材料力學(xué)性能試驗(yàn)類似，主要根據(jù)不同加載速度，選液壓拉伸試驗(yàn)機(jī)或霍普金森桿裝置。設(shè)計(jì)合適的夾具，通過拉伸試驗(yàn)，獲得連接件的極限拉伸載荷及剪切載荷［7－9］。

本文針對飛機(jī)平尾內(nèi)段前緣所用鋁合金材料，用電子液壓伺服拉伸試驗(yàn)機(jī)進(jìn)行準(zhǔn)靜態(tài)及中應(yīng)變率下材料拉伸試驗(yàn);利用分離式霍普金森拉桿(SHTB)裝置獲得鋁合金在高應(yīng)變率下應(yīng)力－應(yīng)變關(guān)系，擬合出反映鋁合金不同應(yīng)變率本構(gòu)關(guān)系的Johnson-Cook模型方程。通過SHTB拉伸試驗(yàn)給出平尾結(jié)構(gòu)連接所用鉚釘?shù)臉O限拉伸及剪切載荷。基于PAM-CRASH軟件，結(jié)合鋁合金的Johnson-Cook方程及連接件斷裂載荷，耦合SPH方法與有限元方法建立民機(jī)平尾前緣的鳥撞數(shù)值計(jì)算模型，并通過試驗(yàn)對數(shù)值分析結(jié)果進(jìn)行驗(yàn)證。

1 鋁合金動(dòng)態(tài)力學(xué)性能測試及連接件失效試驗(yàn)

1.1 鋁合金動(dòng)態(tài)力學(xué)性能試驗(yàn)

飛機(jī)平尾內(nèi)段包括前緣組件(前緣蒙皮、輔助隔板組件)、前后梁組件及盒段肋組件。由于鳥撞過程中承受撞擊載荷的主要部位為前緣組件、肋組件及前梁組件，材料為2024－T3及7075－T6鋁合金，故對兩種材料進(jìn)行應(yīng)變率相關(guān)的材料動(dòng)態(tài)力學(xué)性能試驗(yàn)。

利用電子萬能試驗(yàn)機(jī)進(jìn)行鋁合金準(zhǔn)靜態(tài)拉伸試驗(yàn)，高速拉伸試驗(yàn)機(jī)進(jìn)行材料中應(yīng)變率拉伸試驗(yàn)，用SHTB裝置獲得鋁合金高應(yīng)變率下力學(xué)特性。不考慮外界溫度變化對材料力學(xué)特性影響，所有試驗(yàn)均在室溫下進(jìn)行。將試驗(yàn)所得材料工程應(yīng)力σE、應(yīng)變?chǔ)臙轉(zhuǎn)化為真實(shí)應(yīng)力 σT、應(yīng)變 εT，即

三種動(dòng)態(tài)拉伸試驗(yàn)對應(yīng)的應(yīng)變率見表1。

表1 動(dòng)態(tài)拉伸試驗(yàn)對應(yīng)的應(yīng)變率Tab.1 Strain rate for differentdynamic tensile tests

圖1為采用CSS－88000電子萬能試驗(yàn)機(jī)進(jìn)行的準(zhǔn)靜態(tài)拉伸試驗(yàn)，板材試件通過楔形夾具夾持在萬能試驗(yàn)機(jī)上，利用引伸計(jì)測量試件標(biāo)距段應(yīng)變。

圖1 準(zhǔn)靜態(tài)拉伸試驗(yàn)Fig.1 Quasi－ static tensile test

鋁合金材料中應(yīng)變率拉伸試驗(yàn)通過Instron VzHS 160/100－20高速液壓伺服試驗(yàn)機(jī)配合高速相機(jī)，采用非接觸測量方法完成，見圖2。在試件標(biāo)距段隨機(jī)噴灑非均勻散斑，見圖3，應(yīng)變采用數(shù)字圖像(Digital Image Correlation，DIC)技術(shù)測量。由Photron－1高速相機(jī)對高速拉伸階段進(jìn)行高頻采集，獲得每個(gè)負(fù)載階段標(biāo)距段表面圖像，通過分析軟件測量獲得標(biāo)距段在每個(gè)時(shí)間步晶格片的坐標(biāo)值，進(jìn)而獲得標(biāo)距段全場應(yīng)變。

圖2 中應(yīng)變率拉伸試驗(yàn)Fig.2 Tensile test under medium rate

圖3 試驗(yàn)件標(biāo)距段散斑Fig.3 Speckles on gauge length of the specimen

高應(yīng)變率拉伸試驗(yàn)借助SHTB裝置。由于材料供應(yīng)商未能提供合適的鋁合金棒材試件，只能通過設(shè)計(jì)專用接頭將板材試件安裝于加載桿，見圖 4。為保證試件拉伸過程保持一維狀態(tài)，其外形參照準(zhǔn)靜態(tài)試驗(yàn)標(biāo)準(zhǔn)設(shè)計(jì)。加工中嚴(yán)格控制尺寸精度、平行度，并反復(fù)調(diào)試、校核試件與霍普金森桿同軸度。工程應(yīng)力、應(yīng)變及應(yīng)變率基于一維應(yīng)力波理論分別求得，即

圖4 SHTB拉伸試驗(yàn)Fig.4 SHTB tensile test

式中:E為加載桿彈性模量;A，As分別為加載桿、試件截面積;L為試件標(biāo)距段長度;C0為波速;εI(t)，εR(t)，εT(t)分別為入射應(yīng)變、反射應(yīng)變及透射應(yīng)變。處理試驗(yàn)數(shù)據(jù)時(shí)，認(rèn)為應(yīng)力波在試件中來回傳播3個(gè)循環(huán)后應(yīng)力應(yīng)變達(dá)到均勻狀態(tài)，取此刻后的應(yīng)力應(yīng)變數(shù)據(jù)作為有效數(shù)據(jù)。

示波器記錄的2024－T3、7075－T6鋁合金試驗(yàn)典型原始波形見圖5，不同應(yīng)變率下應(yīng)力－應(yīng)變曲線見圖6。由兩圖看出，兩種鋁合金力學(xué)特性隨應(yīng)變率變化規(guī)律基本一致，即在同一應(yīng)變率區(qū)間內(nèi)動(dòng)態(tài)力學(xué)性能相近。從準(zhǔn)靜態(tài) 0.001 s－1到中應(yīng)變率 10 s－1，材料的屈服及抗拉強(qiáng)度僅有小幅提高，表明兩種材料在準(zhǔn)靜態(tài)、中應(yīng)變率下的力學(xué)特性差別不大;進(jìn)入高應(yīng)變率區(qū)間，鋁合金應(yīng)變率強(qiáng)化效應(yīng)明顯體現(xiàn)，屈服、抗拉強(qiáng)度較準(zhǔn)靜態(tài)及中應(yīng)變率有大幅提升;相同應(yīng)變率下，7075－T6的屈服、抗拉強(qiáng)度均高于2024－T3。

圖5 SHTB動(dòng)態(tài)拉伸試驗(yàn)典型波形Fig.5 Typical experimental waves of SHTB dynamic tension test

圖6 鋁合金不同應(yīng)變率下應(yīng)力－應(yīng)變曲線Fig.6 Stress-strain curves of aluminum alloy under different strain rates

1.2 Johnson-Cook本構(gòu)模型擬合

Johnson-Cook模型能反映金屬等材料應(yīng)變硬化效應(yīng)、應(yīng)變率強(qiáng)化效應(yīng)及溫度軟化效應(yīng)［10］。工程上廣泛用于描述材料的力學(xué)特性，表達(dá)式為

不考慮材料溫度軟化效應(yīng)，只需擬合獲得Johnson－ Cook模型方程中 A，B，C，n參數(shù)。取0.001/s為參考應(yīng)變率。應(yīng)變率敏感系數(shù)C為常數(shù)時(shí)，擬合的Johnson－Cook模型方程往往不能精確反映所有應(yīng)變率下的力學(xué)特性［11］，修正C并擬合獲得兩種鋁合金Johnson－Cook方程表達(dá)式為

擬合后的Johnson-Cook模型方程計(jì)算所得應(yīng)力－應(yīng)變曲線與試驗(yàn)結(jié)果對比見圖7。由圖7看出，兩種材料擬合結(jié)果均與試驗(yàn)值吻合較好，說明Johnson-Cook方程能較準(zhǔn)確反映鋁合金不同應(yīng)變率下本構(gòu)特性。

圖7 Johnson-Cook模型擬合結(jié)果與試驗(yàn)對比Fig.7 Comparison between Johnson-Cook model results and experimental results

1.3 連接件動(dòng)態(tài)失效試驗(yàn)

飛機(jī)平尾連接件為鉚釘、螺栓，受撞擊載荷的失效形式有拉斷及剪斷。連接件模擬及失效判據(jù)選擇對整個(gè)結(jié)構(gòu)尤其撞擊區(qū)域的動(dòng)態(tài)響應(yīng)至關(guān)重要。利用體單元能模擬實(shí)際連接，但對含大量連接件結(jié)構(gòu)而言，建模工作量較繁重。目前對連接件處理主要采用共節(jié)點(diǎn)、MPC單元或特定連接單元(參考連接件靜態(tài)力學(xué)參數(shù))解決。共節(jié)點(diǎn)及MPC單元未考慮連接件失效;而特定連接單元能有效模擬真實(shí)情況，若采用動(dòng)態(tài)失效力學(xué)參數(shù)替代靜態(tài)力學(xué)參數(shù)作為失效判據(jù)，則可進(jìn)一步提高結(jié)構(gòu)撞擊區(qū)域的計(jì)算精度。

平尾前緣與前梁組件，肋組件與蒙皮等區(qū)域主要采用鋁合金鉚釘(MS20426AD3－5、MS20426AD5－5、MS20470AD5－5、MS20470AD6 －5、NAS1097AD4 －5、NAS1097AD5－5及NAS1097AD6－5等)連接。利用SHTB裝置進(jìn)行鉚釘?shù)膭?dòng)態(tài)拉伸、剪切試驗(yàn)。鉚釘連接的母材選擇高強(qiáng)度鋼，通過轉(zhuǎn)接頭安裝于加載桿上，試件安裝見圖8。鉚釘橫截面拉伸、剪切載荷表達(dá)式為

式中:FI(t)，F(xiàn)T(t)分別為入射桿上力、透射桿上力;F(t)為鉚釘橫截面平均拉伸載荷或剪切載荷。

據(jù)鳥撞試驗(yàn)撞擊區(qū)域鉚釘失效大致時(shí)間設(shè)置SHTB裝置中子彈入射速度，使鉚釘從加載到斷裂時(shí)間盡可能與鳥撞試驗(yàn)一致。由于連接件動(dòng)態(tài)失效的母材為高強(qiáng)度鋼而非真實(shí)結(jié)構(gòu)件中鋁合金材料，會(huì)使試驗(yàn)結(jié)果有一定誤差。經(jīng)多次預(yù)試驗(yàn)，最終取入射桿速度10 m/s。每種規(guī)格鉚釘拉伸、剪切試驗(yàn)均分別進(jìn)行至少5次，取平均值作為最終極限載荷。

圖8 鉚釘動(dòng)態(tài)失效試驗(yàn)Fig.8 Dynamic failure test of the rivets

MS20470AD6－5鉚釘拉伸、剪切的力－位移曲線見圖9。鉚釘拉伸試驗(yàn)中主要失效為釘帽拉脫，見圖10(a);剪切試驗(yàn)中斷面則基本與上下母材搭接面平行，見圖10(b)。鉚釘動(dòng)態(tài)極限拉伸、剪切載荷及供應(yīng)商提供的靜態(tài)極限載荷對比見表2。由表2看出，隨加載速度提高極限載荷小幅增加。

圖9 MS20470AD6－5鉚釘力－位移曲線Fig.9 Force-displacement curves of MS20470AD6 －5 rivet

表2 鉚釘極限拉伸及剪切載荷Tab.2 Limit tensile load and shear load of the rivets

圖10 鉚釘失效Fig.10 Failure of the rivets

2 飛機(jī)平尾前緣鳥撞數(shù)值計(jì)算及試驗(yàn)驗(yàn)證

2.1 數(shù)值計(jì)算

鳥體按適航要求取3.6kg，形狀為中間圓柱、兩端半球柱體，據(jù)鳥體密度950kg/m3及長徑比2計(jì)算獲得鳥體長度及直徑。用SPH單元模擬鳥體可避免大變形導(dǎo)致的網(wǎng)格畸變，更好模擬高速?zèng)_擊過程中的鳥體破碎飛濺。引入Mornaghan EOS狀態(tài)方程［12］定義鳥體本構(gòu)模型，即

式中:P0為初始壓力;ρ0為鳥體初始密度;B，γ為常數(shù)，由數(shù)值計(jì)算結(jié)合鳥撞試驗(yàn)優(yōu)化反演獲得，本文取B=128 ×106MPa，γ =7.98［13］。

內(nèi)段平尾前緣撞擊點(diǎn)位置見圖11，該處靠近固定前緣及可卸前緣對接處，屬鳥撞薄弱區(qū)域。

圖11 鳥體撞擊位置Fig.11 Bird impact location

飛機(jī)內(nèi)段平尾前緣鳥撞數(shù)值計(jì)算模型見圖12。鳥體單元數(shù)2440;內(nèi)段平尾結(jié)構(gòu)件均用Beltyschko－Tsay殼單元模擬，數(shù)量42641;本構(gòu)模型采用材料動(dòng)態(tài)力學(xué)性能試驗(yàn)所得Johnson－Cook方程，當(dāng)單元等效塑性應(yīng)變達(dá)到材料極限應(yīng)變值時(shí)單元失效被刪除。由于鳥撞過程中材料應(yīng)變率主要集中在中應(yīng)變率區(qū)間，取2024－ T3 及7075 －T6 鋁合金在 10 s－1、100 s－1及 500 s－1三組應(yīng)變率下極限拉伸應(yīng)變平均值作為失效應(yīng)變，兩種材料等效失效應(yīng)變分別為0.181及0.148。鉚釘用1923個(gè)Plink單元模擬，失效判據(jù)所需極限拉伸、剪切載荷則由連接件動(dòng)態(tài)失效試驗(yàn)(表2)提供。計(jì)算中拉伸或剪切載荷達(dá)到各自極限載荷時(shí)鉚釘失效。將翼根處盒段上下壁板邊緣節(jié)點(diǎn)固支約束模擬真實(shí)試件安裝夾。鳥體撞擊速度方向?yàn)轱w機(jī)逆航向，速度大小取試驗(yàn)中兩次激光測量平均值160.935 m/s。

圖12 鳥撞數(shù)值計(jì)算模型Fig.12 Numerical calculation model of bird impact

2.2 鳥撞試驗(yàn)

平尾鳥撞試驗(yàn)采用空氣炮法，將包扎好的鳥彈放入彈殼裝進(jìn)空氣炮管，啟動(dòng)空氣壓縮機(jī)，當(dāng)壓力達(dá)到所需值時(shí)打開空氣釋放機(jī)構(gòu)，在壓縮空氣作用下將鳥彈發(fā)射，并在炮口處將彈殼剝離，僅將鳥彈射出，由激光測量鳥彈飛行速度，使鳥彈按預(yù)定速度、方向撞擊試件指定部位。

試驗(yàn)設(shè)備包括空氣炮系統(tǒng)、控制系統(tǒng)、測速系統(tǒng)及高速攝像系統(tǒng)等。鳥撞試驗(yàn)設(shè)備框圖見圖13，空氣炮系統(tǒng)見圖14。測量設(shè)備包括高速激光測速儀、高速數(shù)碼攝像機(jī)、電子秤、手持式激光測距儀。

圖13 鳥撞試驗(yàn)裝置Fig.13 Bird impact test equipment

圖14 空氣炮系統(tǒng)Fig.14 Gas gun system

2.3 結(jié)果討論

本次飛機(jī)內(nèi)段平尾鳥撞試驗(yàn)主要據(jù)鳥撞適航條款規(guī)定，檢查平尾結(jié)構(gòu)是否能滿足適航要求。未測量鳥撞過程中結(jié)構(gòu)的應(yīng)變、位移。故計(jì)算結(jié)果與試驗(yàn)結(jié)果對比主要針對結(jié)構(gòu)的破壞形式，見圖15?？梢娪?jì)算所得撞擊點(diǎn)處總體破壞形式與試驗(yàn)結(jié)果基本吻合;試驗(yàn)中蒙皮沿展向變形340 mm，而數(shù)值計(jì)算展向變形為336 mm。不同時(shí)刻內(nèi)段平尾等效塑性應(yīng)變云圖見圖16。由圖16看出，前緣蒙皮等效塑性應(yīng)變超過失效應(yīng)變被擊穿，一部分鳥體進(jìn)入前緣內(nèi)部，撞擊肋板及輔助隔板，另一部分鳥體則緊貼前緣蒙皮沿翼展方向從翼根向翼梢擴(kuò)散。結(jié)構(gòu)破壞主要集中在撞擊點(diǎn)附近的前緣組件，連接帶板被撕裂，輔助隔板變形破壞嚴(yán)重，前梁腹板及緣條未發(fā)生可視變形及破壞，見圖17。

圖15 內(nèi)段平尾前緣破壞對比Fig.15 Comparison of damage forthe root horizontal tail wing leading edge

圖16 內(nèi)段平尾等效塑性應(yīng)變云圖Fig.16 Equivalence plastic strain fringe of the root horizontal tail wing

圖17 前緣內(nèi)部變形破壞情況Fig.17 Damage ofthe inner leading edge

檢查鳥撞試驗(yàn)后結(jié)構(gòu)連接發(fā)現(xiàn)，固定前緣蒙皮在與對接帶板連接的釘孔處撕裂，撕裂18個(gè)釘距;可卸前緣與對接板連接的釘孔拉脫10個(gè)×2排;4號肋與蒙皮有6個(gè)釘脫落。圖18為連接件破壞模擬分析結(jié)果，可見與試驗(yàn)結(jié)果基本相符。

圖18 連接件破壞模擬結(jié)果Fig.18 Simulation result of the failure for rivets

通過計(jì)算結(jié)果與試驗(yàn)結(jié)果對比，無論前緣組件、肋組件、輔助隔板變形及破壞，或連接件破壞等均吻合較好，表明所建數(shù)值計(jì)算模型可靠，可有效指導(dǎo)飛機(jī)平尾復(fù)雜結(jié)構(gòu)的抗鳥撞設(shè)計(jì)。

3 結(jié)論

基于積木式試驗(yàn)分析方法，進(jìn)行元件級材料力學(xué)性能試驗(yàn)、連接件動(dòng)態(tài)失效破壞試驗(yàn)及部件級結(jié)構(gòu)飛機(jī)平尾前緣鳥撞分析、試驗(yàn)驗(yàn)證，結(jié)論如下:

(1)平尾前緣所用2024－T3、7075－T6鋁合金材料在中應(yīng)變率至高應(yīng)變率區(qū)間，均表現(xiàn)出較明顯的應(yīng)變率強(qiáng)化效應(yīng)，在瞬態(tài)動(dòng)力學(xué)分析中須考慮。

(2)提高加載速度，7種鋁合金鉚釘?shù)臉O限剪切載荷僅有小幅增加，即極限載荷對加載速度不敏感。

(3)民機(jī)平尾前緣鳥撞數(shù)值結(jié)果與試驗(yàn)結(jié)果基本吻合，表明據(jù)鋁合金材料拉伸試驗(yàn)擬合的Johnson-Cook方程、連接件破壞試驗(yàn)建立的失效判據(jù)及所用鳥撞數(shù)值計(jì)算方法合理、可靠。

［1］Georgiadis S，Gunnion A J，Thomson R S，et al.Bird-strike simulation for certification of the Boeing 787 composite moveable trailing edge［J］.Composite Structures，2008，86(1/2/3):258－268.

［2］Guida M，Marulo F，Meo M，et al.SPH-Lagrangian study of bird impact on leading edge wing［J］.Composite Structures，2011，93(3):1060－1071.

［3］Heimbs S.Computational methods for bird strike simulations:a review［J］.Computers and Structures，2011，89(23/24):2093－2112.

［4］劉軍，李玉龍，劉元鏞.基于SPH方法的葉片鳥撞數(shù)值模擬研究［J］.振動(dòng)與沖擊，2008，27(9):90－93.LIU Jun，LI Yu-long，LIU Yuan-yong.Numerical simulation study of bird-impact on a blade using SPH method ［J］.Journal of Vibration and Shock，2008，27(9):90－93.

［5］李娜，李玉龍，郭偉國.3種鋁合金材料動(dòng)態(tài)性能極其溫度相關(guān)性對比研究［J］.航空學(xué)報(bào)，2008，29(4):903－908.LI Na， LI Yu-long， GUO Wei-guo. Comparison of mechanical properties and their temperature dependencies for three aluminium alloys under dynamic load ［J］.Acta Aeronautica ET Astronautica Sinica，2008，29(4):903－908.

［6］王宏欣，查曉熊.3004鋁的動(dòng)態(tài)力學(xué)性能及本構(gòu)模型［J］.華中科技大學(xué)學(xué)報(bào):自然科學(xué)，2011，39(5):39－42.WANG Hong-xin，ZHA Xiao-xiong.Dynamic mechanical behavior and constitutive model of 3004 aluminum ［J］.Journal of Huazhong University of Science and Technology:Natural Science Edition，2011，39(5):39 －42.

［7］Porcaro R， Langseth M， Hannssen A G， et al.Crashworthiness of self-piercing reveted connections［J］.International Journal of Impact Engineering，2008，35(11):1251－1266.

［8］Hoang N H，Hannssen A G，Langseth M，et al.Structural behaviour of aluminium self-piercing riveted joints:An experimental and numerical investigation［J］.International Journal of Solids and Structures，2012，49(23/24):3211－3223.

［9］Han L，Chrysanthou A，Young K W.Mechanical behaviour of self-piercing riveted multi-layer joints under different specimen configurations［J］.Materials ＆ Design，2007，28(7):2024－2033.

［10］Pothnis J R，Perla Y，Arya H，et al.High strain rate tensile behavior of aluminum alloy 7075 T651 and IS 2062 mild steel［J］.Journal of Engineering Materials and Technology，2011，133(2):1－9.

［11］林木森，龐寶君，張偉，等.5A06鋁合金的動(dòng)態(tài)本構(gòu)關(guān)系實(shí)驗(yàn)［J］.爆炸與沖擊，2009，29(3):306－311.LIN Mu-sen， PANG Bao-jun， ZHANG Wei， et al.Experimental investigation on a dynamic constitutive relationship of 5A06 Al alloy［J］.Explosion and Shock Waves，2009，29(3):306 －311.

［12］Monaghan J J，Kocharyan A.SPH simulation of multiphase flow［J］.Computer Physics Communication，1995，87(1/2):225－235.

［13］McCarthy M A，Xiao J R，McCarthy C T，et al.Modeling of bird strike on an aircraft wing leading edge made from fibre met al laminates part 2:modeling of strike with SPH bird model［J］.Applied Composite Materials，2004，11(5):317－340.