地震作用下群桩体系桥梁连续倒塌模拟研究

陈皞， 周敉， 安风明， 朱国强

(1.中交第一公路勘察设计研究院有限公司， 陕西 西安 710075； 2.长安大学 旧桥检测与加固技术交通行业重点实验室)

1 引言

由于中国地理位置的特殊性，位于环太平洋地震带和欧亚地震带的交汇处，中国地震频繁震害严重。据相关数据显示，1900年至今，发生在中国的6级及以上地震就有800余次。桥梁基础是桥梁的5大构件之一，对确保桥梁的抗震性能和震后的快速修复具有重要意义，分析大量震害资料可以得出：群桩基础是在地震灾害发生时比较容易破坏的结构，发生病害有桩-土分离出现间隔、桩基础移动较大、桩和柱相接的地方开裂、桩身发生弯曲剪切破坏、桩基础和承台连接处断裂等，就具有群桩基础的高桩承台桥梁而言，目前的研究主要是在静力和稳定分析方面，而对于桥梁的地震失效原理以及相应的冗余储备研究仍略显单薄，很难保证在符合抗震要求的前提下桥梁不发生破坏倒塌。

当前研究结构的连续倒塌过程相对成熟的方法和手段包括：结构连续倒塌数值仿真、基于地震风险评估的结构连续倒塌研究和结构连续倒塌模型试验研究。桥梁的地震倒塌研究是很复杂和前沿的课题，对主流的桥梁动力有限元程序而言，功能上有欠缺，还无法精确实现对桥梁连续倒塌的数值模拟，部分软件基于散体有限元可以计算简单结构的倒塌问题，但在计算效率上仍有很大的改进空间。目前国内外在地震作用下的桥梁连续倒塌分析领域对桥梁性能评价和地震动与动力响应的关系等方面研究较多。范振华等对一座连续变截面箱梁进行计算分析，通过算例对钢构桥抗倾覆倒塌性能进行深层次的论证分析，得到相关的评价衡量准则；叶爱君等将群桩作为整体构件进行研究，得到最后结论：对群桩基础施加不同的作用，得到的群桩基础的相关评价指标数值都比较低；彭卫兵等使用OpenSees对圬工拱桥的连续倒塌分析研究，得出如果把倒塌构件模拟为发生小应变的单元会使结果和实际结果不太相符。同时对TLFD8nBrick单元进行了二次开发，弥补几何非线性对结构模拟的不足，最终对倒塌过程进行了有效的模拟；Uwe STAROSSEK指出了目前设计方法在预防结构发生连续倒塌方面的短板，同时提出了一种基于分类隔震的对桥梁结构非常适用的设计方法。

该文在常规有限元动力分析的基础上针对倒塌问题进行扩展研究。首先清楚界定结构构件的临界倒塌状态，然后使用OpenSees程序对群桩基础连续倒塌过程进行仿真模拟，对OpenSees中的Remove命令可以更新分析模型功能进行深入研究，以解决结构的大变形和大转动相关问题，同时解决传统有限元数值不易求解的问题，完成基础失效准则的建立，对桥梁的连续倒塌过程进行更加合理的模拟。该文把桥梁构件模拟为OpenSees中的非线性梁柱单元，非线性梁柱单元都带有附属节点，在模拟群桩基础释放内力和结构逐步失效破坏过程中需要对附属节点进行处理和控制，同时对破坏杆件进行删除并实时更新模型，完成更新后有限元模型的地震响应分析。

2 基于OpenSees的桥梁结构连续倒塌分析模型

2.1 结构概况

该文选用的桥梁结构见图1，结构选用二类场地，结构承受上部结构质量3 000 t和承台质量400 t，桥梁基础的桩基布置选用3×1的钻孔灌注桩，每根桩基直径为2 m，结构中使用HRB335钢筋，选用砂性土，桩和土之间的相互作用通过等效嵌固模型模拟(图2)，通过输入 Saratoga - Aloha Ave波、El-Centro波以及 Gilroy Array 4#波3种不同地震波对结构的地震响应进行对比，结果如图3～5所示。

图1 桥梁结构模型(尺寸单位：m)

