劉信超， 徐亞芳， 王露晨， 陸曉華， 左洪福

(1. 南京航空航天大學 民航學院，南京 211100； 2. 中國民航上海航空器適航審定中心，上海 200335)

作為迎風部件，風擋與尾翼、機翼、縫翼、發(fā)動機唇口和吊掛等飛機部位一樣，容易遭受飛鳥撞擊。據(jù)統(tǒng)計，全世界每年大約發(fā)生1萬次鳥撞飛機事件，國際航空聯(lián)合會已把鳥害升級為“A”類航空災難[1]。根據(jù)CCAR 25.775條規(guī)定，位于正常執(zhí)行職責的駕駛員正前方的風擋玻璃及其支承結構，必須能經受住1.8 kg(4磅)的飛鳥撞擊而不被擊穿，此時飛機的速度(沿飛機航跡相對于飛鳥)等于按第25.335(a)條選定的海平面VC值。因此，運輸類飛機風擋必須經過地面試驗或仿真驗證，從而達到抗鳥撞損傷的設計要求。

根據(jù)適航審定經驗，鳥撞適航審定以試驗驗證為主。但由于地面鳥撞試驗的高昂費用，試驗次數(shù)十分有限。因此有必要在適航驗證之前，通過仿真研究手段對風擋材料、風擋厚度、自由螺栓邊和夾緊邊緣的不同安裝方式、撞擊姿態(tài)、撞擊位置以及撞擊速度等參數(shù)進行參數(shù)化研究，從而對試驗進行指導。Doubrava等[2]研究了總厚度為14、18、20 mm的多層風擋在1.8 kg鳥體以83-125m/s速度撞擊時的動態(tài)響應，發(fā)現(xiàn)風擋的臨界穿透速度與風擋厚度成正比。Kangas等[3]對不同構型以及不同材料的風擋進行了研究，發(fā)現(xiàn)對多層風擋影響最大的因素是聚合物層的厚度。Dar等[4]對PMMA材質的風擋進行了鳥體重量、鳥體形狀、速度、撞擊角度、撞擊位置的參數(shù)化研究。Wang等[5]對PMMA材質的風擋進行了環(huán)境溫度、撞擊位置以及撞擊速度的影響研究。Mohagheghian等[6]對SGP、PVB、PU層合風擋進行了構型、溫度以及厚度等參數(shù)的分析研究。

以上的參數(shù)化研究均針對有機玻璃風擋，有關無機玻璃風擋的參數(shù)化研究較少，而運輸類飛機的風擋通常為無機玻璃，因此需要加強對無機玻璃風擋的研究。Grimaldi等[7]對無機玻璃風擋對于撞擊角度、撞擊目標尺寸以及撞擊目標尺寸的敏感性進行了參數(shù)化研究，但撞擊目標比較簡化并忽略了風擋邊界的影響。為了充分考察風擋結構的抗鳥撞能力，本文使用真實機頭模型以及無機玻璃風擋模型，對撞擊位置進行參數(shù)化研究，旨在為適航驗證試驗中的撞擊點選取工作提供參考。并從損傷、接觸力、位移、應力以及能量變化等方面來綜合比較不同撞擊點的動態(tài)響應，確定不同撞擊點的嚴酷性，并最終為適航工作提供了撞擊參考點。

1 鳥體計算模型

1.1 鳥體建模

本文所建的鳥體幾何模型為兩端帶半球帽的圓柱體，縱向總長度為截面直徑D的兩倍，如圖1、2所示。具體建模過程為：首先在CATIA中建立兩端帶半球帽的圓柱體的幾何體，隨后導入HyperMesh環(huán)境中用Hexa單元離散化鳥體得到圖1，最終將圖1鳥體導入PAM-CARSH中利用Convert功能轉化為SPH模型，如圖2。圖1、2中的鳥體模型重1.8 kg，尺寸D為0.114 m。

圖1 鳥體Hexa單元網格模型

圖2 SPH鳥體計算模型

1.2 SPH算法及參數(shù)

