AP1000安全殼在大型商用飛機撞擊下的動力響應分析

汪大洋， 汪秀清， 楊 強， 張永山， 吳成清

(1. 廣州大學 土木工程學院，廣州 510006; 2. 悉尼科技大學 土木與環境工程學院，悉尼 2007)

核能作為一種清潔且高效的重要能源，已經成為電力工業發展一種必然趨勢。然而，隨著恐怖襲擊的日益猖獗和現役飛機數量的迅速增加，核電廠正面臨越來越多的蓄意和意外飛機撞擊威脅。特別是“9·11”事件以來，大型商用飛機對核電站的撞擊成為關注焦點，我國國家核安全局發布HAF 102—2016《核動力廠設計安全規定》[1]明確要求核電站考慮商用飛機惡意撞擊。核電安全殼是核燃料反應堆最后一道安全屏障，其重要性不言而喻。

在試驗研究方面，1993年美國Sandia國家實驗室同日本Kobori研究所聯合進行F4“鬼怪式”戰斗機原型撞擊鋼筋混凝土厚板試驗[2-4]，得到飛機沖擊平面靶板的撞擊力時程曲線，將其同利用Riera模型的理論計算結果進行對比，驗證了Riera模型的正確性，該試驗也是全球目前唯一一個足尺模型試驗。孔建偉等[5]將CPR1000核電站安全殼作為研究對象，設計了1/20安全殼模型，采用16 kg混凝土塊模擬質量為13 t的新舟600客機，用質量為140 kg的鐵制圓筒模擬質量為40 t的波音737-800客機，研究了飛機模型撞擊下安全殼結構的動力響應和損傷破壞現象及其過程。

在仿真研究方面，左家紅等[6-7]采用ADINA軟件對秦山核電站進行了動力響應分析研究，討論了安全殼撞擊區域的非線性行為。Kukreja[8]對PHWR核電站的雙層安全殼進行了簡化非耦合數值模擬，結果表明該安全殼有足夠的承載能力來抵御B707-320和A300B4-200飛機撞擊。Iqbal等[9]和Sadique等[10]采用ABAQUS軟件建立了1.2 m厚鋼筋混凝土安全殼有限元模型，分別將5種不同飛機的撞擊荷載時程曲線直接加載到安全殼結構上，研究分析了結構的動力響應。曹健偉等[11-14]建立了新舟MA600飛機、空客A320飛機和預應力鋼束鋼筋混凝土安全殼有限元模型，開展了耦合與非耦合有限元仿真研究。本文課題組針對AP1000核電站結構的地震動性能展開了相關研究[15-18]，對其安全殼結構及其力學性能有較深入理解。

由于試驗復雜、費用巨大，試驗研究難以對大型商用飛機撞擊下的核安全殼的整體動力響應進行詳細分析。隨著有限元快速發展，數值模擬成為一條重要研究途徑，并取得了較多的成果。本文以波音737 MAX 8大型商用飛機和AP1000安全殼為研究對象，進行安全殼在不同初始撞擊速度和撞擊高度下的動力響應分析，研究安全殼的損傷破壞機理，其結果可為后續安全評估和結構設計提供參考。

1 有限元模型

1.1 AP1000安全殼

根據西屋電氣公司關于AP1000的設計文件和團隊現有AP1000安全殼結構[19]進行建模，AP1000安全殼是一個自由直立的圓柱形鋼制容器，帶有橢球形的上封頭和下封頭。鋼制安全殼容器被完全包容在一個混凝土屏蔽體中。直徑39.62 m，總高度65.63 m，其中圓柱形筒壁高度54.16 m(地下部分為12 m)，圓形穹頂高度11.47 m。圓柱形筒壁鋼板厚度44.45 mm，圓頂鋼板厚度41.27 mm，如圖1所示。

圖1 安全殼有限元模型及網格劃分(m)Fig.1 Finite element model of containment and meshing(m)

