陳賀賀， 原梅妮， 李立州， 史明東， 何小晶， 韓玉杰， 姜 波,3

(1. 中北大學 機電工程學院， 山西 太原 030051； 2. 中國航空制造技術研究院， 北京 100024；3. 北京科技大學， 北京 100083)

0 引 言

飛鳥與飛機在空中相撞所產(chǎn)生的事故稱為“鳥撞”. 飛鳥的撞擊會給飛機造成嚴重后果， 甚至發(fā)生機毀人亡的事故. 風擋是飛機上重要且易遭鳥撞擊的部件之一， 其抗鳥撞研究一直是每個國家在飛機設計過程中一個必不可少的過程， 風擋的設計準則要求風擋在飛機巡航速度下能夠抵抗 1.8 kg 的飛鳥撞擊而不擊穿[1-4].

風擋鳥撞問題的研究一般采用試驗或數(shù)值模擬等方法. 朱書華[5]等通過某型飛機全尺寸風擋鳥撞試驗， 研究了鳥撞風擋的動響應過程， 所獲結果為建立準確的鳥撞風擋有限元模型提供了重要的試驗依據(jù)和驗證算例. 姚小虎[6]從試驗和數(shù)值模擬兩方面研究了鳥撞飛機風擋的動力響應問題， 數(shù)值模擬結果與試驗吻合較好， 驗證了有限元模型的正確性. 另外， 由于鳥撞試驗成本極為昂貴， 所以鳥撞問題的研究廣泛采用數(shù)值模擬方法. 目前國內(nèi)外鳥撞問題的仿真研究， 在分析方法上多運用拉格朗日(Lagrange)、 光滑粒子流體動力學(SPH)法定義鳥體模型. 王富生等[7]通過Lagrange方法建立塑性動力學鳥體模型， 采用多島遺傳算法等優(yōu)化方法實現(xiàn)了鳥體材料參數(shù)的反演. 賈建東等[8]采用SPH方法建立鳥體有限元模型， 對某飛機圓弧風擋鳥撞過程進行了數(shù)值模擬， 計算結果與試驗基本吻合； 同時與鳥體采用任意拉格朗日法(ALE)和無網(wǎng)格伽遼金方法(EFG)進行對比， 驗證了SPH方法在分析鳥撞問題中的優(yōu)越性. 朱書華等[9]還分別采用耦合解法和解耦解法研究了圓柱形和兩端半球形、 中間圓柱形的鳥體形狀對風擋鳥撞動響應的影響， 結果表明： 兩種形狀的鳥體模型計算結果與試驗結果都基本相符. 王猛等[10]對一種非對稱結構飛機前風擋的鳥撞動態(tài)響應進行了三維數(shù)值模擬， 計算結果表明這種非對稱結構設計并不會使鳥撞風擋臨界速度產(chǎn)生明顯的降低. Uzair Ahmed Dar等[11]還研究了不同撞速和不同撞擊角度下風擋的鳥撞動響應. 但在風擋鳥撞問題研究中， 影響因素眾多， 目前并未對全部影響因素一一展開研究.

本文采用有限元仿真方法開展了飛機風擋鳥撞問題的研究， 重點對不同鳥體力學模型、 風擋上不同撞擊位置以及風擋不同材料等因素進行了研究.

1 鳥撞飛機風擋有限元模型

1.1 Lagrange法和SPH法

在鳥撞問題仿真研究時， 鳥體模型一般采用Lagrange或SPH方法定義， 材料模型一般選擇彈塑性或自定義材料模型. 其中Lagrange[12]法多用于固體力學的有限元計算中， Lagrange網(wǎng)格將節(jié)點固定在分析對象上， 單元由節(jié)點連接形成， 并組成網(wǎng)格. 網(wǎng)格固定在物體上隨物體一起運動， 當分析對象發(fā)生變形時， 網(wǎng)格節(jié)點隨之移動， 并與物質(zhì)點始終保持重合， 同時， 單元也隨之變形. SPH是Lucy等人于1977年最早提出的一種無網(wǎng)格粒子法， 首先被應用于解決無邊界天體問題， 后逐漸在流體動力學、 侵徹、 碰撞等領域得到廣泛應用. SPH把分析對象離散化， 使用固定質(zhì)量的可動點即質(zhì)點或節(jié)點代替網(wǎng)格單元， 減少了有限元法中單元劃分的工作， 也沒網(wǎng)格畸變等問題.

1.2 鳥體有限元模型

