李國舉，張潛銳，侯暉東，張昕喆,4

（1.鄭州航空工業(yè)管理學(xué)院 航空宇航學(xué)院，河南 鄭州 450046；2.中國兵器工業(yè)標(biāo)準(zhǔn)化研究所，北京 100089；3.中國兵器科學(xué)研究院，北京 100089；4.鄭州航空工業(yè)管理學(xué)院 河南省通用航空技術(shù)重點(diǎn)實(shí)驗(yàn)室，河南 鄭州 450046）

0 引 言

現(xiàn)代飛機(jī)的機(jī)頭外部均裝備雷達(dá)天線罩以保護(hù)雷達(dá)天線安全。然而，雷達(dá)天線罩作為飛機(jī)主要的迎風(fēng)面結(jié)構(gòu)之一，在服役過程中不可避免地會受到鳥體、冰雹等的高速撞擊，因此其抗沖擊性能對飛機(jī)安全有著重要影響。

在所有高速撞擊載荷中，鳥撞載荷是最為危險(xiǎn)的極限載荷，因此，開展鳥體撞擊飛機(jī)的研究與測試尤為重要。當(dāng)前，對于飛機(jī)鳥撞研究主要分為鳥撞試驗(yàn)和鳥撞仿真模擬兩類方法[1]?？琢钣碌萚2]采用鳥撞試驗(yàn)裝置對碳纖維復(fù)合材料層合板進(jìn)行撞擊測試，重點(diǎn)對材料的鋪層方式與其耗能之間的關(guān)系進(jìn)行了系統(tǒng)研究。陳琨等[3]采用空氣炮發(fā)射明膠鳥彈撞擊復(fù)合材料蜂窩夾芯板，借助高速攝影研究了夾芯結(jié)構(gòu)在軟體高速沖擊下的損傷情況。但由于高速攝影、空氣炮等試驗(yàn)設(shè)備對場所有一定要求且價格高昂，導(dǎo)致研究成本過高，且操作過程復(fù)雜、耗時。隨著計(jì)算機(jī)性能和計(jì)算力學(xué)的不斷發(fā)展，有限元仿真方法逐漸成為分析鳥撞問題的主要手段[4-6]。鄭涵天等[7]采用有限單元法（finite element method，簡稱FEM），研究了雷達(dá)罩抗鳥撞動力響應(yīng)及結(jié)構(gòu)受損情況。劉玄等[8]采用FEM 方法對無人機(jī)雷達(dá)天線罩在抗鳥撞沖擊過程的耗能進(jìn)行了研究。然而，F(xiàn)EM 方法在結(jié)構(gòu)受到大變形時容易引起單元畸變，難以有效模擬高速鳥撞過程中的鳥體塑性流動和鋪展過程，而光滑粒子流體動力學(xué)(smooth particle hydrodynamics，簡稱SPH)方法可以有效克服FEM 方法因結(jié)構(gòu)大變形引起的計(jì)算不穩(wěn)定問題，能夠較好地模擬鳥體在高速沖擊過程中的飛濺效果，因此已被廣泛應(yīng)用于鳥撞仿真分析中[9-12]。毋玲等[13]采用SPH粒子算法，對雷達(dá)罩抗高速鳥撞過程中的耗能隨撞擊位置的變化情況進(jìn)行了系統(tǒng)研究。于永強(qiáng)等[14]同樣采用SPH 粒子研究了復(fù)合材料層合板受鳥撞擊的損失失效情況。研究發(fā)現(xiàn)，鳥體與天線罩之間的接觸傳力對鳥撞過程分析非常關(guān)鍵，而當(dāng)前單一SPH方法對于接觸模擬尚不夠精確，限制了其在模擬鳥撞過程中的進(jìn)一步應(yīng)用。