采用殼單元(SHELL163)仿真建模，單元厚度為結構各部分實際厚度。單元尺寸統一500 mm×500 mm，單元總數55 808。鋼制安全殼SA738-B型鋼材，密度7 870 kg/m3，彈性模量2.0×105MPa，泊松比0.2，屈服應力415 MPa，切線模量1.0×103MPa[20-21]。鋼材本構模型選用*MAT_PLASTIC_KINEMATIC，在關鍵字文件中打開應變率。在飛機撞擊安全殼過程中使用*CONTACT_AUTOMATIC_NODES_TO_SURFACE模擬飛機機身、機翼和引擎部件與安全殼之間的接觸，罰函數系數取1.0確保合理接觸。通過關鍵字文件*HOURGLASS在不同部件(*PART)上執行沙漏控制。

1.2 波音737 MAX 8飛機

波音 737 MAX 8[22]總質量約80 t，其中兩臺引擎質量約10 t。總長度約39 m，翼展長度約36 m，機身直徑約3.88 m，如圖2和圖3所示。由于文件資料有限，故現根據已經掌握的設計文件資料對飛機模型進行簡化：①幾何形狀與實際飛機相同；②考慮飛機蒙皮等主要結構，包括機身、機翼和引擎，所有尺寸均與實際相同，未考慮飛機內部連接和構造；③僅建立飛機結構模型，未考慮飛機燃料、貨物和乘員狀況。

采用殼單元(SHELL163)對飛機建模，采用MATLAB軟件計算飛機幾何模型機身、機翼及引擎的面積，再根據飛機材料密度，求得機身殼單元厚度約為32 mm，左右大翼殼單元厚度約為36 mm，引擎殼單元厚度取19 mm。Boeing 737 MAX 8飛機引擎為超級鋼材料，機身與機翼為鋁合金材料。飛機材料本構模型均采用*MAT_PLASTIC_KINEMATIC，考慮了材料失效、塑性流動，模擬金屬材料在撞擊過程中的動力響應，其具體的材料參數如表1所示。根據文獻[23-24]對材料失效應變分析，確定機身鋁合金材料的失效應變為0.05，機翼鋁合金為0.01，引擎超級鋼為0.1。單元尺寸統一為500 mm×500 mm，單元總數3 459個，其中鋁合金材料單元個數3 075，超級鋼384個。飛機各部件之間采用共節點處理。

圖2 波音737外形及幾何尺寸(cm)Fig.2 Shape and geometric of Boeing 737 (cm)

圖3 波音737 MAX 8Fig.3 Boeing 737 MAX 8

表1 飛機材料參數Tab.1 Aircraft material parameters

2 基于Riera法的飛機模型驗證

首先通過Riera法計算得到大型商用飛機波音737 MAX 8撞擊剛性墻的撞擊力時程曲線和沖量時程曲線，然后利用ANSYS/LS-DYNA建立飛機和剛性墻模型并進行撞擊過程仿真，最后對比分析撞擊力的理論數值模擬結果來驗證飛機有限元模型的準確性。

2.1 Riera法計算撞擊力

設t=0開始撞擊剛性墻，在t時刻速度為V(t)，質量為M(t)，機頭到壓碎位置距離X(t)，飛機質量密度μ[X(t)]。由動量定理可得撞擊力F(t)計算公式

(1)

假設Pc(X)為飛機壓碎力，在撞擊處由于假定飛機材料為剛塑性，因此撞擊處材料一維桿內力為Pc(X)，并對未撞碎區域取隔離體，由牛頓第二定律可得

Pc(X)=M(t)a(t)

(2)

因此Riera法可寫成

P(t)=Pc[X(t)]+μ[X(t)]V2(t)

(3)

式中：P(t)飛機總撞擊力；Pc[X(t)]飛機在該位置壓碎力；μ[X(t)]飛機單位長度質量，不考慮撞擊過程質量損失；V(t)飛機未壓碎部分剩余速度。

根據文獻[25]，飛機初始撞擊速度取200 m/s，μ為1。Riera法理論計算出的撞擊力曲線，如圖4所示。

圖4 Riera法理論計算結果Fig.4 Theoretical calculation result of Riera

2.2 基于Riera法的飛機模型驗證

