基于LS-DYNA的大空间柱壳结构爆炸波压力场分析

高轩能，刘 颖，王书鹏

(华侨大学 土木工程学院，福建泉州 362021)

大空间结构在候车(机)楼、影剧院、体育馆和博物馆等公共建筑中得到广泛的应用，已成为我国公共建筑的主要形式之一。由于人员活动频繁密集，这类建筑也极易成为恐怖爆炸袭击的目标，造成重大的人员财产损失和恶劣社会影响。如2003年12月5日俄罗斯皮亚季戈尔斯克市车站的恐怖爆炸，造成近50人死亡150多人受伤;2004年1月23日我国郑州火车站售票厅内的炸弹爆炸，导致2人死亡3人受伤;2004年5月9日俄罗斯车臣共和国格罗兹尼市迪那摩体育场的爆炸，造成32 人当场死亡近 100 人受伤［1，2］。因此，对大空间建筑结构进行防爆和抗爆研究，具有迫切的现实意义。

对于大跨空间结构在爆炸作用下的研究，国外迄今仅见Raftoyiannis［3］关于多点悬吊屋盖在爆炸作用下的动力性能研究的1篇论文，国内的研究也不多，尚处于探索阶段。李海旺等［4］对单层球壳结构在球顶集中冲击荷载作用下的动力反应及稳定性进行了试验与数值分析。范峰等［5］对K8型单层球壳在球顶集中冲击荷载作用下的动力特性、失效模式和失效机理进行了数值模拟计算，建立了杆件破坏形式与网壳失效模式之间的对应关系。周岱等［6］对拉索-网壳结构的动力特性和非线性动力反应进行了研究。宋卫东等［7］对机库内爆炸效应进行了3维数值模拟分析，计算了冲击波超压峰值与爆炸点距离的关系。高轩能等［8］对大跨空间柱壳结构在爆炸荷载下的动力效应及防爆措施等进行了研究，探讨了爆炸冲击波在壳体顶部的聚集效应和对结构动力性能的影响以及泄爆方法等。但对大空间结构进行爆炸作用下的动力响应研究，首先需要确定爆炸冲击波压力场在结构上的分布，因大空间结构自振频率分布密集、振型复杂和计算量大等困难，目前的研究尚少。为此，本文基于有限元程序ANSYS/LS-DYNA，对大空间柱壳结构进行爆炸模拟实验，在确定合适的模型及参数后对其进行爆炸作用下的动力响应计算。引入POD法［9］，利用提取的结构表面有限测点的冲击波超压数据，分析大空间柱壳结构的爆炸冲击波压力场分布以及结构几何参数等因素对冲击波分布的影响。

1 爆炸冲击波及其经验计算公式

图1 爆炸冲击波示意图Fig.1 Blast shock wave sketch map

化学炸药爆炸后瞬间便转变成高压(104MPa～3×104MPa)和高温(3.5 ×103℃ ～4×103℃)状态的气态爆轰产物。气体急剧膨胀并迫使周围的空气离开它原来占据的位置，于是在此气体的前沿便形成一压缩空气层，即爆炸冲击波(如图1 所示)［10］。爆炸的瞬时能释放巨大的能量并产生各种破坏效应，但对建筑物破坏力最强、破坏区域最大的还是爆炸冲击波的破坏效应。从一些典型的恐怖爆炸事件看，建筑物的爆炸破坏情况主要与爆炸冲击波的大小、持续时间、爆炸源的位置和建筑物本身的建筑结构特性等因素有关，还与恐怖爆炸袭击的方式等其它诸多因素有关。由于爆炸过程的短暂性和复杂性，爆炸冲击波压力很难通过简单的理论分析得到。理论与实验研究发现，炸药在空气中爆炸形成的空气冲击波，阵面上的压力并不取决于装药的绝对尺寸和质量，而完全取决于离爆炸地点的距离与装药半径的比值、该炸药的比能和周围空气的压力。空气冲击波的峰值超压ΔPf是指冲击波阵面上的峰值压力Pf减去空气中的大气压力P0(标准大气压)，即:

根据我国国防工程设计规范，化爆冲击波超压及正压作用时间按如下公式计算［11］:

