多歐拉域耦合法在平尾鳥撞中的應用

胡文剛，林長亮，王剛，門坤發

航空工業哈爾濱飛機工業集團有限責任公司 飛機設計研究所, 哈爾濱 150066

鳥撞是指飛行器與天空中飛行鳥類相撞造成的飛行事故。目前，對于飛機機頭罩、風擋、機翼、尾翼前緣、發動機葉片等易遭受鳥撞的結構，國內外都開展了大量的試驗和分析研究。

Reza和Saeed[1-2]討論了不同鳥體形狀對沖擊的影響，通過數值計算與試驗對比，確定了仿真的鳥體形狀,并采用野鴨的物理模型，研究了真空、空氣和不同撞擊方向對鳥撞響應的影響。McCallum等[3-4]用天鵝物理尺寸建立鳥體，通過研究發現，含頭、長頸的鳥體模型能正確預測沖擊載荷及靶體破壞，且計算結果沖擊時程更長，沖擊力極值更大。Lakshmi[5]建立了真實鳥體模型，對不同部位取不同密度，通過研究發現，分析小鳥體鳥撞問題時可用傳統的簡化鳥體軀干模型，但在分析大鳥體鳥撞時鳥體中的其他部位不可忽視。Hedayati等[6-7]建立幾何形狀更準確精細的鳥體模型，研究真實形狀鳥體從各部位對結構沖擊的靶體損傷，并將真實形狀鳥體計算結果與替代模型計算結果及試驗結果進行了對比，結果表明真實鳥體模型計算結果與理論計算值差別較大，與試驗值吻合較好。Guida等[8]對設計的3種纖維金屬層合板前緣結構進行沖擊能量吸收性能測試，對比分析了3種結構的抗鳥擊性能;Dar等[9]研究了不同撞速和不同撞擊角度下風擋的鳥撞動響應。Liu等[10]采用有限元與光滑粒子流體動力學(Smoothed Particle Hydrodynamics, SPH)結合的方法對某商用飛機平尾前緣進行了鳥擊仿真，與試驗結果進行了驗證，并提出了在前緣結構中引入三角鋼筋構件，能夠顯著提高前緣結構的抗鳥撞性能。Caprioa等[11]研究了垂尾前緣的抗鳥撞性能，討論了不同材料系統如何滿足前緣結構的減重需要和抗鳥擊性能，采用經過試驗驗證的數值程序，分析了幾種典型前緣結構的抗鳥撞性能。Yu等[12]提出了一種確定鳥類撞擊造成最嚴重損傷的臨界位置的分析方法，利用PAM-CRASH軟件，分析了鳥體垂直于機翼前緣的沖擊損傷，指出傳統的鳥體平行于前緣沖擊得到的臨界速度不是結構抗鳥撞的最小速度。Xie等[13]采用霍普金森拉桿(SHTB)系統進行了平尾前緣常用的2024-T3和7075-T6鋁材不同應變速率下的拉伸試驗，建立了反映鋁合金應變率硬化效應的約翰遜-庫克模型。中國ARJ21垂尾抗鳥撞設計是通過與美國ASC公司合作通過適航審查的。Wilbeck和Barber[14]采用10%孔隙率的明膠，制作2∶1的圓柱體來模擬鳥體，通過大量試驗和計算，指出在50～300 m/s的沖擊速度下，鳥撞可看作非恒定的流體動力學過程。總結了鳥撞過程主要分為初始撞擊、壓力衰減、恒定流動和流動結束。Zhang和Fei[15]采用SPH粒子法，對比了經典的半橢圓球形鳥體和真實鳥體對發動機風扇的沖擊響應，結果表明鳥體形狀和撞擊方向對鳥撞結果有很大的影響。

