基于土-結(jié)構(gòu)相互作用的核電廠房抗飛機(jī)撞擊耦合分析與模型驗證

王友剛， 孫運輪， 李建波， 李曉紅， 陳 巖， 梅潤雨， 牛燕如， 林 皋

(1. 中核能源科技有限公司，北京 100193；2. 大連理工大學(xué) 海岸與近海工程國家重點實驗室，遼寧 大連 116024；3. 大連理工大學(xué) 工程抗震研究所，遼寧 大連 116024)

發(fā)展核電是國家重要戰(zhàn)略需求，確保核電安全是未來內(nèi)陸核電發(fā)展的先決性條件[1]。受“911”事件的重大啟示，大型飛機(jī)撞擊安全殼及其可能的嚴(yán)重后果日益受到人們的關(guān)注，尤其近年來各類大型飛機(jī)在航空業(yè)應(yīng)用廣泛，翼展達(dá)數(shù)十米，與安全殼直徑幾乎相當(dāng)[2]，在以100～200 m/s的速度產(chǎn)生的巨大沖量作用下，可能對安全殼產(chǎn)生極大的沖擊荷載。國內(nèi)外相繼出臺各項法規(guī)[3-4]，正式將抵御大型商用飛機(jī)惡意撞擊列入核電站建造的考慮范圍。

飛機(jī)撞擊安全殼的研究最早開始于1968年Riera公式[5]，主要是聚焦在不同因素對飛機(jī)撞擊荷載曲線的影響。隨著數(shù)值分析技術(shù)的快速發(fā)展，基于飛射物-靶體相互作用的精細(xì)化分析方法，日益在核電結(jié)構(gòu)抗飛機(jī)撞擊中得到廣泛應(yīng)用[6]，這一領(lǐng)域研究的廣度和深度也在逐漸擴(kuò)展和完善[7]，如撞擊位置，安全殼墻體曲率，撞擊過程飛機(jī)速度衰減特征等。

Sadique等[8]使用4種不同型號飛機(jī)撞擊安全殼，并對其最大變形和整體性能進(jìn)行了對比；Jeon等[9]對采用纖維增強(qiáng)混凝土的核安全殼結(jié)構(gòu)抗沖擊性能進(jìn)行了分析；黃濤等[10]建立了精細(xì)化的大型商用飛機(jī)與核電站安全殼(廠房)有限元模型，并對飛機(jī)撞擊規(guī)律與特點進(jìn)行了總結(jié)；Zou等[11]提出一種將比例邊界有限元和八叉樹技術(shù)相結(jié)合的跨尺度分析方法，對飛機(jī)撞擊核安全殼問題進(jìn)行了精確且高效的模擬。

安全殼等核電廠房結(jié)構(gòu)在飛機(jī)撞擊作用下的動力響應(yīng)，是一種短時、多因素影響的強(qiáng)耦合過程。與采用簡單的撞擊荷載曲線方法相比，精細(xì)化的相互作用分析方法具有顯著的優(yōu)勢[12]，尤其是便于從飛機(jī)-廠房-地基耦合體系的角度開展機(jī)理模擬。內(nèi)陸非巖基廠址條件下，地基對上部結(jié)構(gòu)的剛度約束作用大幅減弱，大飛機(jī)撞擊這一沖擊載荷的耗散過程及其對安全殼-地基相互作用體系的整體動力影響是值得關(guān)注的重要課題，對這一問題的系統(tǒng)性分析，將對核電結(jié)構(gòu)抗飛機(jī)撞擊領(lǐng)域的發(fā)展起到積極的促進(jìn)作用。

土-結(jié)構(gòu)相互作用框架下，飛機(jī)撞擊和地震載荷所采用的地基動力模型，具有明顯的異同，飛機(jī)撞擊對于地基外邊界來說是典型的內(nèi)源動力問題，而地震則是外源動力問題[13-15]，如圖1、圖2所示。從地基中散射波動的單向傳播角度看，兩者地基外邊界的透射條件要求本質(zhì)上是一致的。

圖1 地震荷載

圖2 飛機(jī)撞擊荷載

