土中爆炸作用下箱涵動力響應(yīng)的SPH－FE耦合分析＊

崔 溦，宋慧芳，張社榮

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

爆炸荷載作用下大型箱涵的安全性是工程界關(guān)注的重要問題，主要涉及兩個方面，一是沖擊荷載在土中的傳播，二是土與結(jié)構(gòu)物的動力相互作用。在理論計(jì)算方面，雖然修正后的單自由度分析方法可以近似考慮這種相互作用，但在荷載輸入、模型建立、參數(shù)選擇以至計(jì)算結(jié)果上都存在較大復(fù)雜性和不確定性［1］。數(shù)值分析方法可以精確模擬爆炸沖擊下地下結(jié)構(gòu)物的響應(yīng)，但現(xiàn)有的研究方法，無論是有限元、有限差分還是他們之間的耦合算法，一般將爆破過程與結(jié)構(gòu)物響應(yīng)割裂，直接將爆破荷載以壓力的形式施加到結(jié)構(gòu)物上，計(jì)算結(jié)果的科學(xué)性難以得到保證［2－3］。同時(shí)，由于爆炸引起的土體大變形的影響，已有的數(shù)值算法會因網(wǎng)格扭曲而導(dǎo)致模擬整個爆炸過程存在困難［4］。精確模擬地下結(jié)構(gòu)物爆破響應(yīng)的主要難點(diǎn)，在于數(shù)值模型必須包含藥包周圍的大變形土體與鋼筋混凝土箱涵的復(fù)雜體型，并且在計(jì)算結(jié)果與計(jì)算效率上都能滿足工程要求［5］。

本文中，以SPH法模擬藥包周圍土體，以拉格朗日有限元法模擬鋼筋混凝土箱涵，通過二者之間耦合研究大型箱涵在爆炸荷載作用下的響應(yīng)，擬為工程設(shè)計(jì)與維護(hù)提供一定參考。

1 SPH與FE耦合算法

在采用常規(guī)有限元方法分析爆炸引起的大變形問題時(shí)，由于網(wǎng)格的劇烈扭曲經(jīng)常會導(dǎo)致計(jì)算困難。為了克服這種困難，光滑粒子流體動力學(xué)（SPH）作為一種重要的無網(wǎng)格法得到了較廣泛的應(yīng)用［6］。但SPH法在模擬復(fù)雜體型類似于墻的結(jié)構(gòu)物時(shí)，由于相對于結(jié)構(gòu)尺寸時(shí)墻體厚度較小，需要較小顆粒和時(shí)間步長，所以在計(jì)算效率與精度上，該方法受到一定限制。采用SPH與FE耦合技術(shù)，在材料小變形和結(jié)構(gòu)復(fù)雜區(qū)域采用有限元模擬，在爆炸近域采用SPH模擬，可以克服網(wǎng)格的劇烈扭曲而無法計(jì)算的困難。

SPH顆粒和FEM網(wǎng)格有兩種耦合方法，一種是通過可以相對滑動的接觸面連接，另一種則是固接。在本文中，由于SPH顆粒和FEM網(wǎng)格的交界面位于土介質(zhì)內(nèi)，顆粒與網(wǎng)格之間不存在相對變形而固接在一起。因此根據(jù)運(yùn)動方程，其他顆粒傳遞到交界面顆粒上的力與有限元網(wǎng)格傳遞到交界面顆粒上的力相同。

2 材料模型

2.1 炸 藥

炸藥產(chǎn)物狀態(tài)方程用JWL方程描述

式中：v＝1／ρ，為比容，A、B、R1、R2和ω為常數(shù)，對于常規(guī)爆破，他們可以由動力實(shí)驗(yàn)確定。

TNT材料參數(shù)分別為：ρ＝1.63g／cm3，D＝6.93km／s，A＝374GPa，B＝3.73GPa，R1＝4.15，R2＝0.9，ω＝0.35，e＝6GJ／m3，pCJ＝21GPa。

2.2 土

土用沖擊狀態(tài)方程和基于D－P準(zhǔn)則的彈塑性模型描述，并定義拉伸極限。為了避開土的剪脹問題，采用非關(guān)聯(lián)的流動法則。在沖擊荷載下，即使沖擊速度為初始聲速c0的約兩倍、沖擊壓力達(dá)到百吉帕量級，Rankine－ Hugoniot狀態(tài)方程都可以采用下式

式中：p為靜水壓力，ρ0為初始密度，e為比內(nèi)能，Γ為Grüneisen參數(shù)。

土體參數(shù)為：ρ0＝1.92g／cm3，G＝220MPa，ft＝－100kPa，Γ＝0.11，c0＝1.614km／s，s＝1.5。

2.3 混凝土

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