由于鳥體在高速運動下展現(xiàn)出良好的流體性，故鳥體計算模型一般采用較為成熟的光滑粒子動力學(Smoothed Particle Hydrodynamics，SPH)方法進行建模。SPH算法最初是用在天體物理學領域，是一種插值方法[8]。Liu等[9]以及Monaghan[10]對SPH方法進行了詳細的總結。在SPH算法中，計算鳥體被離散成有限個粒子(即插值點)，各粒子均攜帶速度、密度及應力等物理特性，并根據(jù)控制方程的定義隨著其他粒子一起運動。每個粒子的物理特性均可以通過對相鄰粒子的物理特性進行插值計算得到，可以用規(guī)則的內插函數(shù)計算全部粒子的場函數(shù)值，從而近似描述整個問題的場分布[11]。對于流體動力學問題，相關的控制方程包括連續(xù)性方程(即質量守恒方程)與動量方程(即(Navier-Stokes方程)[12]，兩個方程的Lagrangian形式見式(1)、(2)。

連續(xù)性方程

(1)

動量方程

(2)

Tait狀態(tài)方程

(3)

忽略大氣壓力并將其變成Tait方程中的壓強常量p0，因此Tait方程變成Monaghan狀態(tài)方程[14]

(4)

式中：p為現(xiàn)時壓強；ρ/ρ0為鳥體當前密度與初始密度比值；b為體積彈性模量；γ為指數(shù)，文獻[15]中通過平板試驗擬合了b、γ，分別為2.8 GPa、7.99。初始密度ρ0為950 kg/m3。Tait方程的使用對象是弱可壓流體，因此b、γ的取值必須保證密度相對變化值|Δρ|/ρ0較小，一般要滿足|Δρ|/ρ0≤0.01。同時可知流體速度vf、流體中的聲速cs滿足式(5)

(5)

因此|vf|/cs需滿足|vf|/cs≤0.1，另外，已知體積彈性模量b、指數(shù)γ、流體中的聲速cs以及初始密度ρ0滿足式(6)關系

(6)

且本文中的最大流體速度為150 m/s，最終將上述取值代入式(5)、(6)，得到|vf|/cs≈0.087<0.1，說明文獻[15]中通過平板試驗反演所得鳥體參數(shù)的取值滿足Monaghan狀態(tài)方程的使用條件。

1.3 SPH算法驗證

本小節(jié)給出了依據(jù)上述算法所建立的SPH鳥體撞擊剛體的數(shù)值計算動態(tài)響應結果，通過與Wilbeck的試驗結果[16]、Zhang等[17]以及Hedayati等[18]的仿真結果進行對比，證明本文SPH算法及其參數(shù)的合理性。Zhang等以及Hedayati等的仿真工況條件與Wilbeck的試驗條件一致，在Wilbeck的試驗中，鳥體重量為1 kg，撞擊速度為116 m/s。兩端帶半球帽的圓柱體1 kg重鳥體仿真模型的尺寸D為0.093 m。因此本小節(jié)仿真對鳥體進行了重新建模，鳥體尺寸D變?yōu)?.093 m，密度保持不變?yōu)?50 kg/m3，重量變?yōu)? kg，SPH粒子數(shù)為39 456。經過計算，1 kg重的SPH鳥體撞擊剛體過程見圖3，撞擊過程維持約2 ms，本節(jié)鳥體動態(tài)變形過程與Zhang等的結果基本一致。圖4為沖擊壓強時間歷程曲線的對比情況，沖擊壓強曲線共分為沖擊峰及穩(wěn)流區(qū)兩個區(qū)域。本文SPH算法對沖擊壓強的沖擊峰預測較準，沖擊峰的上升趨勢及峰值與Wilbeck的試驗結果基本一致；雖然較試驗結果而言穩(wěn)流區(qū)過于平滑，但穩(wěn)流區(qū)的總體趨勢的預測較準，見圖4。綜上所述，本文的SPH算法及參數(shù)是比較合理的。

2 機頭計算模型

2.1 機頭建模