鳥體幾何模型[7]采用中間圓柱、 兩端半球體的實體， 長徑比約為2∶1， 球體半徑為0.053 m， 中間圓柱長為0.141 8 m， 鳥體質(zhì)量為 1.8 kg. 采用Lagrange， SPH方法分別定義鳥體有限元模型， 在Lagrange方法中， 鳥體采用Solid164實體單元， 鳥體節(jié)點數(shù)為4 390， 單元數(shù)為3 480. 在SPH法中， 鳥體模型在后處理器LS-PrePost中被離散成光滑粒子， 粒子數(shù)為4 425. 兩種鳥體有限元模型， 如圖 1 所示.

在分析中， 彈塑性模型選用帶有失效模式的Plastic Kinematic材料模型. 該模型采用剪切失效準則， 當鳥體應變達1.25時， 單元失效， 失效單元將從網(wǎng)格中自動刪除. 材料參數(shù)值[7]見表 1 所示.

賈建東、 李冶方[8,13-14]等人采用自定義材料模型， 選用*MAT_NULL材料模型和*EOS_GRUNEISEN狀態(tài)方程定義鳥體材料模型， 研究表明仿真結果與實驗結果一致， 鳥體模型具體參數(shù)[14]， 如表 2 所示.

表 1 鳥體彈塑性模型參數(shù)Tab.1 Parameters of the bird model based on plastic kinematic

表 2 鳥體材料模型參數(shù)Tab.2 Parameters of the bird model based on hydrodynamic model

*EOS_GRUNEISEN狀態(tài)方程能用兩種方法確定壓力與體積之間的關系， 從而判斷出材料為壓縮性質(zhì)還是膨脹性質(zhì). 用來確定壓縮性質(zhì)材料的表達式為

(1)

用來確定膨脹性質(zhì)材料的表達式為

P=ρ0C2μ+(γ0+aμ)E,

(2)

式中：μ=ρ/ρ0-1稱為壓縮系數(shù)；ρ,ρ0為即時和初始鳥體密度；C,γ0,a,S1～S3為與材料沖擊壓縮特征有關的常數(shù)；E為內(nèi)能.

1.3 風擋有限元模型

飛機風擋模型[6,15]為單層等厚度圓弧風擋， 材料為DYB-3航空有機玻璃， 厚度為18 mm， 沿風擋表面對稱線前后跨度為800 mm， 圖 2 給出風擋模型及風擋表面對稱線上3個位置：A點為前1/3點， 距前緣266.7 mm；B點為中點， 距前緣400 mm；C點為后1/3點， 距前緣533 mm.

圖 2 飛機風擋模型及其對稱線上位置示意圖Fig.2 FE model of the windshield and the position in its symmetrical line

風擋模型采用Lagrange法定義， 并采用殼單元Shell163進行網(wǎng)格劃分， 劃分單元時將前后邊緣分成96份， 兩側分成80份， 風擋模型共包含7 856個節(jié)點和7 680個單元. 風擋邊界采用固支約束， 材料模型選為線彈性材料， 材料參數(shù)[8]如表 3 所示. 有限元模型中鳥體與風擋接觸方式選用基于罰函數(shù)的接觸算法.

表 3 風擋材料參數(shù)Tab.3 Parameters of the windshield material

2 鳥體力學模型研究

選擇不同材料模型和方法建立鳥體有限元模型， 研究不同鳥體模型對風擋遭撞擊的動響應影響. 為了敘述方便， 對有限元模型進行編號， 見表 4， 其中風擋模型均一致， 采用線彈性模型和Lagrange方法.

對有限元模型進行仿真計算， 4種鳥體模型分別以515 km/h的速度撞擊風擋對稱線前1/3點， 模擬時間均為15 ms. 表 4 同時給出了4種有限元模型與試驗[6]在4.5 ms時刻的撞擊瞬間圖， 可以看出， 此時鳥體模型滑移至風擋后部， 采用彈塑性模型的模型一和模型三鳥體變形相似， 鳥體未發(fā)生大變形； 采用自定義模型的模型二和模型四中鳥體模型表現(xiàn)出大變形狀態(tài)， 但是采用Lagrange方法的模型二出現(xiàn)網(wǎng)格嚴重畸變， 會對鳥撞后風擋響應數(shù)據(jù)的精度產(chǎn)生影響， 而模型四中鳥體被離散成粒子， 不存在網(wǎng)格畸變問題， 能較好呈現(xiàn)鳥體在高速沖擊下呈流體狀飛濺的狀態(tài)， 與風擋耦合效應明顯， 與試驗中鳥體變形模態(tài)更加接近.

