实测模态下三跨整体桥精细化基准有限元模型

张志勇，何乐洽，薛俊青，陈梦圆，黄福云，Giuseppe Carlo Marano，Bruno Briseghella

(福州大学土木工程学院，福建 福州 350108)

0 引言

基于物理模型的损伤识别，是土木工程结构健康监测的一项重要研究内容. 通过比较结构无损状态和损伤状态下动力特性的差异，可发现定位和定量结构的损伤程度，从而为工程结构的智能化管养提供有力的决策支持[1]. 有限元技术是现代桥梁与建筑结构设计和分析的重要支撑. 通过建立合理的设计模型，采用必要的设计假定和简化，能够解决结构设计阶段所面临的主要问题[2]. 大量的结构长期健康监测结果表明，工作状态下的梁桥，其动力特性(如自振频率、阻尼等)受到各种因素的显著影响[3]. 前人的研究工作揭示，根据结构无损/初始状态下实测模态信息建立的有限元基准模型可为结构的长期健康监测提供精确、可靠的参考信息[4]. 作为能够全面、正确反映结构真实行为的基准有限元模型，除了应保证与实测结果较好的吻合度外，还应保证模型的各项参数真实准确的物理意义.

整体式桥台桥梁(简称整体桥)属于无伸缩缝桥梁，常见于中小跨径桥梁[5-6]. 整体桥将传统桥梁的上下部结构连成一个整体，为高阶超静定结构. 本文以一座三跨整体桥为工程背景，考虑结构-土相互作用，基于环境振动试验实测模态信息，开展了该桥基于梁格法设计简化模型与梁-壳单元精细化模型的模型修正. 将与试验频率和振型的吻合度为目标，通过调整模型的物理参数，建立该桥初始状态的基准动力有限元模型[7]. 最后为整体桥动力分析提出了建议.

1 工程背景

马峦山桥位于深圳市马峦山郊野公园北侧，于2019年年底建成通车(见图1). 该桥分南、北两幅，双向八车道. 由于两幅桥的结构完全一致，故仅以南幅桥作为本文的研究对象. 该桥是一座三跨预应力混凝土箱梁整体桥，采用先成桥后填筑台后土的方式建造. 桥台后填土采用中等密实砂性土. 在桩基深度方向，桩周土分别为素填土、粉质粘土、粘土、残积土与不同风化程度的粉砂岩. 桥梁跨度为87.92 m(28.96 m+30 m+28.96 m)，桥宽为17 m，主梁采用单箱三室等高箱梁. 全桥取消了桥台和桥墩处的伸缩缝、伸缩装置和支座，采用整体式桥台将主梁梁端与桥台固结在一起，墩梁也采用固结. 墩下和桥台下均采用单排混凝土桩基础. 由于温度变形产生的梁体伸缩量引起桥台和桥墩的纵桥向变形，并同时通过柔性桩和桥头搭板引起台后填土和桩周土的相应变形，从而产生结构和土的相互作用.

图1 深圳马峦山整体桥Fig.1 Maluan mountain integral abutment bridge in Shenzhen

2 环境激励下的模态试验

图2 试验测点布置图(单位： cm)Fig.2 Layout of test points (unit: cm)

为了获取该桥实际的动力特性，在风和少量车辆为主的环境激励下进行不同测点布置的两次模态试验. 试验测试仪器为4台瑞士GEOSIG公司生产的三向无线加速度传感器GMSPLUS，可同时记录纵桥向、横桥向和竖向的加速度响应. 试验一与试验二中，分别在桥面上等间距的布置了13个与39个测点. 测点位于道路两侧和主梁横断面中间，测点的布置如图2. 采用随机子空间算法，基于实测加速度信号识别到了该桥前五阶试验模态，包括一阶竖弯、二阶竖弯、三阶竖弯、一阶扭转和二阶扭转. 试验结果详见表1，其中f为自振频率，ξ为阻尼比.

表1 实测模态与修正后计算模态对比

3 有限元模型的建立

3.1 基于梁格法的有限元模型

3.1.1上部结构模拟

图3 箱梁截面图与梁格法划分(单位： cm) Fig.3 Cross section of box girder and division of grillage method (unit: cm)

图3为马峦山整体桥的箱型主梁横截面. 该桥主梁宽17 m，高1.6 m，纵向预应力钢束布置在腹板内部. 主梁基于梁格法设计，采用MIDAS-civil软件建立全桥三维有限元模型. 图3中虚线所示为横断面梁格划分. 虚拟横梁沿纵桥向的间距为1 m. 虚拟横梁截面仅计入主梁的顶板和底板对于刚度的贡献，不计入其质量. 桥面混凝土调平层与沥青铺装层以及防撞护栏仅作为附加质量考虑，忽略其刚度影响. 模型中采用的是基于Timoshenko梁理论的三维梁单元，采用分布质量且不考虑钢筋的影响.

