爆炸荷載作用下大型箱涵安全性的SPH-FEM耦合分析

崔 溦，盧珊珊，宋慧芳，張社榮

(天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

地下結(jié)構(gòu)物的安全性是工程界關(guān)注的重要問題．以南水北調(diào)中線工程的大型輸水箱涵為例，其單孔凈尺寸達(dá)到 4.4,m×4.4,m，一旦受到損壞，將直接影響到沿線20多座大中城市的用水安全．

爆炸沖擊荷載作用下，地下結(jié)構(gòu)物的響應(yīng)與地面結(jié)構(gòu)物有很大不同，其最主要的差別在于必須考慮土與結(jié)構(gòu)物的相互作用．雖然修正后的單自由度分析方法可以近似考慮這種相互作用，但在模型建立、參數(shù)選擇以至計(jì)算結(jié)果上都存在較大復(fù)雜性和不確定性[1]．?dāng)?shù)值分析方法可以精確模擬爆炸沖擊下地下結(jié)構(gòu)物的響應(yīng)，但現(xiàn)有的研究方法，無論是有限元、有限差分還是它們之間的耦合算法，一般將爆破過程與結(jié)構(gòu)物響應(yīng)割裂開來，直接將爆破荷載以壓力的形式施加到結(jié)構(gòu)物上，其計(jì)算精度難以得到保證[2-3]．精確模擬地下結(jié)構(gòu)物爆破響應(yīng)的主要難點(diǎn)在于數(shù)值模型必須包含藥包周圍的大變形土體與地下結(jié)構(gòu)物的復(fù)雜體型，且要在計(jì)算精度與計(jì)算效率上都能滿足工程要求[4]．

筆者以 SPH法模擬藥包周圍土體，以拉格朗日有限元法模擬復(fù)雜結(jié)構(gòu)物，通過二者之間耦合研究大型輸水箱涵在爆破荷載作用下的響應(yīng)，以為工程設(shè)計(jì)與維護(hù)提供一定參考．

1 SPH與FE耦合算法

在采用常規(guī)有限元方法分析爆炸引起的大變形問題時(shí)，由于網(wǎng)格的劇烈扭曲經(jīng)常會(huì)導(dǎo)致計(jì)算困難，為了克服這種困難，光滑粒子流體動(dòng)力學(xué)(SPH)作為一種重要的無網(wǎng)格法得到了較為廣泛的應(yīng)用[5]．

SPH法把連續(xù)的物理量用多數(shù)粒子的集合來進(jìn)行插值求解，連續(xù)體的運(yùn)動(dòng)用有限數(shù)量的粒子運(yùn)動(dòng)來離散化，由于沒有使用網(wǎng)格，因此不會(huì)發(fā)生界面變形大所引起的計(jì)算溢出問題．計(jì)算空間導(dǎo)數(shù)時(shí)不需要使用任何網(wǎng)格，而是通過一個(gè)稱為“核函數(shù)”的積分核進(jìn)行“核函數(shù)估值”近似，將流體動(dòng)力學(xué)基本方程組轉(zhuǎn)換成數(shù)值計(jì)算的 SPH方程，整個(gè)流場(chǎng)離散成一系列“粒子”，所有力學(xué)量(密度、速度、壓力、內(nèi)能等)由這些粒子負(fù)載，這些粒子可以按計(jì)算公式，亦即按流體動(dòng)力學(xué)運(yùn)動(dòng)的規(guī)律任意地運(yùn)動(dòng)．在 SPH法中，每個(gè)粒子 I與在和它相距給定距離范圍內(nèi)(一般假定為2,h，h為光滑長(zhǎng)度)的其他所有粒子J發(fā)生相互作用．SPH法要引入一個(gè)特殊函數(shù)W(x-x′,h)作為核函數(shù)，場(chǎng)函數(shù)經(jīng)過核函數(shù)“光滑化”再在整個(gè)求解域上積分，便得到了場(chǎng)函數(shù)的核估計(jì)．對(duì)場(chǎng)函數(shù) f，其核估計(jì)表示為

