復合材料雷達罩鳥撞破壞流固耦合動響應分析

鄭涵天，王富生，岳珠峰

（西北工業(yè)大學 力學與土木建筑學院，西安 710129）

國內外的專家和學者對于鳥撞問題做過大量的探索和研究。在研究的方式上有實驗和計算機模擬；在鳥撞的計算方法上，主要有解耦和耦合算法［1］；在模擬鳥撞的工具方面，目前用的較多的有LS-DYNA，DYTRAN，ABAQUS，PAM-CRASH 等等。

Meguid等［2］建立了鳥撞發(fā)動機葉片有限元模型，分析了3種不同的鳥體本構模型對數(shù)值模擬結果的影響；Lavoie等［3］在進行鳥撞不同材料平板實驗的基礎上，把鳥體采用ALE，SPH兩種不同處理方式所得到的數(shù)值模擬結果和實驗數(shù)據(jù)從多個角度進行了對比分析；Airoldi等［4］針對采用Lagrange網(wǎng)格的鳥撞問題，對比理論計算探討了不同有限元網(wǎng)格劃分形式等不同參數(shù)對鳥撞問題精度的影響；Smojver等［5］進行了鳥撞飛機副翼的數(shù)值模擬，從不同速度，不同失效模式等多個角度對鳥撞復雜結構進行了對比分析；Georgiadis等［6］基于SPH方法對波音787的可動復合材料機翼后緣進行了鳥撞實驗和數(shù)值模擬，并將兩種結果進行了對比，分析了數(shù)值模擬的精度；Hanssen等［7］也在實驗基礎上采用流體—結構耦合的方法對鳥撞復合材料夾芯板進行了數(shù)值模擬，并分析了夾芯板的動力響應。謝宗蕻等［8］建立了鳥撞蜂窩夾芯結構雷達罩有限元模型，選取具有代表性的不同重量鳥體和鳥撞速度，得到了不同情況下的雷達罩動力響應；李玉龍等［9-10］則對不同飛機結構，采用不同算法對鳥撞問題進行了大量研究；王富生等［1，11-13］也對鳥撞問題的接觸算法，鳥體材料參數(shù)等方面進行了一定的分析和研究。

總之，國內外不同鳥撞問題的研究中，對鳥體網(wǎng)格描述上，ALE與SPH方法在模擬鳥撞這種高度非線性、瞬時大變形的過程上要比LAGRANGE方法更接近真實情況；由文獻［2］等可了解到，在一定速度范圍內，ALE方法盡管在計算上較難控制，但其在描述鳥體與空氣，鳥體與雷達罩的作用上其模擬結果更加可靠，能量的控制上更加準確。而且，基于流固耦合方法的復合材料雷達罩抗鳥撞問題幾乎沒有人進行過系統(tǒng)研究，本文即立意于此，計算工具采用LS-DYNA，基于流固耦合方法對采用殼單元的復合材料雷達罩鳥撞過程進行動響應數(shù)值模擬研究。

1 鳥撞復合材料雷達罩模型描述

1.1 鳥撞復合材料雷達罩幾何及有限元模型

采用ANSYS建立有限元模型（見圖1），并用LSPREPOST修改K文件。整個模型由4個part組成，鳥體周圍為空氣網(wǎng)格。由于LS-DYNA進行流固耦合分析一般分為兩種方法：一種是共節(jié)點，即Lagrange單元和ALE單元的邊界共節(jié)點，但實際上相當于施加了額外的約束，不適合大多數(shù)問題；另外一種即使用關鍵字*CONSTRAINED_LAGRANGE_IN_SOLID關鍵字把流體與固體單元耦合在一起，流體網(wǎng)格和結構網(wǎng)格和可以重疊，但在建立流體網(wǎng)格時需在周圍同時建立一個可以提供流體材料輸運、流體網(wǎng)格變形的ALE空間網(wǎng)格，如本文中鳥體周圍的空氣網(wǎng)格提供了鳥體的ALE流動空間。

