带主梁的简化模型与响应面联合的桥梁船撞易损性分析方法

樊伟，孙洋，申东杰，刘斌

（1.湖南大学土木工程学院，湖南长沙 410082；2.湖南大学风工程与桥梁工程湖南省重点实验室，湖南长沙 410082）

船撞桥事故不但威胁着船舶通行的安全，也严重影响着桥梁的安全运营[1]，一旦发生船撞桥事故，桥梁结构可能承受巨大的侧向冲击荷载，因此在对通航水域的桥梁进行设计时，必须考虑船桥发生碰撞的可能性.

20世纪60年代末，人们开始研究船桥碰撞问题，著名的米诺尔斯基（Minorsky）理论[2]就在那时被提出，该理论为后人研究船桥碰撞奠定了基础.目前，桥梁船撞安全问题已经受到了广泛关注.相关研究主要集中于撞击力估算[3]、船撞桥数值模拟[4-6]、桥梁防撞设施研究[7-10]等.但是，以往的这些研究大多都是基于确定性的有限元仿真计算进行分析，难以反映船撞作用的偶然性和概率性特征，以及不同能量撞击下的桥梁损伤演化特征.对桥梁进行船撞易损性分析可以预测结构在不同类型的船舶以及不同船速撞击下发生各级破坏的概率，对结构的设计、加固和维修决策具有实际工程应用价值.

目前，桥梁的地震易损性研究已经受到国内外学者的广泛关注[11-15].Singhal等[16]采用贝叶斯原理，分析1994年1月17日北岭加州地震混凝土框架建筑的地震损伤数据，采用Park-Ang地震损伤指数表示结构损伤，将地震易损性定义为在给定的地震强度下，损伤指标超过某一定值的条件概率.目前，地震易损性的理论研究已经比较丰富.

与地震易损性研究相比，关于桥梁船撞易损性的研究成果却寥寥无几.首先，船撞领域尚未有人提出一个广泛令人接受和信服的判断桥墩船撞损伤等级的指标.张太磊[17]通过数值仿真计算，得出以墩底转角作为评价桥梁损伤等级的指标，并指出混凝土桥墩在受到船舶撞击后的破坏形式与地震作用下的破坏形式有所不同.但是其仅仅对矩形墩柱在船舶正向撞击方面进行了研究，因此其损伤指标的适用性还有待进一步研究.

此外，前期研究[18]表明，桥梁船撞过程动力效应影响明显，对船撞下桥梁结构的影响进行合理的动力分析是必要的，因此，进行桥梁结构的易损性分析时，就需要大量样本的动力计算.如果采用常规的精细化接触-碰撞有限元技术进行船桥碰撞的非线性显示动力分析，单个模型的计算将会消耗大量的计算时间，计算效率低.若要进行大量样本计算，精细化全尺寸有限元模型必然会制约运算规模，因此，提高计算效率非常关键.近年来也有学者[19-20]将响应面运用到了撞击下桥梁参数分析和可靠度分析中，极大地提高了计算效率.其中就有张军等[19]提出运用响应面法进行船撞桥的可靠度计算，但是其并未进行非线性有限元动力碰撞计算，而是仅仅以规范的计算方法确定船撞力.

由此可见，高精度高效率的有限元简化模型与响应面代理模型的联合是实现大样本分析的前提条件.因此，本文提出有限元简化模型与响应面代理模型联合的桥梁船撞易损性分析方法，为构建基于概率性的桥梁船撞设计与评估方法奠定基础.

1 墩柱剩余承载能力试验与直接模拟方法

1.1 墩柱剩余承载能力试验

Fan和Liu等[21-22]对一组不同参数的受压RC墩柱进行了落锤冲击试验，并对受到落锤冲击后的受压RC墩柱进行了轴向加载试验.试验共对10根不同参数的RC墩柱进行了轴向加载试验，其中有2根未受损伤的试件，还有8根受落锤冲击后的受损试件.这8根受损试件中，有1根被完全砸断，还有2根出现了“反拱”现象[21].Fan等[21]指出，反拱的现象具有偶然性，存在一定的不可重复性.因此，在本次模型验证中排除了这2根“反拱”试件以及1根完全砸断的试件，对剩余的7根试件开展了数值模拟.

1.2 墩柱剩余承载能力直接模拟方法

为了研究墩柱受冲击后剩余承载能力的直接模拟方法，采用LS-DYNA显式非线性动力有限元软件建立了如图1所示的精细化有限元模型.

图1 有限元模型计算过程图Fig.1 Calculation process of FE model

