砌体结构长脉宽爆炸荷载损伤等级评估方法*

2021-11-15肖桂仲冯晓伟

爆炸与冲击 2021年10期

曾繁，肖桂仲,3，冯晓伟，黄超，田荣

（1. 中物院高性能数值模拟软件中心，北京 100088；2. 北京应用物理与计算数学研究所，北京 100088；3. 南京理工大学，江苏南京 210094；4. 中国工程物理研究院总体工程研究所，四川绵阳 621999）

近年来，百千吨级当量爆炸工业事故频繁发生，如2013 年4 月17 日美国德克萨斯州化肥厂爆炸事故、2015 年8 月12 日中国天津滨海新区爆炸事故以及2020 年8 月4 日黎巴嫩贝鲁特港口爆炸事故。这类爆炸产生的冲击波持续时间长、危害范围大，导致方圆数公里内的大量建筑结构受损和人员伤亡。研究建筑结构在长脉宽爆炸荷载作用下的损伤评定方法，对结构的抗爆安全设计以及灾后评估具有重要的理论意义和工程价值。

爆炸荷载作用下的损伤程度评估中常用压力-冲量曲线。如何确立压力-冲量曲线，是建筑结构损伤评估方法中的主要研究内容。在构件级，学者们通过实验、理论、数值模拟等手段获得了构件的压力-冲量曲线。Wesevich 等[1]通过236 次不同爆炸荷载作用下的砖墙实验，研究了不同跨度、厚度、钢筋构造方式砖墙的归一化压力-冲量曲线。Ma 等[2]基于单自由度体系归一化的无量纲分析方法，采用最大位移的破坏准则，给出钢混梁的单自由度体系归一化的压力-冲量曲线。Shi 等[3]基于钢筋混凝土结构柱竖向剩余承载力的破坏准则，通过数值模拟方法获得了柱结构的压力-冲量曲线。

而对结构级损伤评估的研究较少。在抗震领域中，陆新征[4]开展了地震响应与结构损伤程度的规律性研究，以层间位移角为判据获得砌体结构破坏等级准则。在爆炸效应领域中，李翼祺等[5]结合部分实验结果，采用超压破坏准则，通过构件的损伤破坏现象，研究了砌体结构破坏等级与空气冲击波超压的规律。类似地，美国化工过程安全中心（CCPS）[6]研究了一般建筑结构的破坏程度与超压的规律。Ding 等[7]采用单自由度方法评估构件损伤程度与分析移除损坏构件后的钢结构非线性力学行为的两步骤评估方法，研究钢框架结构在意外爆炸荷载下的倒塌概率。陶俊林等[8]以单跨钢筋混凝土框架结构为例，依据内爆后的结构毁伤响应数据，采用神经网格反馈算法对爆炸毁伤作用后的受损结构进行生存能力的研究。

目前，研究多集中在构件级的损伤评估模型，并且取得了一系列有价值的研究成果。在结构级损伤评估中，评估指标与外部荷载特征相关联，还未形成统一的评估指标。同时，有学者采用结构构件的解耦评估模型，研究结构级的损伤破坏情况。实际过程中，强爆炸荷载作用下的结构呈现整体变形特征，而局部构件的损伤对整体结构的损伤程度的贡献尚未见报道。因此，为了准确评估建筑结构在强爆炸荷载作用下的损伤破坏程度，有必要开展建筑结构的评估方法研究。

本文中，以砌体结构损伤评估为切入点，提出一种构件损伤加权评估方法，运用冲击波结构毁伤大规模有限元程序，开展长脉宽爆炸荷载下结构损伤高分辨率数值模拟，将提出的评估方法应用于模拟结果的数据分析，给出结构的损伤超压值，并与文献结果进行对比与讨论，验证评估方法的有效性。

1 构件损伤加权评估方法

建筑结构是指在房屋建筑过程中，由各种构件（门框、玻璃、砖墙、梁、板、柱）组成的能够承受各种作用的体系。在长脉宽爆炸荷载作用下，建筑结构的动力学响应表现出整体变形特性，不同类型构件在结构抗爆安全性分析中的作用不同，同时局部构件的损伤程度也反应了整栋结构的损伤情况。基于此，我们提出一种构件损伤加权的评估方法，即依据构件的损伤程度以及在结构中对抗爆安全性的贡献，综合评估建筑结构的损伤等级。具体的表达式为：