式中x和 x′均為位置矢量．

圖 1給出了核估計(jì)的基本概念．細(xì)節(jié)問題可參見文獻(xiàn)[6]．

SPH法在模擬復(fù)雜體型的類似于墻的結(jié)構(gòu)物時(shí)，由于相對(duì)于結(jié)構(gòu)尺寸時(shí)墻體厚度較小，也就需要較小顆粒和時(shí)間步長(zhǎng)，在計(jì)算效率與精度上，該方法受到一定限制．采用SPH與FE耦合技術(shù)，在材料小變形和結(jié)構(gòu)復(fù)雜區(qū)域采用有限元模擬，在爆炸近域采用SPH模擬，可以克服上述缺點(diǎn)，并保證計(jì)算精度．

SPH顆粒和FEM網(wǎng)格有2種耦合方法，一種是通過可以相對(duì)滑動(dòng)的接觸面連接，另一種則是固接．在本次分析中，由于SPH顆粒和FEM網(wǎng)格的交界面位于土介質(zhì)內(nèi)，顆粒與網(wǎng)格之間不存在相對(duì)變形而固接在一起，因此根據(jù)運(yùn)動(dòng)方程，其他顆粒傳遞到交界面顆粒上的力與有限元網(wǎng)格傳遞到交界面顆粒上的力相同．圖 2給出了 SPH粒子與傳統(tǒng)拉格朗日有限元網(wǎng)格之間的固結(jié)耦合關(guān)系．

圖1 SPH法的核估計(jì)Fig.1 SPH approximation

圖2 SPH與FEM的耦合關(guān)系Fig.2 Coupling of SPH and FEM

2 材料模型

2.1 TNT

沖擊波是爆破過程的重要特征．在本次研究中，炸藥狀態(tài)方程用JWL(Jones-Wilkins-Lee)方程來描述．JWL狀態(tài)方程可表示為

式中：v為比容，1/vρ=；e為比內(nèi)能；1C、1r、2C、2r和ω為常數(shù)，且對(duì)于常規(guī)爆破來說，它們的數(shù)值可以由動(dòng)力試驗(yàn)來確定．本次數(shù)值分析所用 TNT的材料參數(shù)如表1所示．

表1 炸藥的材料參數(shù)Tab.1 Material parameters for TNT

2.2 土

土用沖擊狀態(tài)方程和基于 D-P準(zhǔn)則的彈塑性模型來描述，并定義其拉伸極限．

對(duì)于沖擊荷載，即使在沖擊速度約 2倍于初始聲速0c和沖擊壓力達(dá)到 100,GPa量級(jí)的情況下，Rankine-Hugoniot狀態(tài)方程都可以表示為

采用 D-P準(zhǔn)則和拉伸極限來描述土的強(qiáng)度特性．為了避免土的剪脹問題，采用非關(guān)聯(lián)的流動(dòng)法則，其塑性勢(shì)函數(shù)定義為

式中Y為屈服極限．

本次數(shù)值分析采用土體的具體參數(shù)如表2所示．

表2 土的材料參數(shù)Tab.2 Material parameters for soil

2.3 混凝土

沖擊荷載下混凝土響應(yīng)是一個(gè)復(fù)雜的非線性和率相關(guān)過程，在目前分析中，以 RHT模型應(yīng)用最多，RHT模型可以分為強(qiáng)度模型和狀態(tài)方程模型2部分[7]．對(duì)于壓碎材料來說，強(qiáng)度模型主要包括彈性極限函數(shù)、破壞函數(shù)和殘余強(qiáng)度函數(shù)3部分．

在硬化階段之后，混凝土塑性應(yīng)變導(dǎo)致的材料損傷和強(qiáng)度降低表示為