式中:W为TNT当量，单位kg;R为计算点距爆心的距离，单位m;H为炸药爆炸时离地面的高度，单位m;ΔPf为爆炸冲击波超压峰值，单位MPa，t+为冲击波超压正压作用时间，单位s。

当炸药在地面爆炸(按刚性地面考虑)时，由于地面的阻挡，空气冲击波不是向整个空间传播，而只向一半无限空间传播，被冲击波带动的空气量也减少一半，可看作两倍的装药在无限空间爆炸，此时，爆炸冲击波超压计算公式为:

亨利奇(Henrych)［12］计算公式则为:

为简化计算，本文采用线性冲击波函数［13］:

式中:ΔP-为冲击波负压，t-为负压作用时间。

2 大空间柱壳结构的ANSYS/LS-DYNA爆炸模拟实验

2.1 LS-DYNA简介

LS-DYNA程序是由 Hallquist于1976年在美国Lawrence Livermore National Laboratory主持开发完成，是功能齐全的通用大型非线性程序。LS-DYNA的分析功能非常强大，有140多种材料模式库，包括工程材料、炸药、推进剂、粘性流体及自定义材料等，以及各种单元模型和50多种接触方式。它以Lagrange算法为主，兼有ALE和Euler算法;以显示求解为主，兼有隐式求解功能;以结构分析为主，兼有热分析、流体-结构耦合功能;以非线性动力分析为主，兼有静力分析功能。LS-DYNA独特的算法能十分准确地分析工程结构承受高速碰撞、爆炸、金属成形等高速高度非线性问题，在工程结构抗爆特性研究、抗爆设计以及爆破拆除等方面有着非常重要的用途［14］。

2.2 材料与结构模型

炸药采用MAT-HIGH-EXPLOSIVE-BURN材料模型和JWL状态方程。状态方程定义的爆炸冲击压力为［14］:

式中:E0为单位体积炸药的初始内能，V为相对体积，A，B，R1，R2，Ω 为基于试验的输入参数。TNT 炸药为0.5 m×0.5 m×0.5 m 立方体。炸药的材料输入参数如表1所示，其中ρ为密度，D为爆速，PCJ为爆压，V0为初始相对体积。

表1 炸药的材料参数Tab 1.Material parameters of dynamite

空气采用MAT-NULL材料模型和线性多项式状态方程EOS-LINEAR-POLYNOMIAL。即:

式中:E为单位体积内能，ρ为质量密度，ρ0为参考质量密度。线性多项式状态方程描述空气时遵守Gamma定律，空气的材料输入参数如表2所示。

表2 空气的材料参数Tab 2.Material parameters of air

结构模拟计算模型以公共大跨度建筑中常用的大空间柱壳结构为对象，应用LS-DYNA有限元程序建立模拟大空间结构室内爆炸的计算模型。大空间结构设为三向网格单层柱面钢网壳，如图2所示，结构尺寸为40 m×20 m×15 m，矢高为3 m，ALE空间尺寸为42 m×21 m×19 m，矢跨比F/B=1/5，周边固支，支座下围护墙体高H=10 m。网壳结构的杆件均为Q235钢无缝钢管，焊接球节点，其中轴向直杆截面为φ89 mm×4 mm，斜杆截面为φ140 mm×6 mm。为简化计算，利用对称性取1/2模型进行分析，除对称面外，边界面为透射边界以模拟在无限区域的爆炸。本文着重考虑人体或车载炸药等对大空间结构的爆炸袭击，爆炸点高度取为离地面高1.2 m。爆炸荷载下不考虑屋面围护材料等质量，仅考虑杆件及节点自重。计算模型中，屋面壳体结构、地面和墙体视为刚体，以便考虑空气冲击波压力场分布。考虑5种物质的相互作用，选用3种单元。空气和炸药均采用六面体实体Solid164单元，地面和四周围护墙体采用Shell163单元，结构杆件采用Link160单元。

图2 大空间钢柱壳结构计算模型Fig.2 Model of large-space steel cylindrical reticulated shell

2.3 实验结果及其验证