桥墩桩径为4 m，桥墩长度为15 m。基础柱的直径等于2 m，长30 m，邻柱距离和基础等效嵌固长度都设置为3倍的桩径，桩基根数有3个。桥墩材料密度2 600 kg/m3，弹性模量3.45×107kPa，承台材料密度2 600 kg/m3，桩基础弹性模量3.0×107kPa，截面上钢筋根数40根。

图2 群桩基础等效嵌固和节点示意图

图3 Saratoga-Aloha Ave波(地震波1)加速时程曲线

图4 El-Centro波(地震波2)加速时程曲线

图5 Gilroy Array 4#波(地震波3)加速时程曲线

2.2 群桩基础OpenSees模型

动力分析中采用的力、时间和位移单位为kN、s和m。选用OpenSees软件中的Concrete02和Steel02材料，结构阻尼比为0.05。桩基左到右编号依次为11、22、33。建立模型时用集中荷载模拟上部结构，用弹性梁柱单元模拟桥墩和桩基础连接，用集中荷载和弹性梁柱单元模拟承台，采用塑性铰柱单元模拟桩基础，同时带有附属节点。集中荷载主要包括桥跨结构、桥面铺装、栏杆、照明灯具等重力荷载。

基于柔度法的梁柱单元是把单元划为若干节段，然后根据线性插值函数可以得到节点处的截面力。

(1)

先把单元力通过式(1)转化成截面力，再将截面力通过柔度矩阵转变为截面变形。通过力和变形关系得到截面的抗力和切线刚度，单元柔度矩阵由截面的柔度矩阵沿长度进行积分，得到：

(2)

(3)

式中：b(x)为差值函数矩阵；bT(x)为力差值函数矩阵；r(x)为柔度矩阵。

基于柔度法的梁柱单元在模拟弯曲梁柱构件时可以取得非常好的效果，同时迭代收敛速度比较快。

端部塑性铰单元首先要设定塑性铰的长度，在设定的塑性铰长度范围内限制其非线性，对于基于柔度法的塑性铰单元，单元力和变形在局部坐标系下的关系为：

(4)

其中：e为任意截面的变形向量。塑性铰梁柱单元比其他单元求解速度都要快，因为只要使塑性区收敛就可以。

该文使用带附属节点的端部塑性铰梁柱单元，端部附属节点采用集中质量模型，将单元的质量集中在附属节点上。通过 remove mp、equal DOF 和 Zero-Length Element 的组合命令在OpenSees中对多点约束进行处理，实现构件的逐步失效。过程会出现同位置处有多个节点的问题，使求解变得困难，可以通过多点约束的处理方法进行静力凝聚，如式(5)：

(5)

式中：um为控制节点的位移；Δum为控制节点的位移增量；us为附属节点的位移；Δus为附属节点的位移增量；K、C、M分别为质量矩阵、阻尼矩阵、刚度矩阵；P(t)为节点外荷载矢量，U为节点位移矢量。

2.3 结构失效准则

目前对混凝土构件来说使用最多的就是以形变/能量双重准则为基础的损伤模型。可以通过损伤指数界定构件的抗弯失效。

使用OpenSees分析基础的受力得到计算轴力，再用Midas有限元程序对各桩柱截面的弯矩和曲率进行计算，各桩柱分析得到数据如表1所示。

表1 端部屈服转角和失效转角数值

失效参数如构件端部的屈服转角和极限转角则可以通过下面列出的公式进行计算得出。桩基础屈服转角按照式(6)计算：

(6)

式中：φy为截面屈服曲率；l为构件长度；0.002 5为构件抗弯屈服时的弯曲变形；εy=fy/Es为钢筋的屈服应变；fy为钢筋的屈服强度；Es为钢筋的弹性模量；asl为由于钢筋滑动而引起的端部转动；dh为纵筋的直径；d为截面直径；d′为核心混凝土的直径；fc′为混凝土的抗压强度。

构件端部的极限转角按照式(7)计算：

(7)

式中：φu为截面极限曲率；lp为塑性铰长度；φy为截面屈服曲率，那么φy可由式(8)、(9)计算得出：

(8)

(9)