飛機初始撞擊速度200 m/s，飛機垂直撞擊剛性墻。有限元模型如圖5所示，在墻體中心位置進行撞擊，考慮飛機翼展長度和撞擊區域的長度，墻體長度和寬度均取50 m，厚3.5 m，墻體四周固定約束。采用實體單元SOLID164建模。在墻體長和寬度方向上，單元尺寸取500 mm，厚度方向等分10份。墻體為素混凝土，在ANSYS/LS-DYNA軟件中采用關鍵字文件*DEFORMABLE_TO_RIGID將墻體定義為剛性墻。

圖5 飛機撞擊剛性墻有限元模型(m)Fig.5 Finite element model of impacting the rigid wall(m)

Riera法理論和仿真結果對比，如圖6所示。對比分析可得，兩者的飛機撞擊力時程曲線和沖量時程曲線的趨勢基本一致，且飛機撞擊力峰值最大誤差為1.11%，飛機沖量最大誤差為3.73%。可知，所建模型和相關參數設計是合理的，可以用于后續分析。

圖6 Riera法理論計算結果與數值模擬結果比較Fig.6 Results comparison between Riera’s theoretical calculation and numerical simulation

3 撞擊計算結果及分析

根據目前民航飛機在起飛和著陸的所有速度及現有飛機墜毀記錄，飛行速度約100 m/s。波音737 MAX 8巡航速度約233 m/s。考慮到飛機墜毀、劫機或恐怖襲擊時的蓄意撞擊，將飛機初始撞擊速度設為100 m/s，150 m/s，200 m/s，250 m/s和300 m/s。同時為研究撞擊高度對安全殼動力響應的影響，撞擊高度設為P1(30 m)、P2(39 m)、P3(47 m)、P4(54 m)、P5(65 m)，其中：P1(30 m)和P3(47 m)分別為鋼安全殼內部等效鋼梁筒身處，P2(39 m)為兩個等效鋼梁之間的非等效鋼梁筒身處，P4(54 m)為安全殼筒身與穹頂交界處，P5(65 m)為安全殼穹頂位置，如圖7所示。飛機撞擊工況表如表2所示。撞擊全過程分4個階段：第一階段，機頭和前部機身；第二階段，中部機身、引擎和機翼；第三階段，中部機身和機翼；第四階段，后部機身和尾翼。依次定義為T1、T2、T3、T4，同時將撞擊末時刻定義為T5。為節省計算時間，僅考慮計算單側撞擊面的撞擊作用，由于撞擊速度不同，飛機撞擊單側面的作用時間依次定義為0.40 s，0.35 s，0.25 s,0.20 s,0.15 s。

圖7 撞擊位置及撞擊階段示意圖Fig.7 Diagram of impact position and impact stage

表2 工況表Tab.2 Investigation cases

3.1 飛機動力響應分析

3.1.1 飛機撞擊破壞現象

由圖8(a)可以看出P3V1工況中的飛機撞擊破壞現象。在t=0.08 s時，飛機機頭撞擊到安全殼壓屈破壞，后面的主要機身以及引擎和機翼未發生明顯的變形，機身未進入安全殼內部；在t=0.16 s時，飛機引擎撞擊到安全殼，機身后面的部位未發生明顯的變形，飛機機身亦未進入安全殼內部；在t=0.24 s時，飛機的進一步破壞，引擎變形加大且已與機翼分離，靠近機身的機翼根部已經完全壓屈破壞，飛機機身未進入安全殼內部；在t=0.32 s時，飛機剩余的機身未進一步破壞，但靠近機身的機翼根部已經完全破壞，飛機機翼尾的尾部發生斷裂飛出，分布在安全殼兩側，引擎破壞情況沒有明顯加劇，但是其發生翻轉，剩余飛機機身未進入安全殼內部；在t=0.40 s時，飛機剩余的機身仍未進一步發生破壞且未進入安全殼內部，剩余的引擎和機翼部分在安全殼外部發生了反彈；由圖8(b)可以看出，P4V1工況中的飛機撞擊破壞現象，在t=0.08 s時，機頭已經進入安全殼內部且機身未發生明顯變形，從t=0.16～0.32 s，機身繼續進入安全殼，此時飛機引擎和機翼也完成了撞擊，引擎完全破壞，靠近機身的機翼機翼根部進入安全殼內部，機翼尾部沿安全殼兩側飛行；由圖8(c)可以看出P5V1工況中的飛機撞擊破壞現象。在t=0.08 s,t=0.16 s和t=0.24 s時，飛機未發生明顯變形，都是從安全殼穹頂外側飛過，未進入安全殼內部；在t=0.32 s時，靠近引擎的部分機翼發生變形，其他部位未發生明顯變形，飛機仍未進入安全殼內部；在t=0.40 s時，飛機引擎發生變形，機翼的變形加大，飛機機身未見明顯變形，飛機未進入安全殼內部，以上三種工況的飛機破壞現象為飛機撞擊的典型現象。