式中：D为建筑结构的无量纲损伤程度，下标i为构件编号，Ni为构件权系数，Si为构件的无量纲损伤程度。该评估方法有构件损伤程度和构件权系数两个关键要素。

1.1 构件损伤程度

目前，用于评定构件损伤程度的指标有剩余承载力[3]、裂纹宽度[8]、支座转角及剪应变[9]。在长脉宽爆炸荷载作用下，构件的变形模式多以弯曲变形为主，根据构件支座转角指标评估构件损失程度较合理。美国的国防部[9]和土木工程师协会[10]分别定义了基于弯曲变形模式的混凝土构件设计准则，给出了不同设计标准所对应支座转角的关联：在文献[9]中，定义12°为钢筋混凝土板梁断裂失效的设计水准；而在文献[10]中，定义10°为单配筋钢筋混凝土构件断裂失效的设计水准。对于无配筋的砌体填充墙结构，Wesevich 等[1]根据砖墙实验结果，归纳得到砖墙倒塌破坏所对应的6°支座转角。结合构件的断裂破坏设计水准的极限支座转角θm和变形过程中最大支座转角θu，定义构件损伤程度S为：

对于无配筋砖墙结构，结合文献[1]中的实验结果，给出了砖墙的损伤等级与角度的关系：轻度损伤等级，0°＜θu≤0.3°；中度损伤等级，0.3°＜θu≤1.5°；严重损伤等级，1.5°＜θu≤6°；倒塌等级，θu＞6°；砖墙断裂破坏设计水准的极限支座转角θm=6°。通过式(2)，给出砌体填充墙损伤程度等级：轻度损伤等级，0＜S≤0.05；中度损伤等级，0.05＜S≤0.25；严重损伤等级，0.25＜S≤1；倒塌等级，S＞1。

1.2 构件权系数

在建筑结构的抗爆安全性分析中，将单个构件的可靠度对整个建筑系统可靠度的贡献称为构件权系数。在爆炸荷载外力做功下，建筑结构可以看作传递荷载的基础网络，其中构件负责传递所承担的一定建筑面积上的荷载，同时通过自身变形来吸收外部输入的能量。因此，可以从能量角度出发，给出构件权系数的计算方法。

对于一个建筑结构系统，外力所做功率为：

1.3 建筑结构损伤等级

在爆炸荷载作用下，建筑结构的损伤程度D由构件的损伤程度和权系数决定。而目前，建筑结构损伤等级的评定与D的关系还不清楚。借鉴文献[5]中的建筑结构损伤等级划分标准，将建筑结构损伤等级划分为基本无损坏、极轻度损伤、轻度损伤、中度损伤、严重损伤和倒塌等6 个等级。前2 个等级描述结构中非承重构件（如门框、玻璃）的损伤，后4 个等级描述由局部承重构件的损伤而诱发结构性的损伤。在建筑结构抗爆安全性分析中，承重构件失效，可能导致它所承担的建筑面积上的其他构件失效，从而影响结构的稳定性甚至引起结构的垮塌。因此，我们重点关注如何由承重构件损伤评估建筑结构的损伤等级。以砌体建筑结构为例，按照文献[1]中定义的砖墙损伤程度等级，外推到建筑结构的损伤程度等级：0≤D≤0.05，轻度损伤等级，其中D=0 对应砌体结构门框玻璃易损件完全破坏，而承重结构基本上处于弹性可回复的变形阶段；0.05＜D≤0.25，中度损伤等级；0.25＜D≤1，严重损伤等级；D＞1，倒塌等级。

2 砌体结构动力学响应大规模数值分析

2.1 计算模型与网格划分