模型中雷達罩的幾何模型是根據(jù)相關設計手冊自行設計的，在雷達罩的有限元模型中，為了避免產(chǎn)生過大網(wǎng)格畸變導致計算困難以及更加準確地描述鳥體和雷達罩的流固耦合作用過程，在雷達罩的撞擊點半徑0.25m范圍內網(wǎng)格加密，其它部分采取較稀疏網(wǎng)格；參照國內外研究人員的鳥撞模型，本文中鳥體的造型采用了中間圓柱兩邊半球的造型，鳥體的質量采取國際上通用的1.8kg，鳥體沖擊方式為正面沖擊雷達罩的中心點，模型的部分尺寸及鳥體與周圍空氣的ALE網(wǎng)格表示如圖2所示。

圖1 鳥撞雷達罩整體模型Fig.1 Global model of bird impact

圖2 鳥撞局部模型Fig.2 Local model of bird impact

1.2 主要參數(shù)說明

各parts所采用單元、材料及網(wǎng)格類型列于表1。

表1 各part屬性Tab.1 Properties of parts

1.2.1 鳥體/空氣參數(shù)

鳥撞雷達罩是一個瞬時高度非線性的碰撞過程，鳥體在碰撞瞬間呈現(xiàn)流體的性質，產(chǎn)生大變形，甚至破碎，因此采用了帶有狀態(tài)方程的空材料模型，以模擬碰撞過程中瞬間鳥體帶有流體特性的變形。參考相關資料［14］可知材料變形一般可分為2種類型：一種是變形過程中體積不變，一種是體積發(fā)生改變。應力張量可以由應力偏量和壓力兩部分組成：

其中：

其中：

對于任何材料都可以用應力偏量和壓力來描述它的應力張量。流體材料的變形過程需要同時使用兩種方式來描述材料，用本構模型來描述Δσ'ij和Δε'ij的關系，用狀態(tài)方程來描述ΔP和Δυ/υ的關系。本文用空材料*MAT_NULL描述具有流體行為的鳥體變形。該材料模式本身提供本構模型來描述材料的偏應力（粘性應力）：

然后使用狀態(tài)方程*EOS_LINEAR_POLYNOMIAL提供壓力行為應力組件，與上述材料偏應力本構模型共同提供材料的應力張量：

鳥體/空氣部分參數(shù)列于表2。

表2 鳥體/空氣參數(shù)Tab.2 Parameters of bird/air

在LS-DYNA中控制鳥體和空氣單元算法都采用11號算法，即具有多物質的單點ALE算法。將鳥體和空氣的使用關鍵字*ALE_MULTI_MATRIAL_GROUP定義在一個物質組中，鳥體周圍的空氣網(wǎng)格作為鳥體的ALE空間網(wǎng)格，使用關鍵字*ALE_REFERENCE_SYSTEM_GROUP定義外圍空氣ALE網(wǎng)格可以平動、旋轉和擴張來匹配鳥體的變形，使空氣ALE網(wǎng)格可以充分包換鳥體的流變或碎片。將鳥體和空氣定義為一個part組，并對其施加150 m/s的相對速度。

1.2.2 復合材料雷達罩參數(shù)

復合材料雷達罩采用蜂窩夾芯結構，夾芯結構的鋪層形式如圖3所示。為了簡化計算模型，設置雷達罩為等厚度形式，每個殼單元的厚度為0.007 6 m，剪切因子設置為1.0，上下表面的材料百分比厚度為0.105%，中間的蜂窩夾層百分比厚度為0.79%。在雷達罩邊緣節(jié)點全部施加了固支約束。殼單元中設置3個積分點，分別和上下蒙皮及夾芯對應，算法采用DYNA默認的Belytschko-Tsay算法。蒙皮和夾芯材料主要參數(shù)列于表3。

表3 蒙皮和夾芯材料參數(shù)Tab.3 Parameters of skin and core materials