式中：fy为钢筋的屈服强度；Es为钢筋弹性模量；fc′为混凝土的抗压强度；εc为混凝土的应变；Ec为混凝土的弹性模量；d为构件截面直径；ky为选取截面的屈服曲率最小值。

其中式(8)适用于以受拉钢筋屈服作为屈服标志的情况，式(9)适用于根据混凝土纤维在受压区的非线性应变是否大于εc=1.8fc′/Ec作为屈服标志的情况。

得出剪切失效判定准则：考虑轴力的影响，根据经典剪切-摩擦理论的计算公式，桩基础的抗剪失效可按式(10)计算：

(10)

群桩基础的剪切失效曲线参数取值和剪切失效曲线表达式取值见表2、3。

表2 群桩基础的剪切失效曲线的参数取值

表3 群桩基础剪切失效曲线表达式取值

注：V为桩端部剪力(kN)。

构件轴向失效准则：通过式(11)计算构件轴向上的失效情况：

(11)

其中：Δa为轴向极限状态的侧向位移；L为钢筋混凝土柱的高度；θ=65°为剪切破坏面的方向；P为结构在自重荷载下的轴力；st为构件中箍筋的间距；Ast为箍筋的横截面面积；fst为箍筋的屈服强度。式(11)中的参数取值见表4。

表4 式(11)中相关参数数值

2.4 工况选择

该文选择5个工况。工况1是不对过程干预的常规非线性的处理方法，工况2～5则对应于该文中论证的对程序进行控制的结构逐步失效的分析方法。工况2中释放首个破坏杆端的约束和弯矩，然后对结构响应进行分析。工况3的主要目的是对弯矩、剪力、轴力变化进行分析，要求对首个两端都破坏的杆件实施约束处理。工况4则在工况3之后，去除两端都发生塑性铰破坏的桩柱，然后对剩下的结构进行动力响应研究。工况5是根据轴力失效准则，对所有杆件的约束和弯矩进行释放，找出首个达到轴向失效标准的杆件。具体工况见表5。

表5 工况设置

3 连续倒塌分析结果

该文选择了3条地震波作为地震作用输入群桩体系模型，先对比说明3种地震波作用下群桩体系桥梁倒塌行为的异同点，再分别对结构模型的响应进行分析。3种地震波下倒塌行为异同点如下：

相同点：① 塑性铰出现顺序相同，都是边桩先发生屈服，在桩下端先出现塑性铰，然后桩上端屈服，产生塑性铰；② 在相同的工况下各桩杆端的位移、弯矩、剪力变化趋势大体相同。

不同点：① 由于地震波的随机性，在相同的工况下各个桩出现屈服、塑性铰的时间不同；② 在相同的工况下各个桩柱杆端抗弯、剪切、轴向失效顺序和时间会根据地震波的不同而发生变化。

3.1 Saratoga - Aloha Ave波输入下连续倒塌分析

3.1.1 结构动力时程响应

输入Saratoga - Aloha Ave波，处理表5中各种工况的结构响应，得到曲线见图6～10。

图6 输入地震波1：工况1和2时11桩底弯矩时程曲线

从图6可以看出：6.115 s时，工况1、2出现明显区别，工况1显示构件这时还有很高的承载力，而工况2显示构件的承载力变为0。此时11桩下端已经发生了弯曲失效，所以可以很明显看出工况1与构件损坏后端部抗弯承载力消失的真实情况不匹配，而工况2释放了失效构件端部的弯矩最后得到的结果也更符合实际情况。

图7 输入地震波1：工况1、2和5中的桩顶位移时程曲线

由图7可知：6.115 s前工况1、2和5中的桩顶位移变化趋势基本相同，位移量都不是很大，但是在6.115 s时，桩底达到弯曲失效标准，位移量开始变大，由于3种工况对杆件的控制不同，3条曲线开始出现区别，工况1变大后慢慢趋于平稳，工况2释放了失效杆端的内力，幅度比工况1稍有增大，而工况5由于对所有杆件都释放了内力，位移出现大幅度变化，基础大变形，导致桥梁倒塌。