在該模型中，狀態(tài)方程采用 P-α模型描述，其主要假定是在相同的壓力和溫度條件下，多孔介質(zhì)材料的比內(nèi)能和相同密度實(shí)心材料是一致的，該模型重點(diǎn)關(guān)注高應(yīng)力下的混凝土性態(tài)，但同時(shí)也可以提供低應(yīng)力狀態(tài)下壓縮過程的合理描述．混凝土模型的主要參數(shù)如表3～表5所示.

表3 混凝土的RHT模型物理參數(shù)Tab.3 Physical parameters for RHT model of concrete

表4 混凝土的RHT模型強(qiáng)度參數(shù)Tab.4 Strength ntensity parameters for RHT model of concrete

表5 混凝土的RHT模型狀態(tài)參數(shù)Tab.5 State parameters for RHT model of concrete

2.4 鋼 筋

在本次分析中，鋼筋采用率相關(guān)的John-Cook彈塑性模型模擬．屈服應(yīng)力表示為

在本次分析中，直接將鋼筋和混凝土單元共結(jié)點(diǎn)連接，鋼筋材料參數(shù)如表6所示．

表6 鋼筋材料參數(shù)Tab.6 Parameters used for reinforcement steel bar

2.5 接觸面模型

根據(jù) Huck的試驗(yàn)成果，當(dāng)土與結(jié)構(gòu)物的接觸面比較粗糙時(shí)，接觸面應(yīng)力超出土的最大剪應(yīng)力時(shí)發(fā)生破壞，且高壓力時(shí)接觸面的強(qiáng)度特性與土體屬性非常接近[8]．基于上述原因，本次分析將土與混凝土箱涵之間接觸面按完全連接考慮．

2.6 邊界條件

為滿足輻射條件，數(shù)值模型中土為透射邊界．透射邊界條件允許應(yīng)力波通過物理邊界傳輸而沒有反射[5]．假定邊界的法向速度矢量為 un，邊界處壓力 p采用計(jì)算式如下：

式中：pref和 uref分別為相關(guān)壓力和速度分量；I為材料阻抗．

3 數(shù)值模型

圖3給出本次分析的SPH-FEM耦合模型，考慮無水工況，模型共包含鋼筋混凝土箱涵、土體、炸藥3種介質(zhì)，并近似按二維軸對(duì)稱平面應(yīng)變問題分析．炸藥及周邊土體采用 SPH模擬，箱涵及周邊土體采用FEM模擬，為了提高計(jì)算效率，近炸藥及箱涵周邊土體 SPH和 FEM 網(wǎng)格較密，距離較遠(yuǎn)則網(wǎng)格逐漸稀疏，F(xiàn)EM 網(wǎng)格大約 8萬個(gè)，采用二維四結(jié)點(diǎn)平面單元，高斯積分法；SPH粒子大約2萬個(gè)．

圖3 SPH-FEM耦合模型Fig.3 SPH-FEM coupled model

圖 4給出本次分析模型的具體尺寸及計(jì)算中監(jiān)測(cè)點(diǎn)的位置．TNT質(zhì)量為 80,kg，埋深為 d(1,m，2,m)，距箱涵最小水平距離為 h(6,m，8,m)．箱涵單孔凈尺寸為 4.4,m×4.4,m，鋼筋布置及其他尺寸如圖4所示，鋼筋保護(hù)層厚度為50,mm．

圖4 模型布置Fig.4 Configuration of model

4 結(jié)果分析

4.1 爆坑形成

圖 5給出不同埋藥深度下土中爆坑形成過程．從圖 5中可以看出，隨著埋藥深度的增加，爆坑形態(tài)發(fā)生明顯變化，并逐漸由可見爆坑向爆腔過渡．?dāng)?shù)值模擬很好地揭示了爆破鼓包的運(yùn)動(dòng)過程和表面土體的層狀剝離現(xiàn)象．在埋藥深度1,m情況下，以拋擲爆破為主，形成爆坑最大直徑約為 3.8,m．在埋藥深度2,m情況下，以松動(dòng)爆破形式為主．在地面形成比較明顯的松動(dòng)破壞區(qū)，主要為自由面的反射拉伸波所引起，最大空腔直徑約為 4.2,m．根據(jù)Henrych[9]針對(duì)不同屬性土、不同爆破方式下的爆腔經(jīng)驗(yàn)公式計(jì)算得到爆腔直徑分別為 3.67,m和4.18,m，與數(shù)值模擬結(jié)果吻合較好．