3.1.2下部结构模拟

桥台采用梁格法模拟，主梁与桥台采用共节点的方式连接，桥墩、承台和桩基均用梁单元模拟，相互之间采用主从节点的刚性连接.

模型中考虑了桥台-台后土与桩-桩周土的相互作用，分别采用非线性与线性拉压弹簧模拟桥台-土与桩-土的相互作用. 其中桥台处仅在纵桥向设置了土弹簧，而桩周在纵、横桥两个方向上均设置了相同的土弹簧，如图4所示. 下面简要介绍桥台土弹簧与桩周土弹簧刚度的计算方法[5, 8].

该桥台后土为中等密实砂性土，桥台与台后土仅考虑纵桥向的相互作用，忽略桥台与土之间的竖向和水平向摩擦力. 模型中，将台后土弹簧作用在桥台梁格节点上，按照该节点与相邻节点的间距和土层的性质来确定其刚度. 例如对于桥台上某节点，当桥台位移为零时受到静止土压力Po=Koγzbt的作用.模型中，将其作为附加力单独施加到桥台上； 当桥台挤压填土时，土压力变化至被动土压力Pp=Kpγzbt，此时桥台土弹簧受压； 当桥台远离填土时，土压力变化至主动土压力Pa=Kaγzbt，此时桥台土弹簧受拉.其中Ko、Kp、Ka分别为静止、被动、主动土压力系数，γ为土的容重，z为该节点处台后填土深度，b为该节点处左右相邻单元长度和的1/2，t为该节点处上下相邻单元长度和的1/2.由此，该节点处土弹簧受压时的刚度为k1=(Pp-Po)/Δ1； 受拉时的刚度为k2=(Pa-Po)/Δ2.其中，Δ1与Δ2为桥台节点位移. 上述各参数的取值均可参考NCHRP规范曲线的相应数值[9]. 相应地，桥台梁格模型上各节点的土压力与桥台位移关系可用图5来表示.

图4 马峦山整体桥结构模型南立面图(半幅，单位： cm)Fig.4 Structural model of Maluan mountain bridge: south elevation view(half span, unit: cm)

图5 桥台各节点土压力与桥台位移关系(H为桥台高度)Fig.5 Relationship of earth pressure at each node of abutment and its displacement(H: abutment height)

整体桥中桩-土的相互作用不可忽略，通常可按照温克尔假定，把桩视作弹性地基上的梁，将桩侧土离散为线性弹簧. 该桥桩周土包括素填土、粉质粘土、粘土、残积土与不同风化程度的粉砂岩(从上到下). 模型中，桩基梁单元按土层厚度进行取整划分. 在桩基节点处分别设置横桥向与纵桥向的土弹簧. 对桩基某节点处的土弹簧刚度，取为Kz=Cz·az·bz.其中，地基系数Cz=m·z，m为地基土的比例系数，z为桩基节点所在的地基土深度，az为该节点处土层的厚度，bz为该节点处桩基计算宽度. 上述各参数的取值均可参考《公路桥涵地基与基础设计规范(JTG 3363—2019)》[10]的附录P. 桩周土弹簧施加在桩基单元的节点处. 考虑到该桩为摩擦桩，且环境激励条件下竖向位移微小，与结构位移比较可以忽略不计. 故将桩基单元节点竖向固结.

3.2 基于梁-壳单元的精细化有限元模型

3.2.1上部结构模拟

利用通用有限元软件ANSYS建立了全桥精细化的三维梁-壳单元模型. 桥梁上部结构采用了壳单元SHELL181，将箱型主梁横截面分上顶板、中腹板、下底板三个部分进行模拟. 单元厚度分别取为相应截面厚度的平均值. 主梁在与桥台、桥墩连接处设置有横隔，横隔单元厚度按实际厚度取值.

将桥面混凝土调平层与沥青铺装层作为附加质量加到主梁顶板上，忽略其刚度贡献. 采用梁单元BEAM188模拟防撞护栏，考虑到护栏结构的非连续性和不确定性，根据参考文献[11]将其初始弹性模量做一定量的折减. 与梁格模型中一样，所有单元均采用了分布质量且不考虑钢筋的影响.

3.2.2下部结构模拟

