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

崔 溦，宋慧芳，張社榮

（天津大學水利工程仿真與安全國家重點實驗室，天津 300072）

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

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

1 SPH與FE耦合算法

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

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

2 材料模型

2.1 炸 藥

炸藥產物狀態方程用JWL方程描述

式中：v＝1／ρ，為比容，A、B、R1、R2和ω為常數，對于常規爆破，他們可以由動力實驗確定。

TNT材料參數分別為：ρ＝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 土

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

式中：p為靜水壓力，ρ0為初始密度，e為比內能，Γ為Grüneisen參數。

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

2.3 混凝土

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

在硬化階段后，混凝土附加塑性應變導致的材料損傷和強度降低采用下式

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

在該模型中，狀態方程采用p－α模型描述

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

混凝土模型的主要參數分別為：ρ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 鋼 筋

鋼筋采用率相關的John－Cook彈塑性模型。屈服應力為

式中：Y0為初始屈服強度，εp為有效塑性應變為歸一化的有效塑性應變，B、C、n、m為材料常數。TH為相應溫度，TH＝（T－Tr）／（Tm－Tr），Tm為熔化溫度，Tr為環境溫度。式（8）中的第一部分表示應變硬化效應，第二部分和第三部分分別表示應變率和溫度效應。

鋼筋材料參數分別為：ρ＝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 接觸面模型

根據實驗結果［8］，當土與結構物的接觸面比較粗糙時，接觸面應力超出土的最大剪應力時發生破壞，并且高壓力時接觸面的強度特性與土體屬性非常接近。因此，將土與混凝土箱涵間接觸面按完全連接考慮。

2.6 邊界條件

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

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

3 數值模型和結果分析

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

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

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

3.1 混凝土損傷

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

3.2 箱涵響應

不同爆炸距離下箱涵變形如圖4所示?？梢钥闯?，不同測點變形時程差別明顯。由測點1、2的水平變形可見，不同爆炸距離下變形趨勢基本一致，邊墻中間位置水平變形最大，達12.5mm，出現在最近爆炸距離情況，殘余變形為不可回復的塑性變形，而邊墻與頂板交界處變形較小。由測點1、3的豎向變形可見，不同爆炸距離下變形規律基本一致。在爆炸作用下，箱涵整體呈向爆炸位置傾斜的趨勢。

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

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

圖5 （a）箱涵測點1的應力Fig.5 （a）Stresses at point 1of box culvert

圖5 （b）箱涵測點2的應力Fig.5 （b）Stresses at point 2of box culvert

圖5 （c）箱涵測點3的應力Fig.5 （c）Stresses at point 3of box culvert

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

4 結 論

（1）采用SPH法模擬爆炸近域土體，FE法模擬遠域土體及鋼筋混凝土箱涵結構，對土中爆炸作用下箱涵動力響應進行了數值分析，爆后40ms計算時間約14h，具有較好的計算效率，計算模型較好地揭示了鋼筋混凝土箱涵的動力響應，計算結果可滿足工程要求。

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

［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］王吉，王肖鈞，卞梁.光滑粒子法與有限元的耦合算法及其在沖擊動力學中的應用［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］杜義欣，劉晶波，伍俊，等.常規爆炸下地下結構的沖擊震動環境［J］.清華大學學報（自然科學版），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］王政，倪玉山，曹菊珍.沖擊載荷下混凝土動態力學性能研究進展［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.