表 4 4種鳥撞風擋有限元模型及4.5 ms撞擊瞬間圖Tab.4 Four FE models of the bird-impact windshield and the picture of impact in 4.5 ms

取4種鳥體模型中風擋對稱線上中點B的位移和應變數(shù)據(jù)與試驗實測數(shù)據(jù)[6]進行對比， 結果見表 5. 可以看出， 當選擇彈塑性模型時， Lagrange方法的結果最優(yōu)， 與試驗數(shù)據(jù)相比， 誤差在5%以內(nèi)； 當選擇自定義材料模型時， SPH方法的結果最優(yōu)， 與試驗數(shù)據(jù)相比， 平均誤差在5%以內(nèi). 所以， 模型一和模型四均能準確預測風擋的鳥撞動態(tài)響應.

表 5 計算結果對比Tab.5 Comparison of calculation results

綜上， 采用SPH方法和自定義材料的模型四可以模擬鳥體高速撞擊風擋呈流體狀飛濺的過程， 而且仿真結果接近試驗結果， 因此模型四最適合模擬鳥撞風擋過程， 接下來選擇模型四做進一步分析.

3 風擋不同撞擊點處鳥撞響應

對鳥撞過程進行仿真計算， 鳥體的撞速范圍是450～650 km/h， 對風擋的撞擊位置分別為風擋對稱線上前1/3點、 中點、 后1/3點， 對應圖 2 中的A，B，C點， 模擬時間均為15 ms.

3.1 不同速度下的撞擊力變化

圖 3 為鳥體以不同速度撞擊風擋A點時， 中間B點位置的撞擊力時程曲線. 從圖 3 中可以看出， 隨鳥體撞擊速度增大， 風擋所受撞擊力也增大， 而且達到峰值的時間提前， 速度650 km/h的撞擊力峰值比速度450 km/h的撞擊力峰值提前0.5ms左右， 并且從圖中可以看出整個鳥撞時間約4 ms. 從650 km/h的撞擊力曲線可以看出， 0.9～3.4 ms為鳥體剛撞上風擋， 速度急劇減小， 風擋撞擊力急劇增大， 在3.4 ms撞擊力達到最大值71.8 kN； 但由于鳥體流變， 風擋表面撞擊力又迅速衰減， 6.2 ms開始風擋撞擊力又重新增大， 在8.4 ms左右達到峰值后又開始衰減， 這是由于鳥撞載荷沖擊波回彈造成風擋所受撞擊力二次增大. 速度為550 km/h和450 km/h的撞擊力曲線變化趨勢也類似. 這說明， 風擋遭受鳥撞之后， 一方面撞擊力會迅速增大， 并向風擋周圍擴散， 另一方面應力波還會在風擋中來回振蕩傳遞， 從而對風擋不斷地產(chǎn)生作用.

圖 3 不同撞擊速度下風擋的撞擊力時程曲線Fig.3 Impact force history curves of windshield under different impact velocities

3.2 不同撞擊點處發(fā)生失效的臨界條件

有限元模型中風擋采用最大主應力失效模式， 當最大主應力達到失效強度78 MPa時[8]， 風擋單元失效. 仿真計算了風擋不同撞擊位置所發(fā)生失效破壞的臨界撞速及對應的臨界撞擊力， 見表 6.

表 6 不同撞擊點的臨界撞速和臨界撞擊力Tab.6 Critical impact velocities and force at different impact points

可以看出， 鳥撞風擋A點時， 撞擊速度超過595 km/h風擋會發(fā)生失效破損， 撞擊點為B點時， 風擋發(fā)生失效的臨界速度為545 km/h， 撞擊點為C點時， 臨界撞速為480 km/h. 在風擋的3個撞擊點處， 風擋發(fā)生失效破壞的臨界撞速和臨界撞擊力在前1/3點最大， 中點次之， 后1/3點最小.

3.3 不同撞擊點處的位移分析

圖 4 為鳥體以562 km/h的速度撞擊風擋A,B,C點時對稱線上最大位移隨風擋玻璃表面位置變化曲線.

圖 4 風擋對稱線上最大位移隨風擋表面位置變化的曲線Fig.4 The history curves of maximum displacement change with thewindshield surface position at windshield symmetrical line