陳賀賀等[16]研究了4種鳥體本構模型對風擋的鳥撞，提出SPH粒子法和自定義本構鳥體的計算精度較高。潘春蛟等[17]對直升機主槳變距拉桿鳥撞進行了數值模擬，提出高速鳥撞造成變距拉桿功能失效的危害較大。朱貝蓓和蔡景[18]采用馬爾科夫鏈蒙特卡洛法，對中國運輸機鳥撞沖擊能量開展了研究，建立了鳥撞沖擊能量的概率分布曲線。陳靜和蔡景[19]提出了一種基于SPH粒子的鳥撞修正算法。張海洋等[20]使用PAM-CRASH軟件，對發動機吞鳥導致的葉片損傷進行了數值模擬，給出了葉片的主要損傷形式和部位。龐華華和韓全民[21]建立了全尺寸復合材料平尾有限元模型，通過對比試驗結果，對模型做出修正。幕琴琴等[22]研究了旋轉離心力對鳥撞葉片響應的影響。劉軍等[23-24]對鳥體本構模型進行了反演，研究表明鳥體與平板之間存在強烈的耦合效應，中心點最大變形是平板厚度的7.9倍，鳥體撞擊速度越高，流體特性越明顯。王計真和劉小川[25]提出了一種鳥體本構參數識別方法，并通過鳥撞試驗驗證了該方法的準確性。

綜上，已有的研究大都限于數值計算，鳥體本構以及簡單試件或結構的鳥撞試驗，分析區域多是鳥體的第一撞擊部位。對大型全尺寸結構以及鳥體穿透撞擊點后，對其他結構撞擊的分析研究較少。本文建立了全尺寸平尾模型，使用多歐拉域耦合的方法，在考慮空氣影響的前提下，分析了鳥撞對平尾整體結構的影響和前緣損傷，并與試驗結果進行了對比。

1 鳥撞試驗

根據CCAR29部第29.631要求，旋翼飛行器在速度等于最大極限速度(VNE)或者最大平飛速度(VH)(取較小者)受到1.0 kg的鳥擊后能夠繼續飛行或者安全著陸。根據要求，試驗的鳥彈質量為1.0 kg，撞擊速度為90.6 m/s。撞擊位置選擇在鳥撞擊平尾后能夠使平尾產生最大彎矩的位置。因此，撞擊點選擇在平尾最外側前緣2個小肋之間。平尾通過鋼制支撐結構與試驗平臺連接，試驗平臺裝置及工裝見圖1和圖2。試驗設備由發射系統、試驗靶架系統、速度測量系統和高速攝像系統組成，見圖3。

圖1 鳥撞試驗平臺示意圖Fig.1 Diagram of bird strike test platform

圖2 平尾鳥撞工裝Fig.2 Bird strike with flat tail worksuit

圖3 鳥撞試驗平臺Fig.3 Bird strike test platform

試驗中，將鳥體和彈托預先放入發射裝置中，待控制氣罐壓力達到預定值并穩定后，開啟壓力閥，高壓氣體從控制罐進入發射腔，推動鳥體和彈托在炮管內滑行，經過彈托分離機構分離鳥體和彈托[26]，鳥體撞擊到固定在試驗平臺上的平尾。

2 鳥撞數值模擬

2.1 鳥體本構模型

鳥撞擊靶結構的過程見圖4[27]，圖中：V為鳥體速度；Vx為激波速度。鳥體撞擊靶體過程會經歷初始撞擊、壓力衰減、恒定流動和流動結束4個階段，根據這種鳥撞特點，本文采用了文獻[14]的經典鳥體本構模型，鳥體使用圓柱體模擬，長徑比2∶1，鳥體參數見表1。空氣參數見表2，在鳥撞沖擊中，采用ESPOL多項式狀態方程模擬鳥體：

(1)

式中：p為壓力;e為單位質量的比內能;ρ為密度;ρ0為參考密度;a1～a6和b0～b1為自定義系數;μ=η-1，η=ρ/ρ0。

空氣采用γ律狀態方程定義：

p=(γ-1)ρe

(2)

式中：γ為比熱比。

圖4 鳥撞擊靶體的4個階段[27]Fig.4 Four stages of bird impacting target[27]

表1 鳥體參數
Table 1 Bird body parameters

參數密度/(kg·m-3)彈性模量/MPa直徑/mm長度/mm速度/(m·s-1)數值9302 2008817690.6

表2 空氣參數Table 2 Air parameters

2.2 平尾有限元模型

平尾結構主要由小翼和平尾組成，通過接頭與機體連接，如圖5所示，平尾總體尺寸見圖6。前緣和蒙皮是2024鋁蜂窩結構，前后梁為機加7050鋁梁，肋是2024鈑金件鋁板，各結構材料見表3。