綜合考慮機頭結構件的尺寸比、計算精度以及計算效率，用殼單元離散化機頭金屬鈑金件，用實體單元離散化風擋各層結構。在CATIA環(huán)境中提取幾何模型的外形特征，在HyperMesh環(huán)境中修理模型細節(jié)，并將模型導入PAM-CRASH后進行計算模型前處理。所有殼單元均為4節(jié)點縮減積分單元，殼單元厚度方向設置5個積分點，單個單元邊長約14 mm，總單元數(shù)為357 599。所有實體單元均為八節(jié)點縮減積分單元，單個單元邊長約7 mm，總單元數(shù)為310 100。所有鉚釘用Plink單元模擬，總單元數(shù)為8 226。機頭結構計算模型及尺寸如圖5所示，風擋計算模型如圖6所示，風擋各層由外至內為：3 mm厚無機玻璃層、4 mm厚PU膠層、8 mm厚無機玻璃層、1.5 mm厚PVB膠層、6 mm厚無機玻璃層，其中每層無機玻璃的上下表面均有0.5 mm厚的鋼化處理層。模型的邊界條件為圖5(a)中機頭側壁板的所有緣條結構受固支約束。

(a) t=0 ms(b) t=0.7 ms(c) t=1 ms(d) t=2 ms

圖3 剛體撞擊驗證(1 kg，116 m/s)

Fig.3 Verification by impacting rigid target(1 kg，116 m/s)

圖4 沖擊壓強時間歷程曲線對比

2.2 結構材料及屬性

選取PAM-CRASH中102號彈塑性材料模型來描述機頭金屬結構件的力學行為，其塑性響應用可以反映金屬等材料應變硬化效應、應變率強化效應及溫度軟化效應的Johnson-Cook本構模型描述[19]

(a) 機頭整體模型

(b) 連接件模型(Plink單元)

圖6 風擋結構計算模型

σe=

(7)

圖5(a)中所示的座艙蓋頂部板、側壁板材料為2524鋁合金，風擋下部板材料為2024鋁合金，緣條、鈑金件的材料為7075、7050鋁合金。2024、7075鋁合金材料在中應變率至高應變率區(qū)間，均表現(xiàn)出較明顯的應變率強化效應，在瞬態(tài)動力學中須考慮[20]，同時考慮金屬材料的應變硬化效應，但忽略溫度軟化效應。材料的失效行為由等效失效塑性應變εfail進行判定，即當材料的等效塑性應變到達εfail時，認為材料失效，刪除其單元。文獻[21-24]對上述金屬材料的JC本構模型參數(shù)進行了測量，具體數(shù)值如表1所示。

表1 Johnson-Cook本構模型參數(shù)

Plink單元的失效準則為

(8)

式中：N與T分別為拉力和剪切力；Tmax與Smax分別為最大拉力、最大剪切力。本文中所有Plink單元的最大拉力、最大剪切力設為5.1 kN、3.2 kN，指數(shù)常量m、n分別設為1.5和2.1。

文獻[25]中對航空無機玻璃力學性能進行了試驗研究，如圖7所示，玻璃材料為典型的脆性材料；應力-應變曲線由彈性加載段和失效段組成，近似為一條曲線；隨著應力水平的提高，試樣斷裂，應力水平急劇下降。因此本文選取PAM-CRASH中16號彈塑性本構來描述其力學性能，其中失效刪除塑性應變取0.001。由圖7可知，無機玻璃的強度對應變率較為敏感，400 s-1應變率下的強度(1.1 GPa)幾乎為準靜態(tài)時強度(0.517 GPa)的兩倍。同時，本文考慮了無機玻璃表面化學鋼化處理的影響。文獻[26]中表明，視窗玻璃經過化學鋼化后，表面離子交換層的準靜態(tài)強度可達到0.8 GPa。由于缺乏表面離子交換層的高應變率下力學性能數(shù)據(jù)，本文較為保守地采用準靜態(tài)下的玻璃強度，其中中間未鋼化層強度取0.517 GPa，表面鋼化層強度取0.8 GPa。文獻[27-28]中給出了PU膠層以及PVB膠層的動態(tài)應力-應變曲線，如圖8所示，本文同樣采用選取PAM-CRASH中16號彈塑性材料模型近似地描述其力學行為。

3 撞擊位置選取