參考 LS-DYNA 理論手冊［15］，22號材料*MAT_COMPOSITE_DAMAGE，其破壞模型為 Chang-Chang Composite Failure Model［16］，參照具體模型對 22 號材料的失效準則及破壞模式進行了非常詳盡的研究。在這種破壞模型中，共有3種破壞準則對應不同情況下的破壞模式；5個參數(shù)用于這3個破壞模式，分別是：S1為縱向拉伸強度、S2為橫向拉伸強度、S12為剪切強度、C2為橫向壓縮強度、α為非線性剪切應力參數(shù)。

圖3 復合材料鋪層Fig.3 Layers of composite

S1，S2，S12和C2由相應材料的強度試驗數(shù)據(jù)獲取，α由材料的剪應力應變試驗數(shù)據(jù)取得。應變由以下公式得到：

在此基礎上得到的3個失效準則：

基體開裂的失效準則：

壓縮失效準則：

纖維的斷裂失效準則：

1.2.3 流固耦合參數(shù)

流固耦合的計算過程主要涉及到兩個方面：流體的處理方式以及流體和結構間的耦合算法。本文中鳥體作為流體采用了ALE網(wǎng)格，ALE算法的控制方程可以由下列守恒方程給定：

（1）質量守恒方程：

（2）動量守恒方程：

控制固定域上的牛頓流體流動問題的增強形式由控制方程和對應的初始及邊界條件組成，控制流體問題的方程是Navier-Stokes方程的ALE描述：

應力張量σij可表述為下列方程：

（3）能量守恒方程：

本文中流體與結構間的耦合算法采用Euler-Lagrange耦合算法中的罰函數(shù)耦合方式，結構為從物質，流體為主物質，通過罰函數(shù)耦合，將結構的相關參量傳遞給流體單元。罰函數(shù)耦合系數(shù)追蹤拉格朗日節(jié)點和歐拉流體物質位置間的相對位移d，檢查每一個從節(jié)點對主物質表面的貫穿，如果從節(jié)點不出現(xiàn)貫穿，就不進行任何操作；如果發(fā)生從節(jié)點對主物質表面的貫穿，界面力F就會分布到歐拉流體的節(jié)點上。界面力的大小與發(fā)生貫穿數(shù)量成正比：

式中ki表示基于主從節(jié)點質量模型特性的剛度系數(shù)。由于在每一個時間積分步上都要對界面力F進行求解，所以可以認為F是一個外部體力。因此，每一時間步上都可以對總節(jié)點力Fn進行求解，并由Fn引起了結構的加速度、速率和位移的變化。

ALE控制方程配合時間積分、平滑算法、對流算法以及罰函數(shù)流固耦合算法等共同處理流體與結構的耦合過程，在LS-DYNA程序中，上述算法的運用均通過關鍵字定義的形式實現(xiàn)。

在LS-DYNA的流固耦合關鍵字定義中，除了用關鍵字*ALE_MULTI_MATRIAL_GROUP將鳥體和空氣定義在一個物質組中以及用關鍵字*SECTION_SOLID_ALE定義鳥體及空氣的ALE網(wǎng)格外，還需要用關鍵字*CONSTRND_LAGR

ANGE_IN_SOLID定義流體和結構間的耦合，即鳥體和雷達罩之間的耦合。耦合類型為4號，即殼單元和實體單元的罰函數(shù)耦合，主部為鳥體，從部為雷達罩，同時定義時間積分、對流算法等參數(shù)。

另外，計算時長取0.002 s，為防止計算過程中網(wǎng)格畸變過大導致的計算中止，時間步設置為2e-7 s，并加入沙漏控制。

2 計算結果及分析

鳥體以150 m/s的相對速度撞擊雷達罩，鳥體材料在空氣網(wǎng)格中流動，鳥體與空氣兩種不同物質在鳥撞過程中持續(xù)發(fā)生作用，圖4給出了計算開始時刻鳥體與空氣流體互相作用所產(chǎn)生的應力變化。

圖4 t=0 s時刻鳥體/空氣Von Mises應力圖Fig.4 The Von mises stress of bird/air at t=0 s