自主研发的冲击波结构毁伤程序是一款大规模流固耦合并行程序[11-12]，能够高效高精度模拟爆炸冲击类问题，其中流场分析模块基于JASMIN 框架研发[13]，结构动力学响应分析模块基于JAUMIN 框架研发[14]。运用冲击波结构毁伤模拟程序，建立两榀两开间两层砖混结构的三维有限元模型，如图1 所示。该房屋结构水平跨度均为3 m，层高3 m，包含18 根立柱、16 面承重墙、8 块楼板，立柱和承重墙均由砖块和砂浆堆砌而成，楼板由C30 混凝土和单层双向HRB235 钢筋构成。柱结构截面尺寸均为400 mm×400 mm，砖墙厚度为240 mm。钢筋混凝土楼板厚度为150 mm，水平尺寸为3 m×3 m，体积配筋率为1%；房屋结构共设置了16 个门窗，为了简化计算，建模时不考虑门窗玻璃的影响，仅设置洞口，允许爆炸冲击波直接从门窗洞口中传播。

图1 砌体结构的有限元模型Fig. 1 Finite element model of a masonry structure

在流场冲击波传播分析中，远场建筑结构的爆炸冲击波具备近似平面特性及百毫秒级长脉宽正压区作用时间的特点，因此我们不对炸药进行建模模拟，而在流场边界施加三角形冲击波荷载近似等效百千吨级当量炸药爆炸的远场冲击波环境。流场计算域内的空气采用Euler 网格，精细的网格尺寸100 mm足以保证流场域内冲击波传输的精确模拟。为避免边界位置处的冲击波效应对建筑结构的影响，空气域不仅包裹整个建筑物，而且足够大。因此，在建模过程中，空气域三维空间尺寸为[-25 m, 35 m]×[-35 m,35 m]×[0 m, 30 m]，网格规模达到1.26 亿个，砌体结构放置于流场区域的中心。

在结构动力学分析中，砖墙材料、混凝土、钢筋、地面均采用单点积分实体单元，网格尺寸为25 mm，总共582 万个单元。砖墙材料的数值模型一般有两种建模方法，一类将砖和砂浆作为独立的构成部分，分别用不同的材料建模，另一类按照砌体单元结构特性的等效原则建立等效砌体材料的模型。在爆炸荷载作用下，砖墙材料和混凝土材料将处于大应变、高围压及高应变率的状态。在众多混凝土动态本构模型中，因简明合理的描述和计算程序的适应性，Holmquist-Johnson-Cook（HJC）本构模型在混凝土爆炸问题中应用广泛[15]。针对砌体材料，采用基于HJC 的等效材料模型。针对钢筋混凝土楼板，采用钢混分离式建模，混凝土材料采用HJC 本构模型，钢筋采用等向强化弹塑性材料模型。假定地面为刚性。砖墙和混凝土单元的大变形会导致网格扭曲，为解决网格大变形问题，模型中采用单元侵蚀算法，即当单元的最大主应变达到0.15 或最大剪应变达到0.9 时删除单元。

2.2 材料模型

建立等效砌体材料模型的技术手段一般有3 类。第1 类用简化的几何体表征几何复杂的基本求解单元，在此基础上通过理论推导获得均质化问题的封闭解[16]。第2 类针对周期性复合结构，采用基于代表性体积胞元的平均化理论。考虑砌体基本单元的复杂性，有必要采用有限元方法获取该类问题的数值解[17-18]。第3 类采用细观力学和微观结构模型，考虑砌体结构基本单元的不同内部成分之间的相互作用，对砌体结构进行高度精细化的模拟[19]。本文中采用第2 类方法，先根据砖块和砂浆排布方式建立砌体结构的代表性体积胞元及其精细有限元模型，再根据代表性体积胞元面内的应力平均和应变平均的等效原则建立等效砌体材料模型。

HJC 模型包含屈服面方程、损伤演化方程、状态方程三个方面。其中屈服面方程为：

考虑砂浆和砖容易脆裂、塑性变形小，假定材料的状态方程为线性状态方程。

建立砖墙构件的代表性体积胞元有限元模型，其中砖块的弹性模量、泊松比、屈服强度分别为11 GPa、0.25 和19 MPa，砂浆的弹性模量、泊松比、屈服强度分别为5.8 GPa、0.33 和5 MPa。通过面内单轴压缩测试，获得代表性体积胞元的水平和垂直方向的材料参数，分别为：垂直方向屈服强度为14.8 MPa，弹性模量为10.1 GPa，泊松比为0.3，水平方向屈服强度为14.6 MPa，弹性模量为9.98 GPa，泊松比为0.26。通过代表性体积胞元的面内应力应变关系，给出等效砖墙材料模型的参数，分别为：初始密度为2.4 g/cm3；弹性模量为10.04 GPa；泊松比为0.28，静态压缩强度为14.7 MPa；参考应变率为1 s-1，屈服面参数A、B、N和C分别为0.28、1.85、0.84、0.007，损伤演化方程中D1为0.04，D2为1，断裂起始应变εf,min为0.002 5，静水截断压力1 MPa。