为验证模型和材料参数取值的可靠性和正确性，应用LS-DYNA对炸药在无限空中爆炸进行了数值模拟，空气尺寸为12 m×12 m×12 m，炸药尺寸为0.2 m×0.2 m ×0.2 m，位于空气中心。算法采用多物质ALE算法，有限元网格大小取为0.1 m ×0.1 m ×0.1 m，计算模型共划分8节点六面体单元216 000个，其中炸药为1个单元，空气为215 999个单元。边界面为透射边界。模拟计算结果与经验公式［11，12］计算结果的比较如图3所示。图3表明，本文模型得到的冲击波超压与3个经验公式结果都很接近，而与亨利奇公式［12］的结果最为相符，在超压持续时间上，则本文模拟计算结果与亨利奇公式几乎完全吻合，同时还可看出，我国设计规范公式与亨利奇公式计算结果基本符合，略为偏大。

另一方面，炸药爆炸后，因爆炸点距离地面较近，冲击波碰撞地面，地面处空气质点速度骤然变为零，使质点急剧堆积，压力和密度骤然升高，达到一定程度后，质点就要向相反方向运动，即形成反射冲击波。反射冲击波超压大大增强，所以大空间结构的一部分构件会受到冲击波的两次冲击作用，并且第2次对结构的冲击作用远远大于第1次。图4(a)为爆炸点上方h=13 m处的冲击波压力曲线，可以看出该点受冲击波2次冲击作用，第1次为入射冲击波，根据文献［12］该处的入射冲击波超压为0.182 MPa(亨利奇公式)和0.205 MPa(我国规范公式)，本文有限元结果为0.217 MPa，两者分别相差19%和6%，表明本文结果与经验公式符合较好。第2次为地面反射冲击波，由图4(b)质点加速度图可以看出，该处质点在入射冲击波作用后，加速度由最大变为零再到负值，随后受到反射冲击波的作用，加速度再次升到最大，由于反射冲击波较入射冲击波大，质点第2次获得的加速度也会比第1次大(3.25倍)，之后加速度从最大值变小，期间冲击波压力升到最大，然后恢复到标准大气压。

如果结构是封闭的且在冲击作用下不破坏，则冲击波受到阻挡，空气质点会再次在屋盖表面区域堆积，屋盖表面压力会随之增大，但速度较小如图5所示。

图6为马赫反射区水平距爆炸点10 m处的冲击波压力曲线，可以看出，该处只承受马赫反射波1次冲击作用。参考文献［11］，地爆超压仍采用文献［12］的公式计算，TNT当量分别取为1W和2W，计算结果分别作为地爆超压的上下限值，分别为0.271 MPa和0.447 MPa，本文结果为 0.341 MPa。

图5 屋盖表面处空气质点压力及速度Fig.5 Air pressure and velocity near the roof

图6 马赫反射区距爆心10 m处的空气压力Fig.6 Air pressure in Mach reflection region(10 m away from the explosion center)

从以上模型计算结果的分析可知，计算结果与理论分析相符，各测点的冲击波超压和空气质点速度等数据与经验公式的计算结果也符合较好，表明模型及参数选取可靠和正确，计算结果可信。

3 数值模拟结果与分析

对于图2所示大空间柱壳结构，由于炸药距离地面较近，爆炸后冲击波受到地面的阻挡，反射冲击波超压大大增强，因此，结构的一部分构件会受到冲击波的2次冲击作用，并且第2次的冲击作用会远远大于第1次。由于作用于结构表面各点的爆炸荷载时间不同步，大小不相同，使结构表面各点上的冲击波压力形成了时空差异。冲击波压力场的分布是进行大空间结构防爆和抗爆研究的基础，为克服计算上的困难，本文引入本征正交分解法(POD法)［15］，提取结构表面n个点处的冲击波超压数值并将其进行POD分解，即可得到结构的爆炸冲击波压力分布图。

按前述建立的有限元计算模型对不同跨度和矢跨比的大空间柱壳结构进行爆炸模拟计算后，对爆炸冲击波超压进行POD法求解，得到的冲击波压力场第1阶本征模态分布及其超压幅值-时间曲线如图7和图8所示。图中x轴表示结构跨度方向(0 m～20 m或40 m)，y轴表示结构长度方向(0 m～20m，对称取半计算)，结构高度H=10 m，长度L=40 m，炸药0.5 m ×0.5 m×0.5 m。对于跨度B=20 m，矢跨比F/B分别取 1/2、1/4、1/5，对于跨度B=40 m，矢跨比F/B分别取 1/4、1/5、1/8，其它参数如前所述。