其中，纵向和螺旋钢筋均采用Hughes-Liu积分梁单元模拟，采用*MAT_PIECEWISE_LINEAR_PLASTICITY材料本构，该材料本构可以采用试验测得的实际应力-应变关系来定义.混凝土采用六面体单点积分实体单元模拟，材料模型采用连续盖帽本构模型*MAT_CSCM，该材料本构已被证明可以合理地模拟冲击荷载下受压RC柱的动力响应[22-23].钢筋的梁单元与混凝土的实体单元之间使用了

*MAT_GENERAL_NONLINEAR_1DOF_DIS CRERE_BEAM非线性弹簧单元相连接，来模拟钢筋与混凝土之间的粘结滑移关系.非线性弹簧单元的加载卸载曲线的定义参照欧洲规范[24]中的规定.

为了模拟出落锤冲击后进行轴向加载的试验全过程，在建模过程中将三个试验阶段合并在一个计算模型中，如图1所示.整个计算过程分为三个阶段，计算时间为0.14 s.

第一阶段为初始轴力的预加载，时间为0 s～0.022 s.模型中采用关键字*CONSTRAINED_NODAL_RIGID_BODY将预应力钢筋的两端与试件两端的钢板固定起来，预应力钢筋由索单元模拟，通过定义材料*MAT_CABLE_DISCRETE_BEAM将初始应力赋给了预应力钢筋；第二阶段为落锤冲击阶段，时间为0.022 s～0.08 s.在初始轴力完成加载后，落锤将以给定的初速度下落并与受压RC柱发生碰撞，在早于0.08 s的某个时间点完成全部碰撞过程，撞击力归零，锤头与试件完全分离；第三阶段为轴向加载阶段，时间为0.08 s～0.14 s.通过关键字*BOUNDARY_PRESCRIBED_MOTION_SET控制墩柱的两端钢板，在0.08 s时采用预先设定好的轴向位移开始压缩墩柱，即位移控制的轴向加载.加载过程持续到墩柱失去轴向承载能力，结构失效.

1.3 模型验证

试件剩余承载能力的有限元结果与试验结果的比值如图2所示.

图2 剩余承载能力对比Fig.2 Comparison of residual bearing capacity

所有工况的比值范围为75%～126%，平均比值约为93%.总的来说，有限元模型对剩余承载力的预测能力达到了合理的精度，说明这种墩柱剩余承载能力的分析方法可用于后续的研究分析中.

2 带主梁的桥梁船撞简化分析模型

2.1 典型连续梁桥的有限元模型

本文以一座典型的四跨连续钢筋混凝土梁桥为例进行具体分析，如图3所示.桥墩及截面布置如图4所示.桥梁跨径为50 m，桥墩形式为双柱式桥墩.被撞桥墩为正中间的二号桥墩.单个桥墩尺寸为2.2 m×1.9 m×18 m.桥墩混凝土抗压强度为30 MPa.桥墩沿长度方向布置22根纵筋，沿宽度方向布置19根纵筋，间距均为100 mm.纵向钢筋直径为32 mm，屈服强度为400 MPa，纵筋配筋率为1.5%.箍筋间距为100 mm，直径为16 mm，屈服强度为335 MPa，箍筋体积配箍率为0.41%.桥墩混凝土保护层厚度为50 mm.

图3 全桥示意图（单位：m）Fig.3 Diagram of the whole bridge（unit：m）

图4 桥墩及截面布置图（单位：m）Fig.4 Diagram of single pier（unit：m）

为了研究桥梁结构抵抗船舶撞击的能力，采用上述模拟类似方法对该桥例建立了船桥碰撞的精细化有限元模型，如图5（a）所示.其中，被撞桥墩混凝土采用经常被用于低速冲击模拟的连续盖帽本构*MAT_CSCM（159）进行模拟.钢筋使用具有应变率效应的理想双折线弹塑性材料本构*MAT_PLASTIC_KINEMATIC（003）模拟.钢筋与混凝土之间考虑了粘结-滑移的影响，处理方式与1.3节模型验证相同.考虑到非被撞桥墩的响应相对较小，变形均在弹性范围内，为节省计算时间，其余非被撞桥墩均采用线性弹性材料*MAT_ELASTIC（001）进行模拟.本模型中所有盖梁、承台和桥墩均采用六面体单点积分实体单元；桩基和钢筋均采用Hughes-Liu积分梁单元；下部结构采用梁单元桩基加土弹簧单元进行模拟；主梁采用Belytschko-Schwer resultant beam梁单元，该梁单元可以自由定义截面特性.此外，本模型采用关键字*LOAD_BODY对全桥施加了重力，以考虑桥墩的初始轴压状态.船模采用的是前期已经进行过标定的精细化有限元船模[7]，如图6所示.