混凝土材料的本构模型选用HJC 本构模型，其状态方程采用经典的三段式，压力与体积应变的关系分为线弹性、空隙压实过程和完全压实3 个阶段，其关系式不再详细阐述。混凝土材料参数分别为[20]：初始密度为2.4 g/cm3；泊松比为0.25，静态抗压强度为30 MPa；参考应变率为1 s-1，参数A、B、N、C、D1、D2和εf,min与砖墙材料一致，静水截断压力为4 MPa，pcrush为13.3 MPa，μcrush为0.000 7，K1为85 GPa，K2为-171 GPa，K3为208 GPa，plock为0.8 GPa，μlock为0.1。

钢筋模型采用等向强化弹塑性材料模型，应变率效应采用Cowper-Symonds 模型，其屈服面方程为：

式中： σ0为材料初始屈服强度，ET为硬化切线模量，εp为等效塑性应变，ε ˙ 为应变率，C、P为应变率参数。钢筋材料的参数分别为：初始密度为7.83 g/cm3，泊松比ν 为0.28，初始屈服强度为235 MPa，切线模量为1.6 GPa，C和P分别为40 和5。

2.3 数值模拟结果分析

在流场分析数值模拟中，将爆炸荷载近似成三角波形的超压时程，作为流场分析的边界条件。固定正压区时间100 ms 不变，改变入射超压峰值pI，选取45、55、60、67、80、95 和142 kPa 等7 种入射超压峰值工况，开展远场冲击波在建筑结构的传播规律研究。以入射超压峰值80 kPa 为例，采用大规模可视分析平台TeraVAP[21]对计算结果进行后处理分析，对爆炸波在空气中的传播以及引起的空气压力分布进行了可视化显示，如图2 所示。在没有接触建筑物前，冲击波以平面波的形式往前传播（见图2(a)～(b)）；爆炸冲击波在遇到建筑物立面时发生反射，并且反射冲击波超压远比入射冲击波提高数倍（见图2(c)）；爆炸冲击波继续向前传播，逐步包裹整个建筑结构，同时冲击波通过门窗进入建筑物内部（见图2(e)～(f)）；由于爆炸冲击波在建筑内部的来回反射，反射波与入射波、反射波与反射波多次叠加，使建筑内部的冲击波非常复杂。在浪潮大规模计算服务器上计算该工况，模拟物理时间0.4 s，计算网格1.26 亿个，采用112CPU 核，耗时4.1 h。其他超压工况的结果不再详细阐述。

图2 超压峰值80 kPa 下流场冲击波传输的模拟结果Fig. 2 Numerical simulation results of shock wave propagation under the overpressure peak of 80 kPa

为了验证流场计算的精度，抽取了不同入射超压下建筑物迎爆面中心点位置的反射超压峰值，并与理论公式结果进行对比。由文献[22]，正反射条件下入射超压与反射超压得关系为：

式中：p0为未扰动空气压力，pR为反射冲击波阵面超压。由图3 可知，建筑迎爆面上反射超压峰值与经验式(12) 的计算结果相当，相对误差在5%以内，基本满足了流场的计算精度。相比经验式，流固耦合计算同时精确地捕捉了冲击波在建筑结构内、外表面的反射、绕射等效应，为结构的变形分析提供准确的冲击波荷载。

图3 迎爆面的反射超压峰值Fig. 3 The reflected overpressure peaks on the blasting face