為彌補(bǔ)上述研究中單一FEM 與SPH 方法的不足并充分發(fā)揮兩種方法各自的優(yōu)勢，本文引入FEMSPH 耦合算法[15-18]，開展對無人機(jī)天線雷達(dá)罩不同位置的抗鳥撞沖擊動態(tài)響應(yīng)過程仿真模擬研究。該算法在鳥體材料失效前以拉格朗日有限元網(wǎng)格單元模擬鳥體與無人機(jī)雷達(dá)天線罩之間的接觸傳力過程；鳥體失效后有限元網(wǎng)格轉(zhuǎn)換為SPH 粒子，以模擬其后續(xù)的飛濺和鋪展等大變形過程。由此得到的仿真結(jié)果相較于單一FEM 與SPH 方法更加真實(shí)可信，可為鳥撞仿真模擬以及飛行器結(jié)構(gòu)設(shè)計(jì)領(lǐng)域提供理論依據(jù)與參考數(shù)據(jù)。

1 雷達(dá)天線罩鳥撞有限元模型建立

1.1 簡化鳥體模型選擇及有限元模型建立

能夠描述真實(shí)鳥體結(jié)構(gòu)特征是開展鳥撞仿真分析的首要前提。在實(shí)際有限元仿真計(jì)算中，通常將真實(shí)鳥體形狀簡化為圓柱與其端部的兩個半球形成的半球端圓柱體結(jié)構(gòu)。其中，圓柱體的高度是半球體半徑的兩倍，半球端圓柱體的特征幾何尺寸和鳥體的質(zhì)量、密度的關(guān)系由式(1)給出。

式(1)中，L為圓柱體的高度，R為半球體半徑，m為鳥體的質(zhì)量，ρ為鳥體的密度。

國際上對鳥撞研究有著較多的規(guī)范，一般選取四類鳥體進(jìn)行表征，這四類鳥體的質(zhì)量和密度參數(shù)如表1 所示。因此，本文同樣選擇家雞為研究對象，將其質(zhì)量m=1.8Kg，密度ρ=977.71Kg/m3帶入式(1)，計(jì)算可得家雞簡化模型中圓柱高度L=112mm，半球半徑R=56mm，如圖1(a)所示。最終通過CATIA 軟件建模獲得鳥體幾何模型如圖1( b)所示。采用LSPrepost 軟件建立的鳥體有限元網(wǎng)格模型如圖1(c)所示，并進(jìn)一步基于有限單元中心點(diǎn)獲得對應(yīng)的SPH粒子模型如圖1(d)所示。鳥體的FEM 單元為六面體實(shí)體單元，其網(wǎng)格尺寸控制在3mm—7mm。最終獲得鳥體有限元網(wǎng)格單元數(shù)量為56 832個，節(jié)點(diǎn)數(shù)量為58 407個。

圖1 鳥體形狀尺寸及其計(jì)算模型

表1 四類鳥體質(zhì)量和密度參數(shù)

本文將分別采用FEM、SPH 和FEM-SPH 耦合方法描述鳥體材料在沖擊過程中的變形失效情況，并對比分析拉格朗日有限元網(wǎng)格模擬碰撞接觸的精確性以及SPH 粒子允許大變形特性在鳥撞問題中的影響。

拉格朗日有限元網(wǎng)格單元采用式（2）給出的彈塑性硬化模型描述鳥體的動態(tài)力學(xué)行為。

式（2）中：σy為塑性應(yīng)力；σs為屈服強(qiáng)度應(yīng)力；EP和εP分別為切線模量和塑性應(yīng)變。鳥體的失效準(zhǔn)則采用最大塑性應(yīng)變進(jìn)行描述。

其中，鳥體的彈塑性硬化模型參數(shù)取值如下：E=1×104MPa，v=0.3，σs=1.0MPa，Ep=5.0MPa，εcr=1.22[19]。