圖5 不同埋藥深度下爆坑形成過程Fig.5 Crater formation under different buried depths

4.2 土中爆炸沖擊波傳播

圖 6給出了爆破后土中速度時(shí)程曲線．圖中測(cè)點(diǎn)位置如圖 4所示，其中測(cè)點(diǎn) 8、9位于 SPH網(wǎng)格內(nèi)，測(cè)點(diǎn) 10、11位于 FEM 網(wǎng)格內(nèi)．從圖 6中可以看出，隨著與藥包距離增加，速度峰值逐漸減?。凰幇虻膲毫ψ兓哂袥_擊波的特征，單峰值且持續(xù)時(shí)間很短，隨著與藥包距離增加，速度變化開始具有應(yīng)力波的特征，幅值減小且持續(xù)時(shí)間增加．通過與文獻(xiàn)[9-10]對(duì)比可以看出，無論是測(cè)點(diǎn) 8、9還是測(cè)點(diǎn) 10、11，速度時(shí)程曲線特征與試驗(yàn)和理論計(jì)算結(jié)果一致，表明本文所采取的分析方法是較準(zhǔn)確的．

圖6 土中沖擊速度時(shí)程Fig.6 Velocity time history in soil

4.3 混凝土損傷

圖7給出不同爆破距離下箱涵損傷分布．從圖7中可以看出，損傷嚴(yán)重區(qū)域主要出現(xiàn)在邊墻與底板和頂板交界處，主要為沖擊荷載引起的壓損傷．爆炸距離對(duì)箱涵累積損傷有較大影響，隨著水平距離h和埋深d的增加，累積損傷區(qū)域逐漸減小．

圖7 不同爆破距離下混凝土損傷分布Fig.7 Distribution of damage in concrete at different detonation distances

4.4 箱涵響應(yīng)

圖8給出不同爆炸距離下箱涵變形情況．

圖8 混凝土變形時(shí)程曲線Fig.8 Concrete deformation time history curves

從圖 8中可以看出，不同測(cè)點(diǎn)變形時(shí)程差別明顯．從測(cè)點(diǎn)1、2、3的水平向變形來看，不同爆炸距離下變形規(guī)律基本一致，邊墻中間位置水平變形最大，最大量達(dá)到 12.5,mm，出現(xiàn)在最近爆炸距離(D＝1,m、h＝6,m)情況，殘余變形量為不可回復(fù)的塑性變形，而邊墻與頂板底板交界處變形較?。畯臏y(cè)點(diǎn) 1、4、6、7的豎向變形來看，不同爆炸距離下變形規(guī)律存在一定差別，爆炸深度越大，箱涵整體將呈現(xiàn)一定向上移動(dòng)趨勢(shì)．

圖9給出不同爆炸距離下箱涵應(yīng)力變化情況．從圖 9中可以看出，各測(cè)點(diǎn)應(yīng)力時(shí)程差別明顯．從測(cè)點(diǎn)

1、2、4的水平向應(yīng)力來看，變化規(guī)律基本一致．在沖擊峰值之后，存在較長(zhǎng)時(shí)間的擺動(dòng)現(xiàn)象，主要為結(jié)構(gòu)自身響應(yīng)引起．水平向應(yīng)力最大值為-33,MPa，出現(xiàn)在最近爆炸距離(D＝1,m、h＝6,m)情況．邊墻中間位置的應(yīng)力時(shí)程與邊墻和頂板交界處的應(yīng)力時(shí)程明顯不同，存在明顯的二次加載現(xiàn)象，主要為箱涵周邊土體流動(dòng)變形導(dǎo)致壓應(yīng)力增加引起[11]．從測(cè)點(diǎn) 1、2、4的豎向應(yīng)力來看，變化規(guī)律存在較大差別．測(cè)點(diǎn) 1存在較為明顯的沖擊峰值，2、4點(diǎn)則不明顯．不同爆炸距離下測(cè)點(diǎn)2的應(yīng)力時(shí)程差別較大，拉壓應(yīng)力的變化與埋深密切相關(guān)．