圖9為本仿真中選取的各個撞擊點，其中A1～A11位于風擋表面，B1～B3位于風擋中央加持結構表面。A1～A11全部位于航向右風擋，由于飛機機身具有對稱性，因此不對航向左風擋進行重復計算分析。本文選取的A1～A11、B1～B3撞擊點所包含的撞擊區(qū)域基本覆蓋了航向右風擋全部表面區(qū)域及中央夾持部件區(qū)域，其中A10、A11位于風擋邊角點，A5位于風擋中心點。本文對上述撞擊點依次進行鳥撞仿真計算，計算結果作為后續(xù)分析的輸入。上述所有撞擊工況中，鳥體撞擊方向均沿著機身方向，鳥體姿態(tài)無偏航無俯仰，撞擊速度均為150 m/s。

圖7 風擋無機玻璃應力應變曲線 Fig.7 Stress-strain curve of windshield unorganic glass(a) PU應力應變曲線(b) PVB應力應變曲線圖8 風擋膠層應力應變曲線Fig.8 Stress-strain curve of PU and PVB

圖9 撞擊點位置示意圖

4 仿真結果分析

從損傷結果對比各工況，發(fā)現(xiàn)仿真中僅在A10點對應工況下出現(xiàn)風擋的稍許單元刪除，且單元刪除僅存在于最外層無機玻璃的表面層，其它層均無損傷。本節(jié)以A10對應工況為例分析鳥撞模擬過程：鳥體以150 m/s的速度飛向預定著彈點，鳥體在接觸風擋玻璃的瞬間產生高沖擊壓力，撞擊點處的個別風擋單元產生應力集中現(xiàn)象并發(fā)生刪除，如0.16 ms時，接觸部位某單元的范式等效應力達到759 MPa，而后該單元刪除。同時撞擊點處出現(xiàn)凹陷，隨著時間的推移，凹陷區(qū)域經歷了加深、轉移、回彈現(xiàn)象，當接觸時間達到3 ms后，鳥體作用力所剩無幾并隨后沿風擋表面逐漸滑出，圖10給出了t=0 ms到t=3 ms過程中風擋與鳥體的耦合作用過程以及層合風擋的受力情況。

(a) t=0.4 ms

(b) t=0.8 ms

(c) t=1.6 ms

(d) t=3.0 ms

圖10 鳥撞過程及風擋Von Mises等效應力云圖(150 m/s)

Fig.10 Bird-strike process and Von Mises stress contour(150 m/s)

除了損傷情況，本文也從接觸力、動能損失程度、風擋結構變形、窗框結構受力等多個角度來比較工況嚴酷程度。圖11為接觸力時間歷程曲線，圖12為鳥體動能損失曲線。由圖11、圖12可知，各撞擊點對應的接觸力及動能變化趨勢基本一致，但接觸力峰值與剩余動能的大小顯著不同。以4 ms時刻的鳥體動能為剩余動能，并提取了各撞擊點對應的剩余動能及接觸力峰值，得到了圖13。

圖11 接觸力時間歷程曲線

圖12 鳥體動能損失曲線

圖13 不同撞擊位置時的動能損失與接觸力對比

Fig.13 Comparison between different impact positions about kinetic energy loss and contact force

根據(jù)仿真結果可知，如圖13箭頭所示，鳥體越靠近飛機機身縱軸，接觸力峰值越大；鳥體越靠近飛機機身縱軸，鳥體動能損失更大。A3、A6、A9、A10、A11、B1、B2、B3各工況下，鳥體均在縱軸附近，因此接觸力均表現(xiàn)為較大值。其次，根據(jù)對比分析發(fā)現(xiàn)，撞擊點的豎直位置對接觸力與鳥體動能損失的影響不明顯，如B1、B2、B3各工況下接觸力與動能損傷的對比情況。

圖14為鳥體動能變化與風擋變形情況對應關系，其中風擋變形情況由風擋上某一特殊點的位移來表征。該點取為對應工況下接觸力峰值出現(xiàn)時風擋背面上位移量最大點，本文以該點代表對應工況下風擋的變形情況。由圖14可知，當撞擊點距離風擋上下夾持邊界較遠時，即A3、A4、A5三個撞擊點，風擋的變形量較大，而其他工況下由于風擋受到邊界夾持部件的影響，風擋變形情況偏小。層合風擋通過變形耗散部分鳥體動能，中心撞擊點A5撞擊工況下層合風擋吸收能量最多并達到342.7 J，邊角撞擊點A7撞擊工況下層合風擋吸能最少且僅有255.08 J，表明當著彈點位于窗框邊界附近時，窗框會通過自身支持剛度減少風擋變形吸能。然而窗框結構的支持剛度對各撞擊點的影響不同，如位于下邊角的A10工況下層合風擋吸能為263.80 J，而位于風擋上邊角的A11工況在接觸力與A10相當且均位于邊角的情況下，風擋吸能較大且為302 J，說明上下側窗框的支持剛度存在差異，會影響風擋的吸能過程。