圖5 平尾結構數模Fig.5 Numerical module of flat tail structure

圖6 平尾尺寸Fig.6 Flat tail size

表3 結構材料Table 3 Structural materials

名稱結構材料前/后緣蒙皮2024盒段[2024/蜂窩/2024]前/后梁7050肋腹板2024

對于金屬，采用各向同性彈塑性本構模擬(DMATEP卡片)，該本構還能夠模擬雙線性或分段線性的彈塑性行為，本文通過輸入最大塑性應變值來控制單元失效，材料本構關系見圖7,圖中:E為彈性模量；σs為屈服極限。

紙蜂窩采用正交各向異性本構，使用蔡-胡失效準則判斷失效：

(3)

圖7 彈塑性本構Fig.7 Elastoplastic constitutive

表4 鋁材力學性能[28]Table 4 Mechanics properties of aluminum[28]

表5 蜂窩力學性能[29]Table 5 Mechanics properties of honeycomb[29]

平尾結構主要使用SHELL元模擬，共計26 000個單元，剖面圖見圖8。

由于圖5中的小翼和接頭對鳥撞影響較小，為了保證計算效率，在有限元模型中，平尾使用半對稱模型，小翼與連接接頭進行了省略，接頭與平尾結構的連接使用約束X、Y、Z這3個方向平動來模擬，見圖9。

圖8 平尾有限元模型剖視圖Fig.8 Section of finite element model with flat tail

圖9 平尾有限元模型Fig.9 Finite element model of flat tail

2.3 流固耦合算法

流固耦合算法的目的是使歐拉網格與拉格朗日網格之間產生相互作用。兩部分域各自的控制方程通過耦合面聯立起來成為整個系統的控制方程組，這個方程組的求解是很困難的。所以在時間域上，程序采用數值方法，將整個時間段分解為一系列微小時間步，對2個系統分別迭代求解，求解時不考慮另一方狀態的變化。

例如：先假定拉格朗日域在tn時刻的狀態已知，在tn～tn+1時間步內，拉格朗日域的狀態不發生變化。在固定邊界條件下，對歐拉域求解，再將歐拉域內與耦合面相鄰單元的壓力作為載荷施加給拉格朗日單元。通過這種近似值的解法，對2個域輪流求解[30]。

由于鳥體在沖擊平尾前緣時，鳥體碎片在穿透平尾前緣后，會對前梁結構造成二次沖擊，故本文使用多歐拉域耦合(圖10)的方法，對沖擊過程進行模擬。在鳥撞的初始狀態的定義中，平尾外部的空間域使用六面體歐拉網格定義，鳥體形狀使用SHAPE卡片定義為圓柱體，前緣蒙皮和前梁之間的封閉空間使用自適應網格ADAPT卡片定義。當鳥體穿透前緣后，飛濺的碎片會通過內歐拉域網格繼續與平尾結構發生作用，2個流體歐拉域使用空氣定義，空氣與鳥體與平尾結構存在的耦合關系見圖10。

圖10 鳥撞歐拉域示意圖Fig.10 Diagram of bird strik Euler domain

3 計算結果討論

3.1 結構損傷

鳥撞過程試驗與仿真過程見圖11和圖12。在撞擊中，鳥體穿透前緣蒙皮，對前梁造成二次鳥撞。鳥體逐漸破碎成細小的骨頭和肉沫，其動態破壞過程表現為流體，因此表明用流體本構模擬鳥體是合理的；在鳥撞中不僅需要考慮撞擊部位，還需要考慮到鳥撞發生后，鳥體是否會穿透撞擊點對其他結構造成二次鳥撞。

圖11 鳥撞試驗錄像截圖Fig.11 Screenshot of bird strike test video

圖12 鳥撞仿真過程Fig. 12 Bird strike simulation process

平尾損傷對比見圖13，前緣蒙皮在鳥撞中出現了撕裂，前梁在鳥撞后出現凹坑，鳥撞點周圍的蒙皮與肋出現嚴重變形。仿真與試驗的失效形式基本吻合，由于在仿真中失效的單元會被刪除，所以只能通過計算得到結構的失效位置和塑性變形，對于鳥撞后結構的損傷形狀和撕裂形式模擬還有一定缺陷。