在结构动力学响应分析过程中，先对结构施加重力荷载，当结构达到静力平衡后，将流场压力映射至建筑结构表面，开展结构的动力响应数值模拟，如图4 所示。当pI=67 kPa 时，砌体建筑结构整体变形不大，未发生垮塌现象，局部区域出现较大的弯曲变形和损伤，如迎爆面处的柱和砖墙（见图4(a)）。当pI=80 kPa 时，砌体结构迎爆面位置一层的柱和砖墙构件发生损伤破坏，其他位置构件出现不同程度的损伤，导致结构整体承载力显著下降，结构在重力作用下发生竖向垮塌失效（见图4(b)）。当pI= 142 kPa 时，柱、砖墙、楼板构件均发生损伤破坏，已破坏的构件在动压作用下发生抛射现象，未破坏的结构在重力作用下发生竖向倒塌破坏，砌体结构表现一种内爆式整体性解体失效的破坏模式（见图4(c)）。在浪潮大规模计算服务器上计算典型工况，模拟物理时间1.5 s，计算网格582 万个，采用280CPU 核，耗时46.6 h。其他超压工况的结构响应结果不再详细阐述。

图4 砌体结构的损伤破坏形态Fig. 4 Damage and fracture morphology of the masonry structure

3 砌体结构损伤评估分析

砌体结构的损伤评估由构件权系数Ni和损伤程度Si决定。在考虑结构完整性前提下，通过面荷载作用下的建筑结构准静力平衡分析，给出每个构件的权系数。在砌体结构损伤数值分析过程中，根据每个构件端部和跨中位置的位移曲线，统计每个构件的损伤程度。

3.1 构件权系数

首先，对结构施加自重荷载，直到结构静力平衡。然后，在建筑结构外表面施加面荷载，维持一定时间后，直至结构达到新的静力平衡，此时结构的动能与应变能相比可以忽略不计。在静力平衡过程中，当结构所有节点的速度低于0.1 m/s 时，便认为结构达到新的静力平衡。通过计算构件所属网格单元的应变能，进而统计该构件在整个变形过程中的总应变能，对每个构件的应变能求和获取结构的总应变能。在结构达到新的静力平衡时，外力荷载做功等于结构总应变能。最终，由式(8)计算每个构件的权系数。

表1 为砌体结构的构件权系数，其中每个构件对应唯一的标识，标识由构件名称、编号和几何中心点坐标组成，例如柱1-(0.2, 0.2, 1.5)。由表可见：（1）同楼层的同种构件权系数相当，如1 层的8 根外柱（柱5 为内柱）；（2）不同楼层的同种构件权系数相差较大，如1 层柱1 与2 层柱9；（3）同楼层不同构件的权系数相差较大，其中柱构件系数最高、砖墙其次、楼板最小。

表1 柱、砖墙、楼板的权系数Table 1 The weight values of columns, brick walls and floors

根据爆炸荷载的正压区作用时间与结构的自振周期，爆炸荷载可分为冲量荷载、准静态荷载和动力荷载等3 类荷载。在冲量荷载及动力荷载区域内，基于静载计算构件变形能构件的权系数的计算方式关键在于如何设计类爆炸的静态荷载。在准静态荷载区域内，爆炸荷载可近似等效为静载。目前，作用于低层砌体结构的长脉宽爆炸荷载对象，属于准静态荷载类型，依据静载作用下构件的变形能计算构件权系数，能客观反映构件对整个建筑系统的可靠度。

3.2 构件损伤分析

在砌体结构损伤数值分析过程中，先提取每个构件端部和跨中位置的位移曲线；根据构件初始构型和当前构型，扣除构件刚体位移，计算构件的变形信息；再依据构件的变形结果，获取构件的端部支座转角；根据构件损伤等级准则，给出构件的损伤程度Si。砖墙的损伤准则采用文献[1]中的无配筋砖墙损伤准则，砖柱构件近似采用无配筋砖墙损伤准则，钢混楼板采用文献[9]中的设计规范。在文献[9]中，钢筋混凝土板构件损伤等级与支座转角的关系为：轻度损伤等级对应0°～2°，中度损伤等级对应2°～6°，严重损伤等级对应6°～12°。