考慮到拉格朗日有限元網(wǎng)格單元變形過大時會因產(chǎn)生負(fù)體積而導(dǎo)致計(jì)算中止，因此本文設(shè)置當(dāng)塑性應(yīng)變超過1.22 時單元被刪除，從而保證鳥撞仿真過程中計(jì)算的穩(wěn)定性。對于采用SPH 粒子算法的鳥體模型，則直接采用空材料（*MAT_NULL）描述其力學(xué)行為。而對于采用FEM-SPH 耦合算法的鳥體模型，單元未失效前采用拉格朗日有限元網(wǎng)格，此時用彈塑性硬化模型描述鳥體動態(tài)力學(xué)行為，而當(dāng)鳥體有限元單元塑性應(yīng)變達(dá)到1.22 時自動轉(zhuǎn)換為SPH 粒子，并將失效單元節(jié)點(diǎn)位移、質(zhì)量、速度、應(yīng)力和聲速等物理量賦給轉(zhuǎn)換后的SPH 粒子，相應(yīng)轉(zhuǎn)換關(guān)系見式(3)—(7)。

式(3)—(7)中，角標(biāo)p 和n 分別代表SPH 粒子和有FEM 的節(jié)點(diǎn)，而角標(biāo)e 和g 則分別單元和高斯點(diǎn)；x和m為鳥體模型的位移和質(zhì)量；Ne為與節(jié)點(diǎn)相關(guān)的FEM 單元總數(shù)，Nn為任意單元的節(jié)點(diǎn)數(shù)量；ρei、Vei、σei分別為單元材料的密度、體積和應(yīng)力張量；wgj和σgj為單元序號i內(nèi)部的高斯積分點(diǎn)g 處對應(yīng)的加權(quán)系數(shù)和應(yīng)力張量；Ng為任一FEM 單元所包含的高斯積分點(diǎn)數(shù)量。SPH 粒子采用空材料（*MAT_NULL）描述其隨后的飛濺和鋪展流動效果，鳥體失效后的偏應(yīng)力由式(8)給出。

式(8)中，σ為鳥體失效后的偏應(yīng)力;υd為是動態(tài)粘度;D為變形率張量;P為壓力;II為張量不變量。

空材料模型一般配以多項(xiàng)式狀態(tài)方程進(jìn)行描述。鳥體體積變形與壓力之間關(guān)系由式(9)給出，材料的壓縮比由式(10)描述，對應(yīng)的鳥體狀態(tài)方程參數(shù)如表3所示。

式(9)—(10)中，U為鳥體單位體積的內(nèi)能;μ為材料的壓縮比;ρ和ρ0分別為材料的當(dāng)前密度和初始密度;C0—C6均為材料參數(shù)。

其中，鳥體的狀態(tài)方程參數(shù)取值如下：C0=0，C1=2.25×103MPa，C2—C6均為0。

1.2 無人機(jī)雷達(dá)天線罩有限元模型建立

本文的無人機(jī)雷達(dá)天線罩采用三明治夾芯結(jié)構(gòu)，其結(jié)構(gòu)示意圖見圖2。天線罩內(nèi)外兩層為0.8mm的2024T3Al 合金板，中間夾芯層為8.4mm 的Nomex蜂窩鋁材料[19]。

圖2 無人機(jī)雷達(dá)天線罩的三明治夾芯結(jié)構(gòu)示意圖

通過CATIA 軟件建立三維實(shí)體模型并采用Hypermesh 對其進(jìn)行網(wǎng)格剖分，獲得無人機(jī)雷達(dá)天線罩有限元網(wǎng)格模型如圖3 所示。為了體現(xiàn)天線罩的三維結(jié)構(gòu)情況，圖中分別給出了如圖3(a)所示的正視圖和如圖3(b)所示的側(cè)視圖。其中夾芯層采用六面體實(shí)體網(wǎng)格單元(Solid 164)進(jìn)行剖分，截面方向的網(wǎng)格尺寸設(shè)置為5mm,厚度方向尺寸為2.8mm，以確保厚度方向至少劃分3 層實(shí)體單元。最終獲得夾芯層結(jié)構(gòu)的六面體單元數(shù)量為115 680 個，節(jié)點(diǎn)數(shù)量為155 848 個。對內(nèi)、外兩層鋁合金面板采用殼單元(shell 163)進(jìn)行面網(wǎng)格劃分，單元尺寸同樣設(shè)為5 mm，最終獲得內(nèi)、外層單元數(shù)量為38 880，節(jié)點(diǎn)數(shù)為39 284,并將內(nèi)、外層殼單元的厚度均設(shè)置為0.8 mm。