在硬化階段后，混凝土附加塑性應(yīng)變導(dǎo)致的材料損傷和強(qiáng)度降低采用下式

式中：D1、D2為損傷因子，εf，min為混凝土達(dá)到破壞最小應(yīng)變，εp為塑性應(yīng)變，p＊＝p／fc，＝p＊（ft／fc），ft、fc分別為混凝土的抗拉強(qiáng)度和抗壓強(qiáng)度。

在該模型中，狀態(tài)方程采用p－α模型描述

式中：A1、A2、A3、T1、T2為常量，v0為多孔介質(zhì)材料的初始比容；αini為初始孔隙率，pe為多孔介質(zhì)中的骨架開始破壞時(shí)的壓力，ps為多孔介質(zhì)完全壓縮時(shí)的壓力，n為常數(shù)。

混凝土模型的主要參數(shù)分別為：ρ0＝2.3t／m3，ρs＝2.75t／m3，cV＝654J／（kg·K）；A＝1.6，N＝0.61，n1＝0.036，n2＝0.032，εf，min＝0.01，D1＝0.04，D2＝1.0，ft＝3.5MPa，fc＝35MPa；A1＝35GPa，A2＝39GPa，A3＝9GPa，T1＝35GPa，T2＝0，n＝3.0。

2.4 鋼 筋

鋼筋采用率相關(guān)的John－Cook彈塑性模型。屈服應(yīng)力為

式中：Y0為初始屈服強(qiáng)度，εp為有效塑性應(yīng)變?yōu)闅w一化的有效塑性應(yīng)變，B、C、n、m為材料常數(shù)。TH為相應(yīng)溫度，TH＝（T－Tr）／（Tm－Tr），Tm為熔化溫度，Tr為環(huán)境溫度。式（8）中的第一部分表示應(yīng)變硬化效應(yīng)，第二部分和第三部分分別表示應(yīng)變率和溫度效應(yīng)。

鋼筋材料參數(shù)分別為：ρ＝7.9g／cm3，E＝200GPa，cV＝450J／（kg·K），G＝82GPa，Y0＝350MPa，B＝275MPa，C＝0.022，n＝0.36，m＝1.0，Tr＝300K，Tm＝1 800K。

圖1 SPH－FEM耦合模型Fig.1 SPH－FEM coupled model

2.5 接觸面模型

根據(jù)實(shí)驗(yàn)結(jié)果［8］，當(dāng)土與結(jié)構(gòu)物的接觸面比較粗糙時(shí)，接觸面應(yīng)力超出土的最大剪應(yīng)力時(shí)發(fā)生破壞，并且高壓力時(shí)接觸面的強(qiáng)度特性與土體屬性非常接近。因此，將土與混凝土箱涵間接觸面按完全連接考慮。

2.6 邊界條件

為了滿足輻射條件，數(shù)值模型中土為透射邊界。透射邊界條件允許應(yīng)力波通過物理邊界傳輸而沒有反射［5］。

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

3 數(shù)值模型和結(jié)果分析

圖1給出了SPH－FEM耦合模型，模型包含鋼筋混凝土箱涵、土體、炸藥等三種介質(zhì)，并近似為二維軸對稱平面應(yīng)變問題。炸藥及周邊土體采用SPH模擬，箱涵及周邊土體采用FEM模擬，為了提高計(jì)算效率，近炸藥及箱涵周邊土體SPH和FEM網(wǎng)格較密，距離較遠(yuǎn)則網(wǎng)格逐漸稀疏，F(xiàn)EM網(wǎng)格約8萬個，SPH粒子約2萬個。

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

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

3.1 混凝土損傷

圖3給出了不同布藥位置下箱涵損傷分布。從圖中可以看出，損傷嚴(yán)重區(qū)域主要出現(xiàn)在邊墻和邊墻與底板、頂板交界處，爆炸距離對箱涵累積損傷有較大影響，隨著水平距離h和埋深d的增加，累積損傷區(qū)域逐漸減小。

3.2 箱涵響應(yīng)

不同爆炸距離下箱涵變形如圖4所示。可以看出，不同測點(diǎn)變形時(shí)程差別明顯。由測點(diǎn)1、2的水平變形可見，不同爆炸距離下變形趨勢基本一致，邊墻中間位置水平變形最大，達(dá)12.5mm，出現(xiàn)在最近爆炸距離情況，殘余變形為不可回復(fù)的塑性變形，而邊墻與頂板交界處變形較小。由測點(diǎn)1、3的豎向變形可見，不同爆炸距離下變形規(guī)律基本一致。在爆炸作用下，箱涵整體呈向爆炸位置傾斜的趨勢。

圖4 （a）箱涵的水平變形Fig.4 （a）Horizontal deformations of box culvert