當鳥體撞擊風擋A點(距前邊緣266.7 mm)時， 風擋位移沿風擋對稱線從前往后先增大， 在距前邊緣350 mm處達到最大值21.75 mm， 隨著鳥體動能的逐漸減小， 風擋位移沿對稱線向后逐漸減小， 由于風擋四邊固支約束， 所以前邊緣和后邊緣的位移為0. 撞擊點為B和C的曲線變化趨勢類似， 撞擊點為B點(距前邊緣400 mm)時， 風擋對稱線上最大位移為27.5 mm， 發(fā)生在距前邊緣500 mm處； 撞擊點為C點(距前邊緣533 mm)時， 風擋對稱線上最大位移為28.38 mm， 發(fā)生在距前邊緣590 mm處. 可以得出， 撞擊點從A點到C點， 越往后， 風擋的位移就越大， 且無論撞擊點在哪， 風擋的最大位移均未發(fā)生在撞擊點處， 而是在撞擊點后57～100 mm處. 風擋位移達到25 mm(約風擋厚度1.4倍)以上的區(qū)域在距前邊緣435～630 mm范圍內(nèi)， 說明風擋的最大變形發(fā)生在其中后部區(qū)域.

綜上， 在風擋3個撞擊點中， 前1/3點處抗鳥撞的臨界速度和臨界撞擊力最大， 沿風擋向后的兩個點依次減小； 撞擊點沿風擋對稱線向后， 風擋變形越大， 且最大變形發(fā)生在風擋中后部. 所以， 可以得到， 風擋前部抵抗飛鳥撞擊能力最強， 中后部較弱.

4 不同風擋材料抗鳥撞性能對比分析

YB-3和DYB-3航空有機玻璃(PMMA)是國內(nèi)飛機應用較廣泛的兩種飛機風擋材料[16]， 本節(jié)對兩種風擋材料進行鳥撞性能對比分析. YB-3 PMMA的材料參數(shù)詳見文獻[5].

表 7 為兩種材料風擋不同撞擊位置的臨界撞速， 可以看出YB-3 PMMA風擋發(fā)生失效的臨界撞速為310～385 km/h， 小于 DYB-3 PMMA風擋臨界撞速.

表 7 兩種材料風擋的臨界撞速Tab.7 Critical impact velocities of windshield by two materials symmetrical line

鳥體以350 km/h的速度分別撞擊DYB-3和YB-3 PMMA風擋A點時， 獲得了B點的位移時程曲線(見圖 5)， 以及B點的應變時程曲線(見圖 6).

圖 5 兩種材料風擋位移時程曲線Fig.5 Displacement history curves of windshield by two materials

由圖 5 可以看出， YB-3 PMMA的位移峰值為15.2 mm， 大于DYB-3 PMMA的最大位移9.7 mm， 但DYB-3 PMMA的位移時程曲線頻率大于YB-3 PMMA； 由圖 6 可以看出， YB-3 PMMA的應變峰值為10.5×10-3， 大于DYB-3 PMMA應變峰值近一倍. DYB-3 PMMA剛度相對較大、 變形較小， 從而風擋位移、 應變較小， 但遭鳥撞后發(fā)生振蕩的頻率較大.

圖 6 兩種材料風擋應變時程曲線Fig.6 Strain history curves of windshield by two materials

綜上可以看出， DYB-3 PMMA風擋在鳥撞仿真中表現(xiàn)出較高的耐沖擊性， 而且DYB-3 PMMA是由YB-3板定向拉伸而成， 分子鏈經(jīng)過拉伸取向后， 使其具備更好的抗銀紋和抗裂紋擴展等性能[16]. 所以， DYB-3 PMMA更適合用在氣動載荷較大的部位或者巡航速度較大的飛機風擋.

5 結 論

1)采用SPH方法和自定義材料的鳥體模型可以模擬鳥體高速撞擊風擋呈流體狀飛濺的過程， 并能準確預測風擋的鳥撞動態(tài)響應.

2)風擋所受撞擊力隨鳥體撞擊速度增大而增大， 且達到峰值的時間提前. 飛鳥撞擊風擋對稱線前1/3點時， 風擋發(fā)生失效破壞的臨界撞速和臨界撞擊力較大， 中點次之， 后1/3點最小.

3)撞擊點從前1/3點到后1/3點， 越往后， 風擋鳥撞后的位移就越大； 無論撞擊點在哪， 風擋的最大位移均未發(fā)生在撞擊點處， 而是在撞擊點后57～100 mm處； 風擋位移達到25 mm(約風擋厚度1.4倍)以上的區(qū)域在距前邊緣435～630 mm范圍內(nèi).

4) DYB-3 PMMA風擋在鳥撞仿真中表現(xiàn)出較高的耐沖擊性， 適合用在氣動載荷較大的部位或者巡航速度較大的飛機風擋.

參考文獻：