圖3 無人機(jī)雷達(dá)天線罩

雷達(dá)罩內(nèi)、外層材料的動態(tài)力學(xué)行為由式(11)所給出的基于連續(xù)損傷理論的104 號材料模型(*MAT_DAMAGE_1)進(jìn)行描述。

式(11)中，r為損傷累積的塑性應(yīng)變判據(jù)為等效塑性應(yīng)變率;σ0為材料初始屈服強(qiáng)度;Q1、Q2、C1、C2為材料各向同性硬化系數(shù)。

其中，天線罩外殼的材料模型參數(shù)取值如下：σ0=364.5MPa，Q1=334.7MPa，C1=6.16，Q2=0，C2=0，r=0.18。

天線罩夾芯層泡沫鋁結(jié)構(gòu)材料采用154 號材料（* MAT_DESH PANDE_FLECK_FOAM）進(jìn)行描述。其中，材料屈服強(qiáng)度σy和材料參數(shù)的關(guān)系由式(12)給出，壓實(shí)應(yīng)變由式(13)定義給出。

式(12)—(13)中，α2、γ和β都為材料參數(shù);σp為平臺應(yīng)力;εD為壓實(shí)應(yīng)變;ε?為等效應(yīng)變;ρf為泡沫材料當(dāng)前密度;ρf0為泡沫材料初始密度。

其中，天線罩夾芯層泡沫鋁結(jié)構(gòu)的材料參數(shù)取值如下：ρf0=300Kg/m3，E=1.5×103MPa，α=2.1，γ=6.10MPa，εD=2.2，α2=38.1MPa，β=3.1，σp=4.41MPa，εcr=0.2。[20]

本文中，鳥體相對天線罩結(jié)構(gòu)的速度設(shè)為150m/s，鳥體與天線罩正面相撞，撞擊點(diǎn)沿高度方向（Z 軸正方向）排列，設(shè)定為5 組，每組撞擊位置在Z 方向的位移間隔為100mm，如圖4所示。天線罩模型邊界采用固支約束，鳥體拉格朗日有限元網(wǎng)格與天線罩結(jié)構(gòu)使用自動面面接觸算法進(jìn)行接觸控制，鳥體網(wǎng)格失效后轉(zhuǎn)換為SPH 粒子與天線罩結(jié)構(gòu)采用自動點(diǎn)面接觸算法進(jìn)行接觸控制。模型整體求解時間設(shè)置為6ms，每0.1ms 輸出一個結(jié)果，時間步長系數(shù)設(shè)置為0.6。

圖4 無人機(jī)雷達(dá)天線罩的鳥撞不同位置有限元模型示意圖

2 天線罩鳥撞過程有限元仿真分析

2.1 天線罩正面中間鳥撞位置處鳥體不同算法仿真結(jié)果對比

天線罩正面中間位置處不同鳥體描述算法對應(yīng)的天線罩撞擊點(diǎn)整體位移—時間曲線及其在T=3ms時刻的鳥體變形鋪展?fàn)顟B(tài)如圖5所示。由圖中可知，在時間為0ms—1.5ms 時，三種算法描述的鳥體撞擊后獲得的天線罩撞擊點(diǎn)整體位移—時間曲線基本一致。在同一時刻對應(yīng)的撞擊點(diǎn)整體位移中，F(xiàn)EM 算法獲得的位移最大，表明鳥體初始接觸天線罩體時FEM 算法接觸最為充分。隨著時間的增加，同一時刻SPH算法獲得的位移最大，F(xiàn)EM最小。觀察T=3ms時鳥體變形鋪展過程可知，隨著FEM 算法中單元塑性變形增大并達(dá)到失效應(yīng)變而被不斷刪除，最終由于質(zhì)量損失過大而無法對天線罩繼續(xù)作用。而SPH和FEM-SPH 算法中，不存在單元刪除，因此兩者獲得位移均高于FEM 算法。對比SPH 算法和FEMSPH 算法在T=3ms時的變形鋪展過程可知，SPH 算法鳥體已經(jīng)完全鋪展開并向四周飛散。而在FEM-SPH耦合算法中，由于未失效拉格朗日有限單元之間的相互連接，使得鳥體模型完全鋪展后并未如SPH 算法一樣快速飛散，這表明FEM-SPH 耦合算法更能精細(xì)描述鳥體撞擊到天線罩上的變形、失效以及鋪展過程。