3.1.2 飛機撞擊力

飛機撞擊力時程曲線如圖9(a)～圖9(e)所示。可知，在整個撞擊過程中，撞擊力呈先增大后減小趨勢。在不同撞擊速度下整個過程都可分為4個階段，①T1-T2撞擊階段，撞擊作用較為穩定，撞擊力在某個值附近波動；②T2-T3撞擊階段，撞擊力突然顯著增大，圖中曲線的斜率接近90°，短時間內達到撞擊力峰值，各個初始撞擊速度下等效鋼梁筒身處和非等效鋼梁筒身處撞擊力峰值相較于第一階段都增大3～4倍；③T3-T4當引擎完成撞擊，僅有中部機身和剩余機翼與安全殼發生撞擊，撞擊力明顯減小；④T4-T5撞擊力曲線保持穩定波動，該時間段內與安全殼發生撞擊的飛機部位為僅為后部機身。值得一提的是，P5位置處飛機與安全殼之間接觸非垂直作用，飛機在撞擊安全殼后飛行方向發生了向頂上偏移，與安全殼間的撞擊作用較小，因此沒有上述現象的發生。

圖9(f)給出了所有工況下的飛機撞擊力峰值。由圖可知，各個位置下撞擊力峰值隨著速度的增加呈線性增加，對于筒身來說，P1位置撞擊力峰值增加幅度為20%,23%,38%,73%，P2位置處撞擊力峰值增加幅度為71%,78%,106%,150%。不同初始撞擊速度下，筒身等效鋼梁和非等效鋼梁處的撞擊力峰值差異有所不同，當速度為150 m/s時，P3處撞擊力峰值比P2處大17%，在初始撞擊速度為300 m/s時，筒身等效鋼梁處P1和P3撞擊峰值最大值達到了221.2 MN較筒身非等效鋼梁P2處大30%，P4處大71%，P2和P4處的差異是由于P4位置處于筒身和穹頂的交界處，當飛機撞擊時接觸情況不同所導致。

圖8 飛機撞擊破壞現象Fig.8 Aircraft crash damage

圖9 飛機撞擊力時程曲線及撞擊力峰值比較Fig.9 Aircraft impact force-time history curve and comparison of maximum impact force

3.1.3 飛機動能分析

飛機動能時程曲線如圖10(a)～圖10(e)所示，對應不同初始撞擊速度和撞擊高度下的撞擊末時刻的動能比較如圖10(f)所示。由圖可知，對于筒身來說，當初始撞擊速度相同時，筒身安全殼鋼梁位置(P1和P3)處的動能變化曲線幾乎重合，在速度小于或等于200 m/s時，動能在撞擊末時刻都幾乎降為0，且動能損失速度要明顯快于非鋼梁位置(P2和P4)，P4位置在速度達到150 m/s時，末時刻仍有3×108N·m動能，當初始撞擊速度達到300 m/s時，筒身等效鋼梁處兩條曲線幾乎重合仍保持較快的動能損失，但筒身非等效鋼梁處幾乎和穹頂P5位置的曲線重合，說明該速度下筒身非等效鋼梁處對飛機的阻礙作用非常小，以上都說明等效鋼梁對于飛機的撞擊有較好的抵抗效果。對于穹頂來說，飛機與安全殼接觸面之間非垂直作用，與安全殼間的撞擊作用較小，飛機有一定的速度損失，但在撞擊末時刻仍能夠離開安全殼上空繼續飛行。由圖10(f)撞擊末時刻T5動能隨初始撞擊速度的變化情況也可以得到，筒身等效鋼梁處對飛機的抵抗效果是最好的，在速度為300 m/s時，T5時刻飛機撞擊筒身等效鋼梁處的動能僅為其他撞擊位置處的1/2。