圖13 鳥撞仿真與試驗損傷對比Fig.13 Comparison of damage of bird strike simulation and test

3.2 鳥撞載荷

Barber和Boehman[31]通過試驗研究了鳥體和剛體靶板的正撞擊，測量了沖擊載荷時程曲線，發現鳥撞載荷不僅與鳥體質量、鳥撞速度有關，與撞擊結構也有關系，根據被撞擊結構的不同響應，給出了理想剛體模型、局部剛體模型和局部變形模型。

理想剛體模型假設被撞擊結構為剛體，在鳥撞中不發生整體的平動與轉動；局部剛體模型假設撞擊中有一定的平動與轉動；局部變形模型假設被撞擊結構在撞擊中存在平動和轉動，結構局部變形形成一個凹坑，鳥體陷入凹坑中，形成貝努利滯流壓力。

通過計算發現，鳥體在沖擊結構瞬間的載荷作用時長很短，在3 ms內就已完成對前緣蒙皮和前梁的沖擊，見圖14，鳥撞的位移測量點見圖15，結構的位移響應對比見圖16，從圖中可知由于應力波傳播速度的影響，在10 s后才會出現明顯變形。

由圖14可以看出，鳥撞前緣后，由于蒙皮被撕裂，載荷迅速衰減。鳥體沖擊前梁時，載荷出現峰值，之后進入平緩區，逐漸衰減為0。

圖14 鳥撞載荷隨時間變化Fig.14 Variation of bird strike load with time

圖15 位移測量位置Fig.15 Position of measurement of displacement

圖16 仿真與試驗變形值對比Fig.16 Comparison of deformation of simulation and test

這表明在第1次沖擊瞬間，除沖擊點外，整體結構并未發生明顯位移變形。所以本文使用理想剛體模型的理論計算公式(見式(4)),假設在第1次沖擊的瞬間，結構未發生整體變形，并與仿真值進行對比，結果見表6。

F=ρ1AV2sina

(4)

式中：ρ1為鳥體密度;A為鳥體截面積;V為鳥撞速度;a為鳥撞角度。F=930×0.0442×3.14×90.62×sin 90°=46 406 N

表6 沖擊載荷仿真與理論值對比

以上分析表明，在鳥撞中，受到沖擊失效的結構沖擊載荷持續時間較短。而未發生失效破壞的結構，會受到連續性沖擊載荷的影響，持續時間較長，危害性更大。

在此基礎上，本文討論了不同鳥撞速度下撞擊載荷與平尾結構損傷情況，不同鳥撞速度下的撞擊載荷和損傷對比如圖17所示。

從圖17可以看出，鳥體速度大于200 m/s時，平尾前梁2個端肋之間的梁腹板和上下蒙皮出現大面積損傷，2個端肋也出現了不同程度的破壞，肋的主要破壞部位在減輕孔周圍，根部結構沒有出現明顯的損傷，但局部蒙皮結構出現了屈服現象，主要屈服位置是平尾前梁部位。不同鳥撞速度下的沖擊載荷與損傷情況統計對比見表7和表8。

通過表7的對比可以發現，鳥體速度在小于50 m/s時對平尾結構的沖擊沒有出現明顯損傷，但沖擊載荷會隨速度的提高明顯增大。經驗公式(4)的計算值明顯小于數值模擬值，表明在小于50 m/s的低速沖擊中，鳥體本構沒有出現明顯的流體性質，與文獻[14]的結論一致。

在大于50 m/s的沖擊中，沖擊載荷的數值模擬結果與經驗公式(4)的計算結果出現了一定的吻合趨勢，50 m/s的單次數值模擬載荷大于經驗公式計算值；90 m/s的單次鳥撞載荷與經驗公式值基本吻合；150 m/s和200 m/s的前梁沖擊載荷與經驗公式值基本吻合；250 m/s的總鳥撞載荷與經驗公式值基本吻合；300 m/s的數值模擬值小于經驗公式值；以上現象表明經驗公式在平尾鳥撞分析中受鳥撞速度影響，當速度小于200 m/s時，式(4)可用于對前梁結構的二次鳥撞載荷評估，當>200 m/s速度<300 m/s時可用于評估平尾鳥撞總載荷，當速度≥300 m/s，式(4)計算值較保守。