图5 精细化船桥碰撞有限元模型及其简化模型Fig.5 Full-size and simplified ship-bridge collision FE model

图6 精细化有限元船模Fig.6 Refined vessel FE model

2.2 简化模型

精细化全尺寸桥梁有限元模型的计算需要消耗大量的计算时间.尤其是在工况数量较多的情况下，简化分析模型在大样本的参数分析或概率分析中显得尤为重要.

本小节使用Fan等[25]提出的带主梁的桥梁船撞简化分析模型，如图5（b）所示.通过将原有支座底端的竖向自由度约束住，使用横桥向的弹簧代替原有的非被撞桥墩的侧向刚度.使用纤维梁单元建立单个桥墩的OpenSees模型，在盖梁顶部支座约束处施加横桥向单位位移，并提取墩顶反力，则可以得到弹簧弹性刚度为1.442×108N/m.

2.3 模型验证与讨论

本次模型验证采用2000DWT驳船以2 m/s的船速撞击桥墩，如图5所示.对比上述简化模型与全尺寸模型的船桥碰撞响应：包括撞击力、被撞点位移、各墩的墩顶位移及被撞桥墩的墩底的剪力，计算结果如图7所示.

使用简化模型可以极大地提高运行效率.使用8cpu工作站计算时，带主梁的简化模型耗时约17 h，而全尺寸模型耗时约50 h.

如图7所示，将本文所提出的带主梁的简化模型的各项计算结果与全尺寸模型对比可知，撞击力与各墩墩顶位移均能与整桥模型达到极高的吻合度，表明该简化模型可以很好地代替整桥模型，既能提高运行效率，又能高精度地反映整桥模型墩柱撞击过程中的动力响应.

图7 简化模型与全桥模型对比Fig.7 Comparison between simplified model and full-size model

因此，后面易损性分析部分均使用带主梁的简化模型进行计算.

3 连续梁桥船撞易损性分析

3.1 桥梁船撞易损性分析方法框架

结构的地震易损性是指在给定强度的地震作用下，结构达到预定极限状态的条件失效概率，基于概率理论对结构的抗震性能做出定量评价，描述输入地震动强度与建筑结构损伤程度之间的关系.本文借鉴结构的地震易损性定义，可将桥墩的船撞易损性定义为：当发生冲击强度A=a的船撞冲击时，结构达到或超过某种极限状态（Ls）的条件失效概率：

改变冲击强度a的数值，计算结构达到或超过破坏状态Ls的船撞失效概率FR，然后采用某种统计方法进行曲线拟合，所得的光滑曲线FR（a）称为“船撞易损性曲线”.

在本文所研究的桥梁船撞易损性分析中，取参数A为船舶的冲击速度，即定义船撞易损性为当发生船速V=v的船撞冲击时，结构达到或超过某种极限状态（Ls）的条件失效概率为：

因此，在进行桥梁船撞易损性分析之前，首先需要根据损伤指标D的范围给桥墩损伤进行分级，从而确定不同的极限状态（Ls）：轻微损伤，中等损伤，严重损伤，倒塌.

目前在船撞领域，对于桥墩受船舶撞击之后损伤情况的界定，没有像地震分析中存在一个统一的令人信服的评估指标[26].考虑到桥墩的主要功能是支撑上部结构并将上部结构的恒、活载传递至基础，因此本文认为将桥墩受船舶撞击后的剩余轴向承载能力作为评估桥墩损伤的指标是合理的.文献[27]给出了基于墩柱残余承载力的损伤指标分级，如式（3）（4）所示：

式中：Pr为桥墩受到损伤后的残余轴向承载能力；Pd为桥墩的初始轴向承载能力.

尽管采用简化分析模型能够提高计算效率，但需要大量计算的易损性分析而言，仍然是不够的.因此，本文将用简化分析模型获得关键样本结果，采用响应面方法建立合理的响应面替代模型.在此基础上，结合蒙特卡洛抽样的方法计算获得不同船撞情况下各等级损伤下的失效概率.

3.2 试验设计

Box-Behnken试验设计法（Box-Behnken Design，简称BBD），是响应面优化法常用的试验设计方法，以三因素为例，试验设计的试验点分布情况如图8所示.

图8 三因素BBD试验点分布情况Fig.8 Sample designs of 3-Factor，3-Level Box-Behnken Design

合理的试验方法是获得良好的响应面代理模型的前提，因此，在进行响应面工况设计之前，对桥墩在分别受到驳船和球艏船不同船速撞击下的结构响应进行了数值模拟，结果如图9所示.