相互作用體系抗大飛機(jī)撞擊的數(shù)值模擬主要涉及以下幾方面問題：上部精細(xì)化飛機(jī)模型和安全殼模型的建立；相互作用模型的建立；非線性材料本構(gòu)的選擇與參數(shù)調(diào)整；沖擊破壞全過程的數(shù)值模擬。上述數(shù)值模型較為復(fù)雜，片面追求精細(xì)化的有限元模型在提高分析精度的同時，也會帶來工作量的巨大壓力，在保證精度和效率的雙重要求下，如何有效平衡兩者，是本文的關(guān)注點之一。

目前從相互作用體系角度探討核電廠房結(jié)構(gòu)的抗飛機(jī)撞擊過程仍缺乏標(biāo)準(zhǔn)化的分析方法。本文提供了一種考慮相互作用條件下飛機(jī)撞擊分析的方法，主要聚焦在精細(xì)模型建立及參數(shù)的有效性驗證、可模擬無限地基輻射阻尼效應(yīng)的人工邊界設(shè)置、以及不同軟硬程度土體對上部結(jié)構(gòu)振動響應(yīng)的影響，并基于此建立了飛機(jī)-廠房-地基全耦合分析系統(tǒng)。

1 有限元模型

1.1 飛機(jī)模型

1.1.1 飛機(jī)模型參數(shù)

本文飛機(jī)模型使用的是某四引擎大型客機(jī)，如圖3所示。機(jī)身長63.6 m，翼展60.3 m，總質(zhì)量為234 t，由190 290個單元組成。在建模時使用MASS166單元為其增加配重，使之滿足真實的質(zhì)量分布。

圖3 飛機(jī)有限元模型

飛機(jī)蒙皮(機(jī)身、機(jī)翼部位)為航空鋁材料，引擎的材料是剛，均使用Jackson-Cook本構(gòu)模型，飛機(jī)內(nèi)部的鋼筋仍考慮采用與核島廠房鋼筋相同的本構(gòu)，材料參數(shù)參考文獻(xiàn)[16]。

1.1.2 飛機(jī)模型驗證

為確保用于計算的飛機(jī)模型真實有效，通過比較采用數(shù)值方法計算的撞擊剛性靶板的撞擊力和采用Riera方法計算的撞擊力對飛機(jī)模型的合理性進(jìn)行比較。分別建立網(wǎng)格尺寸差異明顯的兩種網(wǎng)格模型，對其網(wǎng)格敏感性進(jìn)行分析。

本文分別使用不同大小的網(wǎng)格對模型進(jìn)行剖分，模型1主要采用的殼單元網(wǎng)格尺寸約為(280×155)mm2，梁單元網(wǎng)格尺寸約為155 mm，該模型的節(jié)點數(shù)為110 601，單元數(shù)為190 290，如圖4所示。模型2主要采用的殼單元網(wǎng)格尺寸約為(140×75)mm2，梁單元網(wǎng)格尺寸約為75 mm，該模型的節(jié)點數(shù)為310 663，單元數(shù)為681 586，如圖5所示。

圖4 模型1網(wǎng)格

圖5 模型2網(wǎng)格

采用數(shù)值方法進(jìn)行飛機(jī)撞擊剛性靶板的撞擊力計算，是通過對飛機(jī)進(jìn)行精細(xì)化的有限元建模，以獲得真實的飛機(jī)撞擊剛性靶板過程。在進(jìn)行分析時，假定飛機(jī)以200 m/s的速度垂直撞擊剛性靶板，如圖6所示。

圖6 飛機(jī)垂直撞擊剛性靶板

通過以上兩種方法，分別得到飛機(jī)撞擊剛性靶板的撞擊力時程曲線和沖量時程曲線對比圖，如圖7、圖8所示。

圖7 撞擊力時程對比

圖8 撞擊力沖量時程對比

通過以上結(jié)果可以看出，模型1和模型2與Riera方法的結(jié)果吻合得很好，以此驗證了兩種有限元模型的合理性。同時，通過結(jié)果對比也可以看出，兩種模型的計算結(jié)果相差不大，然而模型1的計算時間約為4 h，模型2的計算時間約為20 h。因此，采用模型1的飛機(jī)模型進(jìn)行后續(xù)分析。

1.2 核島廠房模型