[1] 李玉龍, 石霄鵬. 民用飛機鳥撞研究現(xiàn)狀[J]. 航空學報, 2012, 33(2)： 189-198.

Li Yulong, Shi Xiaopeng. Investigation of the present status of research on bird impacting on commercial airplanes[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(2)： 189-198. (in Chinese)

[2] Mccarty R, Gran M, Baruch M. MAGNA nonlinear finite element analysis of T-46 aircraft windshield bird impact[C]. Aircraft Systems, Design and Technology Meeting, 2013: 263.

[3] 王富生, 岳珠峰, 馮震宙, 等. 鳥撞飛機風擋動態(tài)響應數(shù)值模擬方法研究現(xiàn)狀[J]. 飛機設計, 2008, 28(5)： 39-46.

Wang Fusheng, Yue Zhufeng, Feng Zhenzhou, et al. Present study of numerical simulation methods of dynamic response for aircraft windshield under bird strike[J]. Aircraft Design, 2008, 28(5)： 39-46. (in Chinese)

[4] Mohagheghian I, Wang Y, Zhou J, et al. Deformation and damage mechanisms of laminated glass windows subjected to high velocity soft impact[J]. International Journal of Solids & Structures, 2017(109)： 46-62.

[5] 朱書華. 鳥撞飛機風擋動響應分析與仿真試驗平臺研究[D]. 南京： 南京航空航天大學, 2009.

[6] 姚小虎. 鳥撞飛機圓弧風擋的實驗研究及數(shù)值模擬[D]. 太原： 太原理工大學, 2001.

[7] 王富生, 李立州, 王新軍， 等. 鳥體材料參數(shù)的一種反演方法[J]. 航空學報, 2007, 28(2)： 344-347.

Wang Fusheng, Li Lizhou, Wang Xinjun, et al. A method to identify bird’s material parameters[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(2)： 344-347. (in Chinese)

[8] 賈建東, 李志強, 楊建林, 等. 用SPH和有限元方法研究鳥撞飛機風擋問題[J]. 航空學報, 2010, 31(1)： 136-142.

Jia Jiandong, Li Zhiqiang, Yang Jianlin, et al. A study of bird impact on aircraft windshield using SPH and finite element method[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(1)： 136-142. (in Chinese)

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

Zhu Shuhua, Tong Mingbo. Study on dynamic responses of bird striking an aircraft windshield and virtual test platform[J]. Nanjing University of Aeronautics and Astronautics, 2009： 16-86. (in Chinese)

[10] 王猛, 黃德武, 羅榮梅. 飛機前風擋非對稱結構的鳥撞數(shù)值模擬[J]. 機械科學與技術, 2012, 31(2)： 291-294.

Wang Meng, Huang Dewu, Luo Rongmei. Simulating aircraft’s asymmetrical windshield subjected to bird impact[J]. Mechanical Science and Technology for Aerospace Engineering, 2012, 31(2)： 291-294. (in Chinese)

[11] Dar U A, Zhang W, Xu Y. FE analysis of dynamic response of aircraft windshield against bird impact[J]. International Journal of Aerospace Engineering, 2013(4)： 1-12.

[12] 李利莎, 謝清糧, 鄭全平, 等. 基于Lagrange、 ALE和SPH算法的接觸爆炸模擬計算[J]. 爆破, 2013, 28(1)： 18-22.

Li Lisha, Xie Qingliang, Zheng Quanping, et al. Numerical simulation of contact explosion based on Lagrange ALE and SPH[J]. Blasting, 2013, 28(1)： 18-22. (in Chinese)

[13] 張志林, 張啟橋, 李銘興. 飛機圓弧風擋鳥撞動響應分析[J]. 航空學報, 1992, 13(9)： A538-A542.

Zhang Zhilin, Zhang Qiqiao, Li Mingxing. Bird impact dynamic response analysis for aircraft arc windshield[J]. Acta Aeronautica et Astronautica Sinica, 1992, 13(9)： A538-A542. (in Chinese)

[14] 李冶方. 風擋鳥撞模型的有限元分析[D]. 鄭州： 鄭州大學, 2013.

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

[16] 秦瑞祥, 歐迎春, 張保軍, 等. YB-3與DYB-3航空有機玻璃綜合性能對比[J]. 中國建材科技, 2011(1)： 29-31.

Qin Ruixiang, Ou Yingchun, Zhang Baojun, et al. The contrast of competitive properties between YB-3 and DYB-3 aviation sheet[J]. China Building Materials Science & Technology, 2011(1)： 29-31. (in Chinese)