通过图8中显示的曲线，再查看相关规范可知，工况5下的位移角数值比各国规范中规定的值(美国加州0.054 rad，新西兰0.049 rad，日本0.049 rad)都大得多，可以看出结构会发生倒塌。

图8 输入地震波1：工况1和5中的桩位移角时程曲线

图9 输入地震波1：工况1、3和4中的桩顶位移时程曲线

由图9可知：工况3中对首个两端都失效杆件释放杆端约束和内力，最后得出的位移响应曲线和工况4中拆除失效杆件得到的墩顶位移响应曲线变化形式基本相同。那么可以基于工况3下的位移响应曲线，再综合考虑轴力方面的因素，可以得出拆除构件后的响应，这可以成为一个新课题。

图10 输入地震波1：工况1、2和4时11桩上端弯矩时程曲线

由图10可知：在6.115 s后，桩上端会随着下端的弯矩释放发生较明显的变化，不仅增大弯矩数值而且缩短达到失效标准的时间。

综合上面的分析，如果在分析过程中对程序进行干预，对不同的工况下采用不同的处理，模型进行实时更新，将会增大桥梁基础在强震下的响应。

3.1.2 结构失效路径分析

对讨论的5个工况来说，在6.115 s也就是11桩下端出现弯曲失效之前，杆端曲率和杆端的塑性铰屈服时间、顺序都是一样的。而且出现屈服现象和塑性铰的顺序为先是桩下端，后是桩上端。

弯曲失效分析。图11(a)对应的是工况1中的构件弯曲失效状况，不对杆端进行约束处理，失效杆端依次为11桩底，22桩底，33桩底。图11(b)对应工况2中的弯曲失效情况，首先发生抗弯失效的是11桩底，接着依次是22桩底、33桩底、22桩顶、33桩顶、11桩顶。由图11(b)可以看出在分析过程中释放了11桩下端约束和弯矩，通过对11桩下端的处理，结构的内力重分布，使得剩下的杆端都达到抗弯失效标准，并且使剩下5个杆端的失效时间提前。图11(c)和图11(b)对比看出，工况2的杆端失效顺序可以运用到工况3设定的场景中去，它们的先后顺序是一致的，需要说明的是，工况3释放的杆端约束和弯矩比工况2多，所以工况3的杆端发生失效更早。图11(d)对应的是工况4的失效情况，从图中标的数字可以看出，首先还是11桩下端先失效，然后22桩和33桩底依次发生失效，第4个发生破坏的是22桩顶，对22桩进行去除处理，随后11桩和33桩顶一起失效，可以看出工况4的弯曲失效顺序和快慢跟前面都不一样。图11(e)对应工况5的失效情况，可以看出：虽然失效顺序和前面的工况2、3相同，但是全部桩端都进行了释放约束和弯矩处理，发生失效的时间会更早，同时也导致结构变形过大，发生倒塌。综上可知，对各个工况采用不一样的处理方法，弯曲失效顺序和快慢会改变。

图11 基础在工况中弯曲失效情况

剪切失效分析。剪切失效主要发生在工况4中，分析如下：22桩底在6.250 s时达到弯曲失效标准，释放22桩底约束和弯矩，其弯矩值马上从13 330.9 kN·m变到4 590 kN·m，慢慢地弯矩值变得接近于0，由于弯矩的变化，剪力随之开始突变；22桩上端在6.320 s时也达到抗弯失效标准，和上端弯矩值一样也开始发生瞬变，随着时间增加也变为0；桩两端剪力随着桩两端弯矩的变化开始再次突变，和弯矩值一样会变为0。构件的上下端部都失去抗弯能力后，剪力和弯矩值一样成为0，这时结构会发生剪切失效。

轴力失效分析。轴力失效主要发生在工况5中，分析如下：工况5环境下全部桩的端部都要释放约束和内力，随后桩也都会出现剪切消失，这时对桩柱接着研究，发现22桩会出现轴向失效，接着11桩也出现轴向失效，这时候结构已经不再具有足够的稳定性，桥梁因为变形过大发生倒塌。

3.2 El-Centro地震波输入下连续倒塌分析

3.2.1 结构响应

输入El-Centro波，研究典型的工况1、2和5的桩柱响应，整理出的时程曲线见图12、13。