以某核電站示范項目為背景建立精細(xì)化的廠房模型。主要由反應(yīng)堆廠房、乏燃料廠房、核輔助廠房和電氣廠房構(gòu)成，如圖9(a)所示。反應(yīng)堆廠房細(xì)致模擬了外部安全殼結(jié)構(gòu)與內(nèi)部艙室結(jié)構(gòu)，如圖9(b)、圖9(c)所示。

核島廠房整體由混凝土澆筑，在安全殼穹頂內(nèi)側(cè)與艙室結(jié)構(gòu)外側(cè)有鋼板，安全殼殼體內(nèi)部布設(shè)有環(huán)向鋼筋與縱向拉筋，如圖9(d)所示。利用ANSYS/LS-DYNA軟件建立了如圖9(a)～圖9(d)所示有限元模型，其中混凝土使用SOLID164單元，鋼板使用SHELL163單元，鋼筋使用LINK160單元，鋼筋與混凝土采用分離式建模方法，并通過*CONSTRAINED_LAGRANGE_IN_SOLID關(guān)鍵字耦合在一起，單元尺寸為0.3 m，單元總數(shù)為4 232 520個。

圖9 有限元模型

1.2.1 核島廠房本構(gòu)模型參數(shù)

不同載荷作用下混凝土材料的率敏感性差異較大，如圖10所示。沖擊荷載下應(yīng)變率約100～101s-1。在動力學(xué)問題分析中，混凝土結(jié)構(gòu)的損傷及演化是十分復(fù)雜的過程，若還使用混凝土的靜態(tài)本構(gòu)則計算結(jié)果的可信度無法保證，可有效模擬材料的率相關(guān)動力特性是DYNA軟件的突出特征。基于LS-DYNA中豐富的材料庫，考慮使用混凝土連續(xù)面蓋帽模型進(jìn)行飛機(jī)撞擊問題的分析。該模型考慮了材料的硬化、損傷以及率相關(guān)性，目前在鋼筋混凝土結(jié)構(gòu)低速沖擊領(lǐng)域應(yīng)用廣泛。

圖10 不同荷載下混凝土應(yīng)變率

核島廠房使用C40混凝土材料，根據(jù)GB 50010—2010《混凝土結(jié)構(gòu)設(shè)計規(guī)范》，抗壓強(qiáng)度為27 MPa，密度為2 400 kg/m3，這兩個參數(shù)為混凝土本構(gòu)所需輸入?yún)?shù)。核島廠房中鋼板厚度為20 mm，環(huán)向鋼筋直徑為40 mm，縱向拉筋直徑為20 mm，密度均為7 800 kg/m3,泊松比為0.3，使用可以考慮應(yīng)變率與失效的Plastic-Kinematic(塑性隨動硬化)本構(gòu)模型，材料參數(shù)如表1所示。

表1 核島廠房鋼筋、鋼板材料參數(shù)

1.2.2 核島廠房網(wǎng)格敏感性分析

在數(shù)值模擬之前，需要對有限元模型進(jìn)行網(wǎng)格敏感性分析，以保證計算的效率和準(zhǔn)確性。為了選出最合適的網(wǎng)格尺寸，以網(wǎng)格尺寸和反應(yīng)堆廠房沿徑向劃分網(wǎng)格個數(shù)為標(biāo)準(zhǔn)，構(gòu)建5種不同類型的模型。各個模型的節(jié)點、單元個數(shù)等信息，如表2所示。

表2 不同類型廠房模型網(wǎng)格信息

在進(jìn)行分析時，使用1.1.2節(jié)分析中選擇的飛機(jī)模型1，飛機(jī)以100 m/s速度垂直撞擊反應(yīng)堆廠房墻體，撞擊點坐標(biāo)為：(52 015.0,61 116.9,16 434.0)，如圖11所示。通過對比這5種模型的最大變形位移，進(jìn)而選擇合適的網(wǎng)格尺寸。

圖11 核島廠房網(wǎng)格敏感性分析飛機(jī)撞擊位置圖

將計算得到的不同模型反應(yīng)堆廠房墻體在飛機(jī)撞擊過程中，最大變形位移結(jié)果如表3所示。從表3數(shù)據(jù)可知，300_R006模型和400_R006模型最大變形值相差不大，而500_R006比這兩個模型增加了10 mm多。反應(yīng)堆廠房模型最大變形值沿著徑向6層單元和5層差不多，4層則比這兩個模型增加了50 mm多。從計算效率來看，300_R006的網(wǎng)格數(shù)量超過400_R006約20萬，如果采用300 mm網(wǎng)格尺寸會大大增加計算耗時，而6層比5層則超過5萬多，對計算時間影響不是很大，最終選擇了400_R006模型，整體網(wǎng)格尺寸400 mm，反應(yīng)堆廠房外墻沿著徑向分布有6層單元。