鳥體與雷達罩在t=0.000 45 s時發(fā)生耦合，鳥體的ALE網(wǎng)格節(jié)點在雷達罩的Lagrange網(wǎng)格表面發(fā)生微小貫穿，DYNA程序根據(jù)關鍵字定義的罰函數(shù)耦合算法判斷貫穿量并跟蹤網(wǎng)格節(jié)點，檢查貫穿范圍，并按照關鍵字定義的罰函數(shù)耦合系數(shù)在流體和結構網(wǎng)格之間執(zhí)行傳遞參量、求解及分配界面力等操作。圖5給出了鳥體ALE網(wǎng)格與雷達罩Lagrange網(wǎng)格耦合開始時刻鳥體與雷達罩的Von Mises應力圖。

圖5 耦合開始時刻（t=0.000 45 s）鳥體與雷達罩Von mises應力圖Fig.5 The Von Mises stress at t=0.000 45 s

圖6給出了不同時刻雷達罩破壞響應和鳥體變形圖，圖7給出了相應時刻雷達罩和鳥體的Von Mises應力云圖。從圖6中可以看出，鳥體網(wǎng)格與雷達罩網(wǎng)格發(fā)生耦合后，雷達罩中心區(qū)域在鳥體沖擊的作用下不斷發(fā)生變形，并向周圍擴散，在沖擊過程末段，雷達罩結構出現(xiàn)失效單元，雷達罩發(fā)生破損；同時，鳥體也發(fā)生流變，縱向長度不斷變小，鳥體材料向周圍擴散，產(chǎn)生破碎趨勢，并在雷達罩破損后從破損處滲出。從圖7中可以看出，雷達罩結構在沖擊作用下產(chǎn)生應力波，隨著耦合過程的進行，應力波從中心區(qū)域不斷向周圍擴散，在耦合作用前期，較大應力區(qū)域出現(xiàn)在雷達罩被撞擊的中心區(qū)域，但隨著雷達罩的破壞和鳥體的變形，產(chǎn)生較大應力的區(qū)域向撞擊點的周圍轉移；同時，鳥體由于受到雷達罩結構的影響，在與雷達罩結構發(fā)生耦合的區(qū)域也產(chǎn)生應力波，并隨著鳥體的變形，不斷向鳥體材料的周圍擴散。

圖6 不同時刻雷達罩的破壞響應和鳥體變形Fig.6 The damage response of radome and deformation of the bird at different times

圖7 不同時刻的Von mises應力云圖Fig.7 Von mises stress at different times

圖8給出幾個典型單元的具體位置，XY平面應變時程曲線如圖9所示。圖6中選取了雷達罩上的6個殼單元 37 091（A）、37 061（B）、37 031（C）、37 200（D）、37 193（E）、38 206（F），圖9分別給出所選取單元的上表面蒙皮（upper surface），夾芯（mid surface）以及下表面蒙皮（lower surface）的XY平面應變時程曲線圖。對比圖中曲線可知：

（1）由ABC和DEF兩組曲線的對比可知，ABC三條曲線所對應的單元在整個撞擊過程中產(chǎn)生的應變遠大于DEF三條曲線對應的單元。

（2）由AB兩條曲線和C曲線的對比可知，在雷達罩正面承受撞擊的中心區(qū)域，雷達罩的材料在沖擊過程中需要承受拉伸和壓縮的不斷變化，部分單元則以壓縮為主，尤其上表面蒙皮材料的變化幅度最為劇烈；與中心撞擊區(qū)域相鄰較近的單元無論上表面蒙皮、夾芯和下表面蒙皮均以壓縮為主要變形形式。

（3）由A曲線和B曲線的對比可知，雖然同在中心撞擊區(qū)域，但不同位置的雷達罩材料在同一時間段內產(chǎn)生的變形也不相同；而且對于不同單元，其上中下3層材料的變形趨勢也不相同，如A曲線所對應單元在最后1/4時間段內上表面蒙皮材料產(chǎn)生拉伸變形，夾芯和下表面蒙皮則產(chǎn)生大幅度的壓縮變形。