桥梁下部结构中，桥台采用壳单元SHELL181模拟，桥墩、承台和桩基均采用梁单元BEAM188模拟. 主梁与桥台采用共节点的方式连接. 主梁与桥墩基于约束方程建立连接. 桥台与承台、桥墩与承台、承台与桩基均采用主从节点的刚性连接.

与梁格模型一样，梁-壳单元模型基于土弹簧考虑了桥台-台后土与桩-桩周土的相互作用. 用离散的非线性弹簧单元COMBIN39，模拟台后填土对桥台的作用. 同前，弹簧的刚度按NCHRP所给出的土压力系数曲线推算确定，具体计算方法和梁格模型相同. 同时，用离散的线性弹簧单元COMBIN14模拟桩周土对桩的作用，其刚度计算方法同前. 与梁格模型相同，将桩基单元节点采用竖向固结.

最终，简化的梁格模型与精细化的梁-壳单元模型分别见图6与图7. 其中，梁格模型总计有1 356个节点，1 918个单元，而梁-壳单元模型总计有6 662个节点，6 975个单元. 两类模型方位与图1实桥一致. 两侧桥台不等高，左侧为6.97 m，右侧为6.46 m，且台下与墩下的桩基长度也不等. 梁格模型与梁-壳单元模型各部分的材料特性取值参考了《公路钢筋混凝土及预应力混凝土桥涵设计规范(JTG 3362—2018)》[12].

图6 马峦山整体桥有限元梁格模型Fig.6 Finite element grillage model of Maluanshan bridge

图7 精细化的ANSYS梁-壳单元模型Fig.7 Refined ANSYS beam-shell element model

4有限元模型的参数化修正

4.1 目标函数的确定

采用基于灵敏度有限元模型修正，将有限元模型的计算模态与结构实测模态的频率和振型的残差平方和作为目标函数. 通过最小化该目标函数，可以获得较为精确的动力有限元模型. 该过程可视作一个数值优化求解问题，即如下所示非线性最小二乘问题：

(1)

式中：J(θ)为目标函数，其中θ为模型参数；Wε为权重因子对角矩阵用来平衡目标函数中各项残差的贡献.出于简化考虑，此处Wε取为单位矩阵；εz为模态数据z的计算值与实测值之间的残差.对于无阻尼特征值z=λ，此时残差ελ, i(θ)可表达为

(2)

式中：λi=ωi2=(2πfi)2，i∈{1, 2, …,nλ}.其中上波浪号表示实测值，f或ω表示结构固有频率或圆频率.

对于特征向量(振型)，采用模态置信准则(modal assurance criterion, MAC)计算，公式为：

其中:φa表示计算振型向量;φe表示实测振型向量.此时z=MAC，残差εMAC, i(θ)可表达为

εMAC, i(θ)=1-MAC

(3)

式中：i∈{1, 2, …,nMAC}.

对于梁格模型采用仅考虑频率的目标函数：

(4)

对于精细化的梁-壳单元模型，采用联合频率与振型的目标函数，最终的目标函数如下：

(5)

4.2 参数灵敏度分析与选取

为了选取对计算模态影响较大的模型参数，从而确定优化问题的设计变量，需要进行模型的参数灵敏度分析. 针对计算频率ω参数灵敏度矩阵定义为式(6)，则针对计算振型与实测振型MAC的参数灵敏度矩阵同理可得.

(6)

由无阻尼广义的特征值方程Kφi=λiMφi可知，同时修正计算模型的质量与刚度将无法得到唯一解. 故灵敏度分析时，不考虑模型参数的质量密度，仅分析模型参数的弹性模量. 针对表1中前五阶计算频率与振型MAC值，并考虑到一阶纵向刚体模态对整体桥的特殊性(结构-土相互作用)，分别对梁格模型与精细化梁-壳单元模型进行参数灵敏度分析.

分析结果发现： 在梁格模型中，自振频率对主梁与桥台梁格的纵梁与横梁的刚度均较为敏感，而对桩周与台后土弹簧刚度不敏感. 据此，考虑将主梁纵梁弹性模量、主梁与桥台虚拟横梁弹性模量、桥台纵梁弹性模量3个参数作为梁格模型待修正参数.