表3 不同模型最大變形位移

2 考慮土-結(jié)構(gòu)相互作用效應(yīng)的飛機(jī)撞擊分析

在工程應(yīng)用中，通常將地基假定為剛性，在核島廠房底板底部施加固端約束。這種假定在柔性地基場地條件下，與實際情況存在較大差異。因此，本文考慮使用直接法，對完整的地基場地進(jìn)行建模，并在地基外邊界設(shè)置虛擬人工邊界，以模擬無限地基輻射阻尼效應(yīng)。

2.1 結(jié)構(gòu)地基相互作用模型

該相互作用模型是在核島廠房模型的基礎(chǔ)上，增加了地基模型而構(gòu)建的。地基模型長300 m，寬250 m，深45 m，如圖12所示。

圖12 地基模型(m)

由于核島廠房外輪廓并非規(guī)則幾何圖形，因此在網(wǎng)格剖分時，靠近核島廠房局部范圍內(nèi)，根據(jù)核島廠房外輪廓進(jìn)行合理剖分，其余整體網(wǎng)格尺寸為2 m。整個地基模型共包含約433萬個節(jié)點和493萬個單元。在地基外側(cè)5個立面上，設(shè)置虛擬人工邊界——完美匹配層(perfect matching layer，PML)，以考慮地基無限半空間的影響。

考慮到場地類型的多樣性，選取了AP1000 3種標(biāo)準(zhǔn)設(shè)計的場地類型[17]以及一種滿足剪切波速包絡(luò)性的軟土場地，進(jìn)行算例對比分析。地基土采用理想彈塑性Drucker-Prager本構(gòu)模型，該模型考慮了土的非線性。不同場地土體參數(shù)參考GB 50218—1994《工程巖體分級標(biāo)準(zhǔn)》進(jìn)行設(shè)置，如表4所示。

表4 不同場地土體參數(shù)

同時，還考慮了3種重要撞擊位置對結(jié)構(gòu)響應(yīng)的影響，分別是飛機(jī)水平撞擊核島廠房外殼(撞擊位置1,)，飛機(jī)水平撞擊空冷塔(撞擊位置2)和飛機(jī)斜向垂直撞擊穹頂(撞擊位置3)，如圖13、圖14、圖15所示。

圖13 水平撞擊核島廠房

圖14 水平撞擊空冷塔

圖15 斜向垂直撞擊穹頂

2.2 虛擬人工邊界

PML是在計算區(qū)域外人為添加的一具有一定厚度、含任意吸收介質(zhì)的吸收層。為驗證PML的精度和可靠性，對圖16所示具有解析解的半無限介質(zhì)波動問題，分別通過解析和數(shù)值方法進(jìn)行計算。

圖16 內(nèi)源問題

地基模型在構(gòu)建時，由于基礎(chǔ)底板形狀不規(guī)則，在網(wǎng)格剖分時，難免會出現(xiàn)三棱柱網(wǎng)格。為分析這樣的網(wǎng)格是否對計算結(jié)果有影響，構(gòu)建了如圖17所示兩種模型。

圖17 不同網(wǎng)格PML模型

模型1包含三棱柱單元，模型2全是六面體單元。該內(nèi)源問題的對比結(jié)果，如圖18所示。從圖18可知，模型1和模型2結(jié)果非常吻合，可見三棱柱單元的影響可忽略。雖然PML和黏彈邊界的結(jié)果同解析解有些差別，但這兩種方法的計算結(jié)果相差不大，PML的計算結(jié)果精度同黏彈邊界的計算精度相當(dāng)。

圖18 結(jié)果對比

2.3 結(jié)果分析

2.3.1 不同場地條件的影響

模型中設(shè)置了兩個監(jiān)測點對計算結(jié)果進(jìn)行分析，點A設(shè)置在反應(yīng)堆廠房底板中心處，點B設(shè)置在核島廠房某個艙室的內(nèi)壁上，如圖19所示。