圖4 （b）箱涵的豎向變形Fig.4 （b）Vertical deformations of box culvert

圖5 （a）箱涵測點(diǎn)1的應(yīng)力Fig.5 （a）Stresses at point 1of box culvert

圖5 （b）箱涵測點(diǎn)2的應(yīng)力Fig.5 （b）Stresses at point 2of box culvert

圖5 （c）箱涵測點(diǎn)3的應(yīng)力Fig.5 （c）Stresses at point 3of box culvert

不同爆炸距離下箱涵應(yīng)力變化如圖5所示。可以看出，各測點(diǎn)應(yīng)力時(shí)程差別明顯。由測點(diǎn)1、3的水平應(yīng)力可見，變化規(guī)律基本一致。在沖擊峰值后，存在較長時(shí)間的擺動現(xiàn)象，主要為結(jié)構(gòu)自身響應(yīng)引起。邊墻中間測點(diǎn)2的應(yīng)力時(shí)程與邊墻與頂板交界處測點(diǎn)1的應(yīng)力時(shí)程明顯不同，存在明顯的二次加載現(xiàn)象，主要為箱涵周邊土體流動變形導(dǎo)致壓應(yīng)力增加引起［9］。由測點(diǎn)1、2、3的豎向應(yīng)力可見，變化規(guī)律存在較大差別。測點(diǎn)1存在較明顯的沖擊峰值，測點(diǎn)2、3則不明顯。不同爆炸距離下測點(diǎn)2的應(yīng)力時(shí)程差別較大，拉壓應(yīng)力的變化與埋深密切相關(guān)。結(jié)合損傷分析可以看出，土中爆炸引起的箱涵破壞主要出現(xiàn)在近爆側(cè)的邊墻和邊墻與頂板、底板交接位置。

4 結(jié) 論

（1）采用SPH法模擬爆炸近域土體，F(xiàn)E法模擬遠(yuǎn)域土體及鋼筋混凝土箱涵結(jié)構(gòu)，對土中爆炸作用下箱涵動力響應(yīng)進(jìn)行了數(shù)值分析，爆后40ms計(jì)算時(shí)間約14h，具有較好的計(jì)算效率，計(jì)算模型較好地揭示了鋼筋混凝土箱涵的動力響應(yīng)，計(jì)算結(jié)果可滿足工程要求。

（2）不同爆炸位置下箱涵結(jié)構(gòu)各測點(diǎn)的爆炸響應(yīng)存在較大差異。土中爆炸引起的箱涵破壞主要出現(xiàn)在近爆側(cè)的邊墻和邊墻與頂板、底板交接位置；值得注意的是，邊墻中間位置受土體變形流動影響，存在較明顯二次加載現(xiàn)象，并且爆炸位置對拉壓變化影響明顯。

［1］Weidlinger P，Hinman E.Analysis of underground protective structures［J］.Journal of Structural Engineering，1988，114（7）：1658－1673.

［2］LU Yong，WANG Zhong－qi，Chong K.A comparative study of buried structures in soil subjected to blast load using 2Dand 3Dnumerical simulations［J］.Soil Dynamics and Earthquake Engineering，2005，25（4）：275－288.

［3］WANG Zhong－qi，LU Yong，HAO Hong，et al.A full coupled numerical analysis approach for buried structures subjected to subsurface blast［J］.Computers and Structures，2005，83（4／5）：339－356.

［4］王吉，王肖鈞，卞梁.光滑粒子法與有限元的耦合算法及其在沖擊動力學(xué)中的應(yīng)用［J］.爆炸與沖擊，2007，27（6）：522－528.WANG Ji，WANG Xiao－jun，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.

［5］杜義欣，劉晶波，伍俊，等.常規(guī)爆炸下地下結(jié)構(gòu)的沖擊震動環(huán)境［J］.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），2006，46（3）：322－326.DU Yi－xin，LIU Jing－bo，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.

［6］Attaway S W，Heinstein M W，Swegle J W.Coupling of smoothed particle hydrodynamics with the finite element method［J］.Nuclear Engineering and Design，1994，150（2／3）：199－205.

［7］王政，倪玉山，曹菊珍.沖擊載荷下混凝土動態(tài)力學(xué)性能研究進(jìn)展［J］.爆炸與沖擊，2005，25（6）：519－524.WANG Zheng，NI Yu－shan，CAO Ju－zhen.Recent advances of dynamic mechanical behavior of concrete under impact loading［J］.Explosion and Shock Waves，2005，25（6）：519－524.

［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.Rotterdam：Balkema A A，1981：141－144.

［9］Baylot J T.Effect of soil flow changes on structures loads［J］.Journal of Structural Engineering，2000，126（12）：1434－1441.