在梁-壳单元模型中，针对主梁弹性模量、桥墩弹性模量、桥台弹性模量、横隔梁弹性模量、护栏弹性模量、桩周土弹簧刚度、台后土弹簧刚度7个参数进行灵敏度分析，结果(取绝对值)见表2. 由表2可知，实测的前五阶频率对主梁刚度、桥墩刚度、桥台刚度、横隔梁刚度较为敏感，对桩周与台后土弹簧刚度灵敏度较小，且均小于2%. 桩周与台后土弹簧刚度仅对一阶纵向刚体频率敏感，约10%. 但由于在环境振动条件下未实测到该桥纵向刚体模态，不宜对其进行修正. 同时，由表1中各阶实测振型图结果获知，桥台底部振型各向位移和桥墩底部振型各向位移基本上均在0～0.03之间. 据此认为，在试验条件环境激励情况下，桥台和桥墩承台下的桩基基础刚度对结构整体刚度影响较小. 由于本文关注重点在于整体桥桥面的动力特性，故不再将台后土弹簧刚度与桩周土弹簧刚度作为修正参数. 且文中已将桥台弹性模量和桥墩弹性模量列为模型修正参数，认为可以部分补偿将桩基模型节点竖向固结造成的误差. 综合以上分析，同时考虑到非结构构件的影响，增加护栏刚度作为待修正参数. 最终确定以下5个参数作为梁-壳单元模型待修正参数，分别是： ① 主梁弹性模量； ② 桥墩弹性模量； ③ 桥台弹性模量； ④ 主梁横隔弹性模量； ⑤ 护栏弹性模量.

表2 梁-壳单元模型参数灵敏度矩阵

4.3 修正过程及结果4.3.1 梁格模型修正

参考式(4)定义的目标函数，采用手动调整4.2节中确定的3个模型参数的方式，进行梁格模型的修正. 修正前后梁格模型的计算模态与实测模态的对比见表3. 修正前后各参数的调整情况见表4. 由表3可知修正后的模型各阶频率的计算值与实测值绝对误差均在5%以内. 这表明选取的参数与修正的过程是基本合理的.

表3 修正前后梁格模型模态与试验模态的对比

表4 梁格模型修正前后参数变化

4.3.2梁-壳单元模型的修正

图8 梁-壳单元模型目标函数的迭代收敛曲线Fig.8 Iterative convergence curve of objective function in beam shell element model

根据式(5)定义的目标函数，利用MATLAB和ANSYS的交互式访问，实现了梁-壳单元模型的参数化修正. 其中，式(1)的求解采用了MATLAB优化工具箱提供的非线性最小二乘函数“lsqnonlin”. 该函数采用了基于Trust-region-reflective(TRR)的优化算法.

表5 修正前后梁-壳单元模型模态与试验模态对比

表6 梁-壳单元模型修正前后参数变化

上述梁格模型与梁-壳单元模型的修正均以试验一的模态结果为参考. 实测模态与修正后的计算模态对比见表1.

4.4 修正结果分析

由表3和表5可知，经过参数化修正后，无论是梁格模型或精细化梁-壳单元模型与实测频率的吻合度均较为良好. 但由表4和表6可知，两种模型修正前后参数的变化量差异显著. 具体分析如下.

针对梁格模型，其修正参数特别是主梁及桥台的虚拟横梁弹性模量增大了69倍. 这表明梁格法中横梁采用基于真实截面与材料特性换算得到的抗弯刚度EI，难以准确地模拟结构的扭转刚度[13]. 针对梁-壳单元模型： ①修正后，主梁弹性模量增大了18%，可能补偿了初始模型中将箱型主梁板厚用均值代替引起的误差； ②修正后，桥墩和桥台弹性模量分别增大了34.2%和56.6%，可能补偿了采用土弹簧简化方法模拟结构-土相互作用引起的误差； ③修正后，横隔梁弹性模量降低了78%，可能因为真实结构横隔与图纸存在差异，而修正后的横隔梁弹性模量值补偿了该误差； ④修正后护栏弹性模量增大了42.9%，考虑到初始值为材料真实弹性模量的70%，修正后数值接近结构的真实弹性模量，说明该桥护栏连续性较好、刚度贡献较为突出[7, 14-15].

5 结语

1) 基于梁格法建立的动力有限元模型，对扭转模态存在显著的固有模型结构误差. 经过模型参数修正后，虽然能获得与实测频率较为接近的结果，但虚拟横梁的物理参数显著不合理，难以作为考虑环境因素影响的长期健康监测有限元基准模型.

2) 精细化梁-壳单元模型能较为合理的拟合实测模态结果，包括自振频率和振型. 其上部结构物理参数基本合理，可考虑引入温度效应对材料属性的影响，以此得到消除温度影响的高精度健康监测基准有限元模型. 相关问题有待深入研究.

3) 精细化基准模型的建立，需要综合考虑模型的精度和效率. 对于该整体桥，建立的梁-壳单元模型，通过选取合理的修正参数，采用基于参数灵敏度分析的模型修正方法，其计算结果精确、可靠.