圖10 飛機動能時程曲線及T5時刻動能比較Fig.10 Aircraft kinetic energy-time history curve and kinetic energy comparison at T5

3.2 安全殼位移響應及撞擊區域破壞情況分析

3.2.1 位移響應分析

圖11給出了不同初始撞擊速度和撞擊高度下，飛機撞擊面沿安全殼高度分布的飛機撞擊方向上水平位移峰值。由圖11可知，在不同初始撞擊速度和撞擊高度下AP1000安全殼正面外壁水平向位移峰值均在各工況的飛機撞擊所在位置處出現最大值，且隨著初始撞擊速度的增加而增大，其上下相鄰高度處的正面外壁水平向位移峰值也均大于其他位置處，隨著遠離撞擊中心，撞擊位移逐漸衰減。在初始撞擊速度小于等于150 m/s時，安全殼筒身等效鋼梁部位P1(30 m)和P3(47 m)的撞擊位移較小，最大只有1.62 m，筒身非等效鋼梁處P2(39 m)和P4(54 m)的撞擊位移較大，最大達到6.1 m，為前者的3.8倍。而當速度增加到200 m/s，情況正好相反，等效鋼梁處P1和P3的撞擊位移發生激增，30 m處的撞擊位移最大達到了12.6 m。撞擊穹頂處所有計算工況下的安全殼正面外壁水平向位移峰值在結構高度60 m及以內保持在很小值且均小于其他所有計算工況，可見此時安全殼撞擊位置以下的安全殼筒體外壁沿飛機撞擊方向未發生明顯位移。

3.2.2 撞擊區域破壞情況分析

安全殼在撞擊區域的貫穿尺寸是安全殼動力響應和破壞情況的最直觀反映。安全殼的具體環向和豎向貫穿尺寸結果，分別如表3和表4所示。部分工況撞擊的安全殼撞擊區域破壞情況圖示，如圖12所示。結果分析表明，在撞擊高度為65 m的所有工況下，飛機撞擊安全殼穹頂，穹頂外壁發生變形，但其均未發生貫穿。除上述未貫穿的計算工況外，飛機初始撞擊速度為100 m/s時，安全殼等效鋼梁位置處，安全殼的外壁鋼板均發生內凹變形，但均未發生貫穿破壞，其余計算工況均發生貫穿破壞，且當撞擊高度為54 m時，安全殼環向和豎向貫穿尺寸均大于同初始撞擊速度下的尺寸，最大環向貫穿尺寸達到了29.68 m，最大豎向貫穿尺寸達到了17.86 m，可知鋼制安全殼的最危險位置為安全殼圓筒形外壁與其穹頂的交界處，即高度為54 m處。根據不同速度下筒身等效鋼梁位置和非等效鋼梁位置的撞擊破壞現象可知，筒身非等效鋼梁處在速度為100 m/s時，就會發生貫穿破壞，環向貫穿尺寸為3.95，豎向貫穿尺寸為5.95，而等效鋼梁處在速度達到150 m/s時，才會發生貫穿，且環向貫穿尺寸為2.71，豎向貫穿尺寸只有2.73，貫穿尺寸為前者的1/3；當初始撞擊速度從100 m/s增加到250 m/s時，安全殼的環、豎向的貫穿尺寸都隨之增加，在初始撞擊速度為250 m/s時達到最大，但當初始撞擊速度增加到300 m/s時，貫穿尺寸反而減小。

圖11 AP1000安全殼正面沿撞擊方向的位移峰值Fig.11 Peak displacement of AP1000 containment along the impact direction

圖12 T5時刻的安全殼破壞情況圖示Fig.12 Diagram of containment damage at T5

表3 安全殼環向貫穿尺寸Tab.3 Circumferential dimensions of containment

表4 安全殼豎向貫穿尺寸Tab.4 Vertical penetration dimension of containment

3.3 安全殼等效應力分布