图9 不同船速下的桥墩响应Fig.9 Response of piers under different vessels speeds

可以看出，桥墩在两种不同类型船撞击下的响应特征截然不同.在球艏船的撞击下，桥墩的响应随着船速的增大近似线性上升，而在驳船的撞击下，桥墩的各响应呈现出双折线的形式.当驳船的船速在0～0.8 m/s内上升时，桥墩的峰值响应迅速上升；当船速超过0.8 m/s后，峰值响应基本上不会再有很大变化.产生这种现象的原因可能是，当撞击速度超过0.8 m/s后，驳船的船首在发生船桥碰撞的过程中前肋板屈曲，船首刚度突然下降至低于桥墩的侧向抗推刚度，驳船剩余的动能几乎全部由船首吸收，因此，桥墩响应未见上升而驳船船首变形严重.这也是当船速超过0.8 m/s后，撞击力-时程曲线出现平台段的原因，如图10所示.

在上述工况分析中，可以认为0.8 m/s是使船首屈曲的临界船速v0.根据试算可知，临界船速的大小并不是一个固定的数值，它与船舶总质量、混凝土强度、钢筋直径等因素息息相关.

图10 桥墩受驳船不同船速撞击下的力-时程曲线Fig.10 The impact force under different vessels speeds

考虑到桥墩在受到驳船和球艏船撞击下的不同响应特征，如果继续使用传统的Box-Behnken设计法，仅仅在每种因素（混凝土强度、钢筋直径、船速）的最大值、最小值和中心点进行试验，那么船速这一因素对驳船撞击下桥墩各响应的影响将会被错误地估计，如图11所示.由这样的试验数据点拟合出来的响应面代理模型并不能反映桥墩受驳船撞击后响应的真实特征，代理模型外推能力较差.

对于船速这一因素来说，临界速度v0是一个重要的数据点.在进行工况设计时，必须引入临界速度的影响.因此，本文针对驳船撞击下的工况，提出了一种基于临界速度的分段BBD试验设计法，将原本[0，4]的船速区间分为[0，v0]和[v0，4]两段.对每一个区间分别进行BBD试验设计，如图11所示.

图11 两种工况设计方法对比Fig.11 Comparison of two design methods

当桥梁受到船撞时，响应结果往往受到各种参数的随机性的影响，影响桥墩受冲击后剩余承载能力的主要参数有混凝土强度fc、纵向钢筋直径dl、箍筋直径ds、纵向钢筋屈服强度fy、箍筋屈服强度fys和船速v等.根据试算可知，ds、fy、fys在其服从的随机分布范围内的变化对剩余承载能力的影响较小，故在本次研究中并未考虑这3个参数的随机性带来的影响，并在响应面代理模型的计算中将这3个参数设为固定值（ds=16 mm，fy=400 MPa，fys=335 MPa）.为了对桥梁船撞易损性进行分析，本节将建立如下响应面代理模型来预测桥梁受船撞后的残余承载能力：

参考现有文献的结论[28-30]，fc和dl均服从正态分布，统计参数如表1所示.

表1 桥墩随机变量所服从的统计参数Tab.1 Statistical parameters of random variables

参考两种参数的概率密度分布情况，在进行试验设计时，将钢筋直径的变化范围设为28～36 mm，将混凝土强度的变化范围设为20～40 MPa.设计工况如表2所示.

表2 不同船舶撞击下的BBD设计工况Tab.2 BBD table for the case of different vessels

3.3 建立响应面代理模型

3.3.1 球艏船撞击下的桥墩响应代理模型

多项式响应面代理模型是响应面分析中常用的代理模型形式，本文首先尝试采用多项式代理模型对桥墩受驳船以及球艏船撞击下的剩余承载能力进行拟合.

代理模型采用多项式形式，其函数形式为：

式中：fc为混凝土强度；dl为钢筋直径；v为船速；p1～p13为多项式中各分项的待定系数.

根据试验结果拟合得到桥墩受球艏船撞击后的残余承载力响应面代理模型多项式的各项系数如表3所示.二次多项式代理模型即可达到令人满意的拟合精度，可决系数（R2）为0.98.各个参数对船撞后的残余承载能力的影响如图12所示，图中坐标轴上未显示的参数均取中值：v=2 m/s，fc=30 MPa，dl=32 mm.

表3 桥墩剩余承载能力响应面多项式各分项系数Tab.3 Coefficients of polynomials

图12 球艏船撞击下桥墩残余承载力响应面模型Fig.12 Response surface model of residual capacity under collision of bulbous-bow ship

