蜂窩夾芯雷達罩抗鳥撞沖擊性能仿真分析

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

（1.鄭州航空工業管理學院 航空宇航學院，河南 鄭州 450046；2.中國兵器工業標準化研究所，北京 100089；3.中國兵器科學研究院，北京 100089；4.鄭州航空工業管理學院 河南省通用航空技術重點實驗室，河南 鄭州 450046）

0 引 言

現代飛機的機頭外部均裝備雷達天線罩以保護雷達天線安全。然而，雷達天線罩作為飛機主要的迎風面結構之一，在服役過程中不可避免地會受到鳥體、冰雹等的高速撞擊，因此其抗沖擊性能對飛機安全有著重要影響。

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

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

1 雷達天線罩鳥撞有限元模型建立

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

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

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

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

圖1 鳥體形狀尺寸及其計算模型

表1 四類鳥體質量和密度參數

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

拉格朗日有限元網格單元采用式（2）給出的彈塑性硬化模型描述鳥體的動態力學行為。

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

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

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

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

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

空材料模型一般配以多項式狀態方程進行描述。鳥體體積變形與壓力之間關系由式(9)給出，材料的壓縮比由式(10)描述，對應的鳥體狀態方程參數如表3所示。

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

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

1.2 無人機雷達天線罩有限元模型建立

本文的無人機雷達天線罩采用三明治夾芯結構，其結構示意圖見圖2。天線罩內外兩層為0.8mm的2024T3Al 合金板，中間夾芯層為8.4mm 的Nomex蜂窩鋁材料[19]。

圖2 無人機雷達天線罩的三明治夾芯結構示意圖

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

圖3 無人機雷達天線罩

雷達罩內、外層材料的動態力學行為由式(11)所給出的基于連續損傷理論的104 號材料模型(*MAT_DAMAGE_1)進行描述。

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

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

天線罩夾芯層泡沫鋁結構材料采用154 號材料（* MAT_DESH PANDE_FLECK_FOAM）進行描述。其中，材料屈服強度σy和材料參數的關系由式(12)給出，壓實應變由式(13)定義給出。

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

其中，天線罩夾芯層泡沫鋁結構的材料參數取值如下：ρ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]

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

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

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

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

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

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

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

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

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

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

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

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

2.4 天線罩正面中間鳥撞位置處三明治結構的耗能情況

無人機雷達天線罩內外層以及夾芯層吸收的內能隨時間變化對比曲線如圖8所示。

圖8 工況1條件下天線罩三明治結構內能—時間曲線對比

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

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

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

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

3 結 論

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

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

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