圖5 天線罩正面中間位置處不同鳥體算法對應(yīng)的撞擊點(diǎn)整體位移—時間曲線以及T=3 ms對應(yīng)的鳥體變形鋪展過程

2.2 天線罩正面中間鳥撞位置處鳥體鋪展和天線罩變形過程

天線罩正面中間位置處隨鳥撞不同時刻對應(yīng)的鳥體鋪展和天線罩整體位移云圖如圖6 所示。從圖中可知，在T=2ms 時，鳥體與天線罩開始接觸并產(chǎn)生局部變形，鳥體頭部有限元網(wǎng)格單元已經(jīng)轉(zhuǎn)化為SPH粒子并開始在天線罩四周鋪展，此時，天線罩結(jié)構(gòu)整體的位移峰值為83mm，見圖6(a)。在T=4ms時，天線罩位移已經(jīng)達(dá)到最大值96mm，且鳥體大部分FEM 單元已經(jīng)失效轉(zhuǎn)換為SPH 粒子并向天線罩上部鋪展，此時的鳥體鋪展與飛濺效果已經(jīng)比較明顯，如圖6(b)所示。在T=6ms時，鳥體網(wǎng)格單元基本上全部轉(zhuǎn)換為SPH 粒子，鳥體飛濺效果更加顯著，并開始整體滑離出天線罩表面，如圖6(c)所示。

圖6 天線罩正面中間位置鳥撞不同時刻鳥體鋪展和天線罩位移變形云圖

2.3 天線罩正面中間鳥撞位置處外層鋁合金板的應(yīng)力演變過程

天線罩正面中間位置處鳥撞不同時刻對應(yīng)的天線罩外層鋁合金板Von-Mises等效應(yīng)力云圖如圖7所示。從圖中可知，在T=0.1ms 時，鳥體與天線罩開始接觸，此時撞擊點(diǎn)位置等效應(yīng)力最大值達(dá)到了383.0MPa，如圖7(a)所示。隨著時間的增加，當(dāng)T=0.8ms時，鳥撞接觸位置處產(chǎn)生應(yīng)力集中，最大等效應(yīng)力增加到427.3MPa，此時由于鳥體頭部單元尚未失效，整個接觸區(qū)域呈現(xiàn)出應(yīng)力值集中現(xiàn)象，如圖7(b)所示。當(dāng)T=1.3ms 時，應(yīng)力峰值增加到431.8MPa，此時與鳥體頭部上側(cè)位置轉(zhuǎn)換為SPH 粒子，導(dǎo)致此處應(yīng)力值降低，僅為316MPa，遠(yuǎn)小于其他接觸區(qū)域應(yīng)力值，如圖7(c)所示。在T=3.3ms時，鳥體頭部接觸區(qū)域應(yīng)力峰值達(dá)到整個鳥撞過程最大值452.0MPa，此時鳥體頭部一部分單元已經(jīng)轉(zhuǎn)換為SPH 粒子并呈現(xiàn)出流體的鋪展和飛濺等流動性特征，使得這些接觸區(qū)域應(yīng)力值穩(wěn)定在300MPa 左右，如圖7(d)所示。在T=5.4ms 時，鳥體幾乎完全轉(zhuǎn)化為SPH 粒子，即主要呈現(xiàn)其流動性特征，使得應(yīng)力集中轉(zhuǎn)移到鳥撞形成的坑體周圍區(qū)域，最大值也降至408MPa，如圖7(e)所示。在T=6.0ms時，由于大部分SPH 粒子已經(jīng)從鳥撞過程中形成的坑體內(nèi)滑移出去，因此整體接觸區(qū)域應(yīng)力值較小，僅為200MPa左右，如圖7(f)所示。