圖19 監(jiān)測點

在分析不同場地對上部結(jié)構(gòu)振動響應(yīng)的影響時，飛機(jī)以100 m/s速度撞擊。飛機(jī)撞擊位置1時，點A的三向加速度反應(yīng)譜結(jié)果如圖20所示，點B的三向加速度反應(yīng)譜結(jié)果如圖21所示。

圖20 撞擊位置1點A的加速度反應(yīng)譜

圖21 撞擊位置1點B的加速度反應(yīng)譜

從不同場地的對比結(jié)果可以看出，在約4 Hz頻率處，加速度反應(yīng)譜存在第一個峰值，該頻率與結(jié)構(gòu)基頻相當(dāng)，在此頻率下結(jié)構(gòu)共振響應(yīng)較大。

加速度峰值作為衡量振動響應(yīng)程度的重要指標(biāo)，其值可通過加速度反應(yīng)譜在極高頻處的譜值進(jìn)行估計。由于Y向是飛機(jī)撞擊方向，點A在不同工況下的Y向加速度反應(yīng)譜峰值及1 000 Hz處的加速度峰值，如表5所示。

表5 撞擊位置1點A處撞擊正向Y動力響應(yīng)

從監(jiān)測點的對比結(jié)果可以看出，場地對反應(yīng)堆廠房振動響應(yīng)有明顯影響。工況4場地的加速度反應(yīng)譜普遍高于其他，此時場地條件最為柔軟，因此我們可以判斷場地越軟，加速度反應(yīng)譜將會越大，此時，土結(jié)構(gòu)相互作用效應(yīng)較為顯著。同時，加速度反應(yīng)譜在極高頻處的譜值等于加速度響應(yīng)的峰值，從工況1～工況4逐漸增大反映出，場地條件越柔，結(jié)構(gòu)的振動響應(yīng)越劇烈。因此，在核電廠房抗飛機(jī)撞擊引起的核電結(jié)構(gòu)整體動力響應(yīng)中，需要考慮土-結(jié)構(gòu)相互作用。

2.3.2 撞擊位置影響

飛機(jī)撞擊不同位置時，點A的三向加速度反應(yīng)譜結(jié)果對比如圖22所示，點B的三向加速度反應(yīng)譜結(jié)果對比如圖23所示。

圖22 不同撞擊位置點A的加速度反應(yīng)譜

圖23 不同撞擊位置點B的加速度反應(yīng)譜

由于Y向是飛機(jī)撞擊方向，點A在不同撞擊工況下的Y向加速度峰值，如表6所示。

表6 不同撞擊位置點A處撞擊正向Y動力響應(yīng)

由于撞擊位置3的高度最高，從以上結(jié)果可以看出，撞擊位置高時，結(jié)構(gòu)的振動響應(yīng)增大。因此，飛機(jī)撞擊位置的增高，亦會較大幅度地加大結(jié)構(gòu)的沖擊響應(yīng)。而由于空冷塔的存在，部分削弱了飛機(jī)作用在反應(yīng)堆廠房上的撞擊能量，導(dǎo)致撞擊位置2在多數(shù)情況下高頻區(qū)段加速度反應(yīng)譜值小于撞擊位置1。

3 結(jié) 論

通過以上分析可以看出，土-結(jié)構(gòu)相互作用的影響在核電站抗飛機(jī)撞擊問題中需要被考慮。

(1) 地基剪切波速由1 066.8 m/s降低為200 m/s時，在飛機(jī)撞擊水平Y(jié)向上，結(jié)構(gòu)基底加速度時程響應(yīng)峰值增加了約93%。因此，軟土地基的土-結(jié)構(gòu)相互作用效應(yīng)十分顯著。

(2) 在工況1場地條件下，飛機(jī)撞擊水平Y(jié)向，撞擊位置由核島廠房升高到穹頂，加速度反應(yīng)譜峰值提高了37.5%，結(jié)構(gòu)基底加速度時程響應(yīng)峰值增加了138.2%，說明撞擊位置的增高，加大了結(jié)構(gòu)的沖擊響應(yīng)。

綜上，場地條件與撞擊位置，在核電站抗飛機(jī)撞擊引起的土-結(jié)構(gòu)相互作用中是需要重點分析的兩項因素。