从图7、图8中可以看出，冲击波本征模态分布随结构的矢跨比减小而变得分散，时间坐标的峰值随矢跨比减小而减小。跨度为20 m时，矢跨比F/B=1/2结构的第1阶本征模态相应的超压幅值时间曲线与另两种矢跨比的曲线差别较大，曲线有两个波峰，其第2个波峰较第1个波峰小很多，曲线的负压区不明显，而矢跨比F/B=1/4和1/5的超压幅值时间曲线出现了3次波峰，且负压区明显。跨度为40 m时，第1阶本征模态相应的超压幅值时间曲线均出现两个波峰，第1个波峰随矢跨比减小而减小，而第2个波峰随矢跨比减小而增大。

4 结论

综上所述，基于有限元软件ANSYS/LS-DYNA建立的大空间柱壳结构爆炸作用数值模拟计算模型，与Henrych经验公式计算结果符合较好，表明计算模型和材料参数取值正确而适用，可用于大空间柱壳结构的爆炸动力响应分析。引入POD法，利用提取的结构表面有限测点的冲击波超压数据，解决了冲击波荷载的时空差异性及结构表面压力场分布问题。通过对不同跨度和不同矢跨比的大空间柱壳结构在爆炸荷载作用下的冲击波压力场的数值模拟计算及分布规律的研究，可知结构的矢跨比等几何参数变化对大空间柱壳结构的冲击波压力场分布影响很大，封闭大空间结构的顶盖处可产生冲击波压力场积聚现象，不利于结构的防爆和抗爆，所得结果给大空间柱壳结构在爆炸荷载作用下的动力特性和结构的抗爆设计提供了基础。需要说明的是，由于所探讨问题的复杂性，本文研究仅仅是初步的，对于大空间柱壳结构在爆炸荷载作用下压力场分布的精确分析，尚需考虑材料的非弹性和结构的大变形影响，做更深入的研究。

［1］钱七虎.反爆炸恐怖安全对策［M］.北京:科学出版社，2005.

［2］高轩能，王书鹏.大空间柱面网壳结构在爆炸荷载下的动力响应［J］.振动与冲击，2009，28(10):68-73.

［3］Raftoyiannis I G，Spyrakos C C，Michaltsos G T.Behavior of suspended roofsunderblastloading［J］. Engineering Structures，2007，29:88-100.

［4］李海旺，郭 可，魏剑伟，等.撞击载荷作用下单层球面网壳动力响应模型实验研究［J］.爆炸与冲击，2006，26(1):39-45.

［5］ Fan F，Wang D Z，Zhi X D，et al.Failure modes for singlelayer reticulated domes under impact loads［J］.Transactions of Tianjin University，2008，14:545-550.

［6］周 岱，刘红玉，李春祥.拉索-网壳结构的动力特性和非线性动力反应［J］.振动与冲击，2002，21(1):7-11.

［7］宋卫东，吴开腾，宁建国.机库内爆炸效应的数值模拟研究［J］.兵工学报，2007，28(11):1329-1334.

［8］高轩能，王书鹏，江 媛.爆炸荷载下大空间结构的冲击波压力场分布及泄爆措施研究［J］.工程力学，2010，27(4):226-233.

［9］张建胜，武 岳，沈世钊.单层网壳结构风振响应的模态耦合效应分析［J］.振动与冲击，2006，25(6):39-42.

［10］张守中.爆炸基本原理［M］.北京:国防工业出版社，1988.

［11］ 李翼祺，马素贞.爆炸力学［M］.北京:科学出版社，1992.

［12］Henrych J等，熊建国译.爆炸动力学及其应用［M］.北京:科学出版社，1987.

［13］林大超，白春华，张 奇.空气中爆炸时爆炸波的超压函数［J］.爆炸与冲击，2001，21(1):41-46.

［14］ LS-DYNA keyword user's manual［Z］.Version 971，LSTC，March，2005.

［15］陈 波.大跨屋盖结构等效静风荷载精细化理论研究［D］.哈尔滨:哈尔滨工业大学，2006.