圖17 速度20～300 m/s鳥撞載荷和損傷Fig.17 Load and damage at 20-300 m/s bird strike velocity

從表7、表8的對比結果可以看出：

1) 鳥撞總載荷具有明顯的非線性特點，在本文沒有考慮材料在高速沖擊中的硬化和材料流動狀態的情況下，總載荷并未隨著速度的增大而線性增長。

表7 不同鳥撞速度下的沖擊載荷與損傷情況Table 7 Impact load and damage at different bird strike velocities

表8 不同鳥撞速度下的數值計算鳥撞載荷比值Table 8 Ratio of bird strike loads at different bird strike velocities for simulation

2) 在某些速度下，結構越“剛硬”，在鳥撞中產生的載荷越大。載荷峰值大小與受沖擊的結構厚度有關，如本文中前梁厚度為1.5 mm，前緣蒙皮厚為前梁的40%，150～250 m/s鳥撞速度的前梁沖擊載荷是前緣蒙皮的40%以上，但90 m/s和300 m/s鳥撞速度的載荷相當。

3) 對于第1鳥撞點前緣，在20～300 m/s的速度范圍內，鳥撞載荷比值與速度比值的倍數基本吻合，且呈比例遞增。

4) 對于第2鳥撞點前梁，其峰值大于第1鳥撞點，并具有非線性特點，當速度大于250 m/s時，載荷峰值出現衰減現象。

3.3 平尾根部響應

鳥撞沖擊載荷，對平尾根部會產生很大的彎矩。為了研究沖擊載荷對平尾根部的影響，本文在平尾根部上下蒙皮處位置粘貼了6個應變片，見圖18，并在有限元模型的相應位置選取單元，與試驗數值進行了對比。峰值誤差對比見表9,應變時程曲線對比見圖19～圖21。

圖18 平尾根部應變片位置Fig.18 Position of flat tail root strain gauge

表9 應變仿真與試驗值對比
Table 9 Comparison of simulation and test results of strain

應變片仿真值試驗值誤差/%E7922.8921.70.1E81 117.21 003.111E9922.8784.817

圖19 E7位置應變試驗與仿真值對比Fig.19 Comparison of E7 position strain between test and simulation values

圖20 E8位置應變試驗與仿真值對比Fig.20 Comparison of E8 position strain between test and simulation values

圖21 E9位置應變試驗與仿真值對比Fig.21 Comparison of E9 position strain between test and simulation values

從90 m/s速度的鳥撞載荷計算結果可以得到，在鳥撞產生了2次約為50 000 N左右的沖擊載荷情況下，鳥撞部位至根部長度約為1 500 mm，由此產生的彎矩為75 000 N/m，但從表9的根部應變值對比可以看出，鳥撞后的根部響應應變峰值不高，平尾根部未發生撕裂性破壞，這表明平尾結構在鳥撞過程中前緣蒙皮和前梁結構起到了很好的吸能作用，消耗掉了大部分的沖擊動能。

根據以上分析得到的鳥撞載荷特點，載荷時程曲線主要分為峰值區、平緩區和衰減區，見圖22。

圖22 鳥撞載荷模型Fig.22 Bird strike load model

4 結 論

1) 在平尾結構的鳥撞分析與試驗中，除了鳥撞部位，還應考慮沖擊載荷在平尾根部產生的附加彎矩對結構造成的影響，避免發生根部斷裂。

2) 多歐拉域耦合法模擬鳥撞，能夠準確地模擬鳥體對結構連續沖擊的過程。鳥撞通常在幾毫秒內完成，而結構響應由于應力波傳遞速率的影響，會出現明顯的遲滯效應。

3) 如果鳥體在沖破第一結構后仍然保持完整，對第二結構的沖擊仍會產生較大載荷，危害性較大。

4) 在平尾結構的抗鳥撞設計中，應預先對主要承力部件如：前梁、肋等結構的剛度值和材料選擇進行測試，避免因材料過強或過弱導致在高速沖擊中出現斷裂，確保選擇的材料在鳥撞中能夠發揮一定的抗沖擊和吸能作用。