然而，在拟合桥墩受驳船撞击下的残余承载能力响应面时，多项式代理模型始终无法达到一个令人满意的拟合精度，无论是二次多项式还是三次或者更高次多项式，最后的效果始终差强人意.

3.3.2 驳船撞击下的桥墩响应代理模型

为了提高响应面代理模型的拟合精度，针对桥墩在受到驳船撞击下的残余承载能力，本文提出了一种基于临界船速的多项式分段函数的代理模型，如式（7）所示：

根据试验结果拟合得到桥墩受驳船撞击后的残余承载力响应面代理模型多项式的各分项系数如表3所示.分段拟合后，三次多项式代理模型可以达到令人满意的拟合精度，分段函数f1和f2的可决系数（R2）分别为0.999 9，0.999 3.

桥墩在受驳船撞击后的残余承载能力与各参数之间的关系如图13所示.当船速小于0.8 m/s或大于0.8 m/s时，桥墩残余承载能力的响应特征有着明显差异.

由图13（c）中船速和混凝土强度对桥墩残余承载能力的影响可明显看出，当混凝土强度下降时，临界船速也会变小，具体表现为随着混凝土强度降低，响应面的转折点向船速低的方向移动.

图13 驳船撞击下桥墩残余承载力响应面模型Fig.13 Response surface model of residual capacity under collision of barge

值得一提的是驳船撞击下的桥墩剩余承载能力，如图13（b）（c）所示，当船速超过临界船速并继续提高至约2.5 m/s时，桥墩的剩余承载能力继续平缓下降；但是，当船速超过2.5 m/s时，桥墩的剩余承载能力反而呈现出上升的趋势，这一现象与荷载频谱特性、桥梁结构动力特性相关[25].

3.4 易损性分析

根据墩柱残余承载力的损伤指标，即式（4），进行桥墩船撞的易损性分析.其中Pr可由相应的响应面代理模型计算得到.将fc=30 MPa、dl=32 mm、v=0 m/s代入相应的响应面代理模型，则可得到Pd=f（30，32，0）.

对于每一个确定的船速v∈（0，4）：对服从正态分布的随机变量fc和dl进行蒙特卡洛抽样，样本数量N=107，统计出损伤指标D分别达到轻微损伤、中等损伤、严重损伤和倒塌的次数L1、L2、L3、L4.则在该船速下，桥墩达到相应损伤指标的失效概率为：

式中：n=1、2、3、4分别表示损伤程度达到轻微损伤、中等损伤、严重损伤和倒塌.对船速v从0 m/s到4 m/s，船速每间隔0.02 m/s计算一次结构的失效概率，将得到的数据点连线则可得到桥墩受船撞的易损性曲线，如图14所示.

对比驳船与球艏船的易损性曲线，可知在相同船速及质量的情况下，驳船撞击所引起结构失效的概率普遍比球艏船要大很多.相比较于球艏船撞击工况，驳船撞击下的桥墩在较低的船速撞击下就能造成较大的损伤.

图14 桥梁船撞易损性曲线Fig.14 Vulnerability curve of vessel collision

4 结论与展望

本文面向两类典型船舶，以桥梁墩柱受船舶撞击后的剩余承载能力作为损伤评估指标，提出了基于带主梁简化分析模型与响应面联合的桥梁船撞易损性分析方法，主要结论有：

1）所建立的高精度响应面模型可以替代需要进行复杂非线性计算的结构模型.可以使用所建立的响应面代理模型在桥梁船撞易损性研究中进行大量的样本分析.

2）由于船首构造、外形等差异，不同船舶类型撞击下的桥墩剩余承载力的响应特征区别较大.在球艏船撞击下，桥墩的剩余承载能力会一直随船速的增大而均匀减少；而在受驳船撞击时，桥墩的剩余承载能力与临界船速密切相关，呈现出双折线的特征.因此，在进行样本设计时，为了反映桥墩受驳船撞击后响应的真实特征，需基于临界速度进行分段.

3）在相同船速及质量的情况下，驳船撞击所造成的结构损伤以及失效的概率普遍要高于球艏船撞击.在进行实际风险评估时，应该尤为关注.

4）在本文算例分析中，无论是驳船还是球艏船都未能对桥墩造成严重及以上的损伤状态，但是实际工程中却屡见船桥碰撞的恶性事故，这是因为实际工程中的桥梁船撞事故往往还伴随着钢筋锈蚀以及下部结构冲刷等多种灾害同时作用.此外桩土、水位变化、撞击角度、通航等级与船只吨位等各种不确定因素都会对事故的结果造成影响，这些都是值得进一步研究的问题.在船撞易损性分析中考虑钢筋锈蚀以及冲刷等因素的影响将是本文的下一步研究方向.