土中爆炸作用下箱涵动力响应的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.