圖13為AP1000核島結構安全殼P3和P4在不同初始撞擊速度下的等效應力云圖，選擇這兩個位置的原因是，P3為等效鋼梁位置，P4是非等效鋼梁位置，可以比較兩者的應力發展情況，且由前面的分析可知，P4是貫穿破壞最嚴重的區域。為了更加明確和比較鋼制安全殼的應力集中現象和安全殼的失效，應力云圖統一為最大應力邊界條件。觀察圖12中AP1000安全殼外壁等效應力云圖可以看出，在不同撞擊位置處，安全殼的等效應力主要集中在正面區域，安全殼正面的損傷破壞區域形狀基本相同，等效最大應力最開始出現于撞擊中心區域，機身撞擊區域呈圓形，隨著撞擊的繼續進行，安全殼上最大等效應力范圍先是沿高度方向擴展，當引擎和機翼撞擊安全殼時，安全殼上最大等效應力范圍開始沿水平方向擴展，且撞擊中心區域最大等效應力分布范圍增大。等效應力的最大值都是出現在安全殼撞擊中心區域處，且此處外壁鋼板產生的最大等效應力均已超過安全殼鋼材的屈服應力415 MPa。由等效應力圖的發展情況來看，當初始撞擊速度相同時，P3位置等效應力影響范圍要明顯大于P4位置，P4的撞擊等效應力更為集中，在t2和t3時刻，P3位置對整個安全殼正面都有一定的影響，在P1位置處都有一定的鋼材達到屈服應力，而P4位置都集中在撞擊點附近。

圖13 安全殼等效應力云圖Fig.13 Equivalent stress contours of the containment

4 結 論

基于已驗證的波音737 MAX 8大型商用飛機和AP1000安全殼，采用ANSYS/LS-DYNA軟件進行了25種工況下飛機撞擊安全殼全過程數值模擬，分析了飛機和安全殼結構的動力響應、安全殼的局部破壞和等效應力分布。結論如下：

(1)基于彈-靶體接觸耦合分析方法能有效預測飛機撞擊動力過程與效應，與Riera法理論計算結果相比，飛機撞擊力峰值和沖量的最大誤差僅為1.11%，3.73%。

(2)飛機撞擊力時程曲線先增大至峰值后減小，其峰值出現在第二階段，且該階段的整個撞擊力主要由引擎撞擊貢獻，撞擊力峰值相較于第一階段增大3～4倍。

(3)不同初始撞擊速度下，筒身等效鋼梁處的撞擊力峰值較非等效鋼梁處大，最小為后者的117%(速度為150 m/s)，最大達到了后者的171%(速度為300 m/s)。

(4)等效鋼梁能夠很好地抵御飛機的撞擊作用。在初始撞擊速度為100 m/s時，筒身非等效鋼梁處已經貫穿破壞，環向貫穿尺寸為3.95 m，豎向貫穿尺寸為5.95 m，而等效鋼梁處在初始撞擊速度達到150 m/s時才發生貫穿破壞，且環向、豎向貫穿尺寸分別為2.71 m，2.73 m，貫穿尺寸僅約為前者的1/3，且撞擊等效鋼梁后飛機動能下降最快，在T5時刻飛機撞擊筒身等效鋼梁處的動能僅為其他撞擊位置處的1/2(初始撞擊速度為300 m/s)。

(5)安全殼筒身段與穹頂的交界處為安全殼結構的最危險位置(54 m高度處)，在各個初始撞擊速度下都發生了貫穿破壞，且在此位置處安全殼環向和豎向貫穿尺寸均大于同初始撞擊速度下其他位置的貫穿尺寸，最大環向、豎向貫穿尺寸達到了29.68 m，17.86 m。安全殼穹頂在所有工況的飛機撞擊下均未發生貫穿破壞(撞擊高度為65 m)，撞擊使飛機的飛行方向發生變化，均向安全殼斜上方飛出。

(6)當飛機撞擊速度超過150 m/s時，安全殼撞擊的鋼板等效應力影響范圍隨初始撞擊速度的增加而減少，且撞擊等效鋼梁處的安全殼等效應力分布范圍相對于撞擊非等效鋼梁處的更大。