圖7 天線罩正面中間位置鳥撞不同時刻外層鋁合金板Von-Mises等效應(yīng)力云圖

2.4 天線罩正面中間鳥撞位置處三明治結(jié)構(gòu)的耗能情況

無人機(jī)雷達(dá)天線罩內(nèi)外層以及夾芯層吸收的內(nèi)能隨時間變化對比曲線如圖8所示。

圖8 工況1條件下天線罩三明治結(jié)構(gòu)內(nèi)能—時間曲線對比

從曲線對比可知，在T=4.0ms 時，天線罩內(nèi)外層及夾芯層內(nèi)能均達(dá)到最大值，后續(xù)隨時間變化保持不變。此時，泡沫鋁夾芯層內(nèi)能達(dá)到3.6×103J，而內(nèi)、外鋁合金面板內(nèi)能僅為775J 和682J，這表明鳥撞過程中沖擊主要由泡沫鋁夾芯層吸收，同時也表明泡沫鋁夾芯層材料在天線罩結(jié)構(gòu)中起著吸收鳥撞沖擊能力的重要作用。因此為了有效保護(hù)天線罩內(nèi)部的雷達(dá)和精密設(shè)備，需要對夾芯層材料進(jìn)行選擇和設(shè)計(jì)。

2.5 不同鳥撞位置下天線罩位移形變情況

不同鳥撞位置下無人機(jī)雷達(dá)天線罩對應(yīng)的最大整體位移云圖如圖9 所示。從圖中數(shù)據(jù)可知，天線罩鳥撞位置從上部向下移動過程中，天線罩結(jié)構(gòu)最大整體位移呈現(xiàn)先增加后減小的趨勢。無人機(jī)天線罩最危險(xiǎn)鳥撞位置為中間靠下的位置，此時對應(yīng)的最大整體位移峰值為123.5mm，較之鳥撞最上部位置的67.0mm，位移增加了約84%。鳥撞下部位置最大位移為103.4mm，稍高于中間位置的99.7mm。這表明鳥撞位置變化對天線罩整體變形影響很大，因此需要重點(diǎn)加強(qiáng)危險(xiǎn)位置的結(jié)構(gòu)設(shè)置。

圖9 不同鳥撞位置下對應(yīng)的天線罩模型最大整體位移云圖

3 結(jié) 論

（1）FEM-SPH 耦合方法能夠充分發(fā)揮FEM 法和SPH法各自的優(yōu)勢，能更精細(xì)描述鳥體撞擊到天線罩上的變形、失效以及鋪展飛濺過程。在鳥體撞擊天線罩正面中間位置的過程中，鳥體有限元網(wǎng)格單元不斷失效并轉(zhuǎn)變?yōu)镾PH 粒子，在T=4ms 時出現(xiàn)明顯的飛濺和鋪展現(xiàn)象。進(jìn)一步研究其外層鋁合金板應(yīng)力演化發(fā)現(xiàn)，在T=5.4ms 后由于SPH 飛濺出撞擊凹陷，應(yīng)力主要集中在凹陷周圍，凹陷底部應(yīng)力顯著減小。

（2）在鳥體撞擊天線罩正面中間位置的過程中，泡沫鋁夾芯層材料耗能達(dá)到3666J，遠(yuǎn)大于內(nèi)外層鋁合金面板耗能，表明泡沫鋁夾芯層材料在天線罩結(jié)構(gòu)中起著吸收鳥撞沖擊力的重要作用。

（3）鳥體以相同速度撞擊天線罩不同位置處的有限元仿真結(jié)果表明，隨著天線罩鳥撞位置從頂部向下移動，天線罩結(jié)構(gòu)最大整體位移呈現(xiàn)先增加后減小的趨勢。無人機(jī)天線罩最危險(xiǎn)的鳥撞位置為中間靠下位置，此時最大位移為123.5 mm，為上部位置最大位移值的1.8倍。