圖8 被選取的典型單元位置Fig.8 The location of selected elements

圖9 不同單元的XY平面應變時程曲線Fig.9 The time history of XY-strain for different elements

3 結論

本文采用流固耦合方法對鳥撞復合材料雷達罩的沖擊過程進行了數(shù)值模擬，并給出了數(shù)值模擬的結果，通過分析主要得到了以下結論：

（1）在鳥體周圍附加空氣ALE網(wǎng)格，并結合流固耦合算法所得到的數(shù)值模擬結果使鳥體的變形更加貼近真實情況，從而也使雷達罩的動力響應情況更加真實可靠。

（2）鳥體撞擊雷達罩導致材料直接破壞的情況主要集中在中心撞擊區(qū)域，而在應力的變化過程當中，并非只有中心撞擊區(qū)域會出現(xiàn)應力峰值，應力的峰值區(qū)域也會出現(xiàn)在中心撞擊區(qū)域的周圍。

（3）中心撞擊區(qū)域的材料在整個撞擊的過程當中都要不斷承受拉伸和壓縮的變化，同時部分材料不同層間要承受不同的應變趨勢。針對這種破壞響應，在復合材料的抗鳥撞設計當中應針對不同材料層間所承受的不同應變形式來設計復合材料夾芯雷達罩的夾層形式。

［1］王富生，岳珠峰，王新軍，等.鳥撞飛機風擋的一體化數(shù)值模擬技術［J］.振動與沖擊，2007，26（5），107-111.

［2］ Meguid S A，Mao R H，Ng T Y.FE analysis of geometry effects of an artificial bird striking an aeroengine fan blade［J］.International Journal of Impact Engineering，2008，35：487-498.

［3］ Lavoie M A，Gakwaya A，Ensan M N，et al.Bird＇s substitute tests results and evaluation of available numerical methods［J］.International Journal of Impact Engineering，2009，36：1276-1287.

［4］Airoldi A，Cacchione B.Modelling of impact forces and pressures in Lagrangian bird strike analyses［J］.International Journal of Impact Engineering，2006，32：1651-1677.

［5］ Smojver I，Ivancěvic D.Numerical simulation of bird strike damage prediction in airplane flap structure［J］.Composite Structures，2010，92：2016-2026.

［6］ Georgiadis S，Gunnion A J，Thomson R S，et al.Brid-strike simulation for certification of the boeing 787 composite moveable trailing edge［J］.Composite Structures，2008，86：258-268.

［7］ Hanssen A G，Girard Y，Olovsson L，et al.A numerical model for bird strike of aluminium foam-based sandwich panels［J］.International Journal of Impact Engineering，2006，32：1127-1144.

［8］謝宗蕻，卞文杰，蜂窩夾芯結構雷達罩鳥撞有限元分析與模擬［J］.爆炸與沖擊，1999，19（3）：235-242.

［9］毋 玲，郭英男，李玉龍.蜂窩夾芯雷達罩結構的鳥撞數(shù)值分析［J］.爆炸與沖擊，2009，29（6）：642-647.

［10］ Li Y L，Zhang Y K，Xue P，et al.Study of similarity law for bird impact on structure［J］.Chinese Journal of Aeronautics，2008，21：512-517.

［11］ Wang F S，Yue Z F.Numerical simulation of damage and failure in aircraft windshield structure against bird strike［J］.Materials and Design，2010，31：687-695.

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

［13］王富生，馮震宙，岳珠峰.飛機風擋鳥撞破壞的一種耦合接觸碰撞數(shù)值模擬［J］.振動與沖擊，2008，27（10）：166-169.

［14］ 李裕春，時黨勇，趙 遠.ANSYS11.0/LS-DYNA基礎理論與工程實踐［M］.北京：中國水利水電出版社，2008.

［15］ Livermore Software Technology Corporation.LS-DYNA Theoretical Manual［M］.1998.

［16］ Iannucci L.Progressive failure modeling of woven carbon composite under impact［J］.International Journal of Impact Engineering，2006，32（6）：1013-1043.