图5 为80 kPa 下柱构件的支座转角曲线，其中图5(a)～(b)分别对应1～2 层。1 层柱构件的最大支座转角均超过倒塌破坏的转角6°，而2 层柱构件的大部分最大支座转角小于1°，仅迎爆面上中柱的最大支座转角达到3.3°。因此，通过柱构件的最大支座转角定量判断，1 层柱结构的倒塌破坏等级和2 层柱结构损伤等级介于轻度和中度损伤之间，个别为严重损伤。根据构件损伤等级准则，给出每个柱构件的损伤程度Si，如柱11 的S11=0.55。对于已经破坏的构件，在动压和惯性作用下，其支座转角会迅速增加。为此，我们定义已破坏构件的支座转角的上限值。考虑一种极端情况，当底层柱构件均发生倒塌破坏失效而其他构件完好无损时，该建筑结构发生垮塌（D=1）。由1 层柱权系数和式(1)，计算上限值为15°。

图5 柱构件的支座转角曲线Fig. 5 Rotation curves of the support of column components

进一步地，根据构件支座转角曲线，统计每个构件的损伤程度Si，见表2。由表可见：（1）底层构件的损伤破坏程度高于1 层构件的损伤程度；（2）不同类型构件的损伤程度不一致，大部分柱、砖墙发生严重损伤、甚至倒塌破坏，楼板发生轻度损伤或中度损伤。根据构件的权系数（见表1）和损伤程度（见表2），统计80 kPa 下建筑结构的损伤程度D=1.31，对应倒塌毁伤等级。砌体结构损伤程度分析与数值模拟结果（见图4(b)）基本一致，这说明基于支座转角的评定准则能较好反应结构构件的损伤程度。其他超压工况的构件损伤分析结果不再阐述。

表2 超压峰值80 kPa 下柱、砖墙、楼板的损伤程度Table 2 The damage degree values of columns, brick walls and floors under the overpressure peak of 80 kPa

3.3 砌体结构毁伤评估结果

根据构件的权系数（见表1）和损伤程度（见表2），统计7 种入射超压峰值下砌体结构的损伤程度D，如图6 所示，图中插图为砌体结构的位移云图。在有限元模拟中，很难通过对砌体结构施加合适的荷载使它的损伤程度恰好落在两个破坏等级的分界线上，如D=1.0, 0.25, 0.05。因此，对数据点进行线性回归分析，得到三段式的拟合公式：

图6 砌体结构损伤的评估结果Fig. 6 Results of damage assessment of the masonry structure

由式(13)，给出损伤等级分界线的超压值：当pI=14.13 kPa 时，D=0；当pI=29.80 kPa 时，D=0.05；当pI=63.08 kPa 时，D=0.25；当pI=73.87 kPa 时，D=1.0。对应D=0 情况，砌体结构的承重构件位于弹性变形阶段，而非承重构件（门框玻璃）完全破损。门框玻璃偶尔开裂或震落所对应的超压2 kPa[5-6]，可作为砌体结构基本无损坏的超压阈值。

表3 为损伤等级临界超压值以及与文献[5-6]中结果的对比。由表可见：（1）文献[5]与[6]的损伤等级分类基本一致，超压值上下限范围大体相同，但相同损伤等级的超压值存在一定差异，如中度损伤的超压值的相对偏差为-27.54%；（2）基于数值模拟结果和评估方法预测的损伤等级临界超压值，介于文献[5]与[6]之间，如极轻度损伤、中度损伤、严重损伤及倒塌超压值；（3）预测的6 个损伤等级的超压值与文献[5]、[6]的超压值的相对误差在-16.88%～26.16%之间。对比结果验证了构件加权评估方法的有效性，在一定程度上能预测建筑结构的超压值。

表3 砌体结构损伤程度和空气冲击波超压值Table 3 The damage degrees of the masonry structure and the overpressure values of the air shock wave

Brasie 等[23]获得了砌体结构破坏的超压区间为41～69 kPa。Perry 等[24]提出，34 kPa 超压导致砌体结构75%的外墙倒塌，69 kPa 超压摧毁传统的砌体结构。Crowl[25]描述了厚25～32 cm 砖墙结构失效对应的超压区间为48.3～55.2 kPa。Kinney 等[26]给出了20～30 cm 空心砌块砖墙破碎对应的超压区间为13.8～20.7 kPa。综合文献数据，得知砌体结构的不同程度的损伤破坏的超压区间为13.8～69.0 kPa，与通过构件加权损伤的评估方法获得的砌体结构损伤破坏的超压区间14.3～78.3 kPa 基于一致，进一步验证了评估方法的有效性。