圖14 鳥體動能變化與風擋變形情況對應關系

除了風擋內能外，窗框結構支持剛度對能量分配的影響也體現(xiàn)在自身內能變化過程中，窗框結構也通過變形行為耗散部分鳥體動能?？紤]到各部位支持剛度的不同，圖15中將窗框結構分為上側、下側、外側以及中央窗框四個部位，圖16中對上述窗框各個部位出現(xiàn)的內能峰值進行了統(tǒng)計。由于A1、A4、A7各工況的接觸力較小(見圖13)，且處于風擋外側的位置，鳥體易滑走，威脅較小，因而圖16中未包含該三點。此外，由于窗框各部位與著彈點的距離不同而導致應力波到達各部位的時間不同，故每一個撞擊工況下都分別單獨考慮四個部位的內能峰值大小，即忽略各部件的峰值到達時間。

圖15 窗框結構

圖16中的圓圈表示對應結構出現(xiàn)了較大的內能以及較大的變形，其中A系列撞擊點中，A2對上側窗框影響最大，A9對下側窗框影響最大，A11對中央窗框影響最大；在B系列撞擊點中，B2對中央撞擊點影響最大；所有撞擊點對外側窗框的影響均明顯小于對其他窗框結構的影響，因此忽略對外側窗框的分析。經過應力分析發(fā)現(xiàn)，這些對應結構出現(xiàn)應力集中現(xiàn)象：在A2撞擊工況中，上側窗框Von Mises應力在1.12 ms時達到544 MPa；在A9撞擊工況下，下側窗框Von Mises應力在1.6 ms時達到517 MPa；在A11撞擊工況下，中央窗框Von Mises應力在2.08 ms時達到329 MPa；在B2撞擊工況下，中央窗框Von Mises應力在2.32 ms時達到349 MPa。以A2工況為例，給出上側窗框結構的Von Mises等效應力云圖，如圖17所示，其集中應力已經超過材料的屈服強度473 MPa，局部產生塑性變形。

圖16 窗框結構內能峰值分布

圖17 A2工況下上側窗框結構Von Mises等效應力云圖(背面)

綜上所述，為了充分驗證風擋的抗鳥撞安全性，必須考慮到鳥體距離機身縱軸遠近、邊界夾持部件兩大因素的影響。因此建議在驗證試驗中選取A2、A9、A11、B2四點作為試驗靶點。

5 結 論

本文以飛機無機玻璃層合風擋為研究對象，研究了無機玻璃層合風擋的鳥撞動態(tài)響應特點及撞擊后果嚴重性與撞擊點分布位置的關系。在150 m/s撞擊速度下，該飛機風擋各工況的計算損傷結果均表現(xiàn)出了對適航條款的符合性，然而在接觸力、風擋變形以及能量變化情況上顯示出了明顯的差異。本文通過參數(shù)化研究來了解這些差異性，并以此為參考挑選出具有代表性的撞擊點作為試驗靶點，可以為審定試驗提供參考，避免大量試驗所帶來的高成本。本文的參數(shù)化研究得到了以下結論：

(1) 越靠近飛機機身縱軸，鳥體與風擋之間的接觸力越大，飛機吸收的鳥體動能越大。

(2) 對于不同沖擊點而言，風擋的變形情況明顯受到邊界夾持部件的影響，不可忽視。

(3) 窗框結構不同部位的支持剛度存在差異，會影響風擋及窗框結構的吸能過程，且當著彈點靠近夾持邊界時，對應夾持部件出現(xiàn)應力集中現(xiàn)象。

(4) A2、A9、A11、B2四點適合作為適航驗證試驗靶點。