圖9 混凝土應(yīng)力時(shí)程曲線Fig.9 Concrete stress time histories curves

5 結(jié) 論

(1) 采用 SPH 法模擬爆破近域土體大變形，F(xiàn)E法模擬遠(yuǎn)域土體及箱涵結(jié)構(gòu)，爆破后40,ms計(jì)算時(shí)間約為 14,h，具有較好的計(jì)算效率，計(jì)算模型較好地揭示了土中爆坑形成和沖擊波傳播過程，計(jì)算精度滿足要求．

(2) 不同爆炸位置下箱涵結(jié)構(gòu)爆破響應(yīng)存在較大差異，尤其是邊墻中間位置，受土體變形流動(dòng)影響，存在較為明顯二次加載現(xiàn)象，且爆破位置對(duì)其拉壓變化影響明顯．

［1］ Wdidlinger P，Hinman E. Analysis of underground protective structures[J].J Struct Eng，1987，114(7)：1658-1673.

［2］ Lu Yong，Wang Zhongqi，Chong Karen. A comparative study of buried structures in soil subjected to blast load using 2D and 3D numerical simulations[J].Soil Dynamics and Earthquake Engineering，2005，25(4)：275-288.

［3］ Wang Zhongqi，Lu Yong，Hao Hong. A full coupled numerical analysis approach for buried structures subjected to subsurface blast[J].Computersand Structures，2005，83(4/5)：339-356.

［4］ 杜義欣，劉晶波，伍 俊，等. 常規(guī)爆炸下地下結(jié)構(gòu)的沖擊震動(dòng)環(huán)境[J]. 清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2006，46(3)：322-326.

Du Yixin，Liu Jingbo，Wu Jun，et al. Blast shock and vibration of underground structures with conventional weapon[J].Journal of Tsinghua University：Science and Technology，2006，46(3)：322-326(in Chinese).

［5］ 王 吉，王肖鈞，卞 梁. 光滑粒子法與有限元的耦合算法及其在沖擊動(dòng)力學(xué)中的應(yīng)用[J]. 爆炸與沖擊，2007，27(6)：522-528.

Wang Ji，Wang Xiaojun，Bian Liang. Linking of smoothed particle hydrodynamic method to standard finite element method and its application in impact dynamics [J].Explosion and Shock Waves，2007，27(6)：522-528(in Chinese).

［6］ Liu G R，Liu M B.Smoothed Particle Hydrodynamics：A Meshfree Particle Method[M]. London：World Scientific，2003.

［7］ 王 政，倪玉山，曹菊珍. 沖擊載荷下混凝土動(dòng)態(tài)力學(xué)性能研究進(jìn)展[J]. 爆炸與沖擊，2005，25(6)：519-524.

Wang Zheng，Ni Yushan，Cao Juzhen. Recent advances of dynamic mechanical behavior of concrete under impact loading[J].Explosion and Shock Waves，2005，25(6)：519-524(in Chinese).

［8］ Huck P J，Saxena S K. Response of soil-concrete interface at high pressure[C]//Proceedings of the Tenth International Conference on Soil Mechanics and Foundation Engineering. Stockholm，Sweden，1981：141-144.

［9］ Henrych J.The Dynamics of Explosion and Its Use[M].Amsterdam：Elsevier Scientific Publishing Company，1979.

［10］ Wang Z，Hao H，Lu Y. A three-phase soil model for simulating stress wave propagation due to blast loading[J].Int J Numer Anal Geomech，2004，28(1)：33-56.

［11］ James T，Baylot P E. Effect of soil flow changes on structures loads[J].J Struct Eng，2000，126(12)：1434-1441.