图12 输入El-Centro波：工况1、2和5桩柱顶位移时程曲线

从图12看出：释放杆件的约束和弯矩后，结构的时程响应会发生变化，释放的弯矩越多变化就会越大，工况5中因为释放的弯矩过大，最后发生倒塌。

从图13可得出：2.00 s时，工况1和2出现明显区别，工况1显示桩柱这时还有很高的承载能力，而工况2显示桩柱的承载能力变为0。

图13 输入El-Centro波：工况1、2下11桩底弯矩时程曲线

3.2.2 结构失效路径分析

弯曲失效分析。工况1中不对杆端进行控制，11桩底，22桩底，33桩底依次达到弯曲失效标准。工况2中，首先发生抗弯失效的是11桩下端，接着依次是22桩底、33桩底、22桩顶、33桩顶、11桩顶。在分析过程中释放了11桩下端约束和弯矩，通过对11桩下端的处理，结构的内力重分布，使剩下的柱端都达到抗弯失效，并且使剩下5个桩端的失效时间发生了提前。工况5，虽然失效顺序和工况2、3相同，但是全部桩端都进行了释放约束和弯矩处理，发生失效的时间会更早，同时也导致结构变形过大，发生倒塌。由上可知，对各种工况采用不一样的处理，桩端失效的顺序和时间会发生改变。

3.3 Gilroy Array 4#波输入下结构连续倒塌分析

3.3.1 结构动态响应分析

输入Gilroy Array 4#波，分析工况1、2和4的结构响应，得到的时程曲线见图14、15。

图14 输入Gilroy Array 4#波：工况1、2和

从图14中可以得出：释放杆件的约束和弯矩后，结构的时程响应会发生变化，释放的弯矩越多变化就会越大，工况4中因为拆除了失效杆件，变形过大，最后发生倒塌。

从图15可知：在5.235 s时，工况1、2出现明显区别，工况1显示桩柱这时还有很高的承载能力，而工况2显示桩柱的承载能力已经变为0。工况2释放了失效构件端部的弯矩最后得到的结果更符合实际情况。

图15 输入Gilroy Array 4#波：工况1、2桩底弯矩时程曲线

3.3.2 结构失效路径分析

弯曲失效分析。工况1中不对杆端进行控制，11桩底，22桩底依次达到弯曲失效标准。工况2中的弯曲失效情况，首先发生抗弯失效的为11桩底，接着依次是22桩底、33桩底、22桩底。在分析过程中释放了11桩下端约束和弯矩，通过对11桩下端的处理，结构的内力重分布，使剩下的桩柱端都到达弯曲失效标准，并且使桩柱端的失效时间发生了提前。而工况4中因为拆除了失效杆件，变形过大，最后发生倒塌。

轴力失效分析。当22桩上下两端都发生弯曲失效，拆除22桩，继续对桩柱进行分析，33桩和11桩会发生轴向失效，这时结构失去承载力，因为变形过大，最后发生倾倒。

4 结论

该文采用OpenSees软件，把群桩基础模拟为带有附属节点的塑性铰梁柱单元，通过对附属节点进行处理和控制，实现了桥梁逐步失效的连续倒塌模拟功能，并应用于群桩体系桥梁连续倒塌模拟研究，得到结论如下：

(1) 根据该文采用的连续倒塌模拟算法，考虑弯曲、剪切和轴向失效准则，与常规方法比较可以更加高效地实现桥梁在地震作用下的连续倒塌模拟。

(2) 通过5种工况对桥梁结构失效响应过程进行比较，得到释放桩基单元端部的约束和内力，地震响应有变大的趋势，同时构件杆端的失效顺序和时间也会出现变化。

(3) 当桩的端部出现弯曲失效后，端部的剪力会突然变化，当构件的两个端部都抗弯失效后，桩柱会剪切失效。

(4) 在时程分析过程中对有限元模型进行实时处理并更新，会使结构的失效过程更加准确地反映出来。

(5) 杆件逐步失效过程中的内力释放会使结构的地震响应增大，同时它还影响桥梁整体损伤演变历程、发生倒塌的模式与强震作用下群桩基础连续倒塌机制。