斜拉桥面内竖向固有振动模型及特性影响的有限差分分析

陈柯帆 ，李源 ，2†，贺拴海 ，2，王康 ，殷怡萍 ，宋一凡 ，2

（1.长安大学 公路学院，陕西 西安 710064；2.旧桥检测与加固技术交通行业重点实验室（长安大学），陕西 西安 710064）

斜拉桥美观、经济、跨越能力强，近年来备受桥梁工程师青睐［1］.与此同时，斜拉桥整体结构复杂，柔度大、阻尼低、刚度不足导致其非线性行为极为突出［2］.尤其当拉索局部模态与斜拉桥整体模态频率比处于“1∶2”或“1∶1”等固定比例区间时［3-6］，容易在环境激励下，引发拉索剧烈振动，给桥梁安全运营带来了极大隐患.因此，除了目前仅有的有限元方法外［7-9］，如何建立准确的斜拉桥整体动力模型来便捷而又准确地计算斜拉桥竖向整体模态参数，对推动斜拉桥应用与发展至关重要［10］.

近年来，国内外学者围绕斜拉桥整体动力学建模方法，尤其对于斜拉桥主梁在多点弹性支撑作用下的力学行为模拟与分析做了大量研究工作.吴庆雄等［11］进行了单索-梁结构和二索-梁结构模型固有振动试验，建立了多索-梁结构动力学模型，讨论了斜拉索对索梁结构面内固有振动特性的影响；Cao等［12］和李专干等［13］通过建立主梁的分段函数，将主梁和拉索等效为若干独立梁段，基于拉索锚固处的边界条件得到了刚塔柔梁斜拉桥整体动力学建模的运动方程，讨论了结构对称性对动力特性的影响；赵文忠等［14］通过分段函数求解了三索结构的动力方程，研究了拉索一阶和二阶频率比条件对共振的影响；Cong 等［15］、Kang 等［16-18］、苏潇阳等［19］建立了多梁弹簧动力学模型，运用传递矩阵法给出了不同体系下斜拉桥整体动力学模型动力学微分控制方程，提出了不同体系斜拉桥竖向刚度评估方法.在此基础上，该课题组还对索拱结构［20-21］、悬索结构［22］的面内外固有振动模态参数进行了系统分析与研究.

受限于现有动力学建模方法与整体模态参数计算过程的烦冗，目前关于斜拉桥整体竖向模态的理论计算方法大多高度简化甚至忽略索力水平投影、振动时拉索间的耦合影响作用、主梁弯曲刚度等因素对结构振动特性的影响，或是迭代解析方法烦琐、复杂，不利于在实际斜拉桥工程中推广与应用.

针对此问题，本文通过离散斜拉桥多点弹性支承梁的集中质量参数体系，建立了一种新的斜拉桥面内竖向整体动力学模型.该模型考虑了斜拉索对主梁的竖向弹性支承作用与对水平梁截面的轴力影响，引入微梁段两侧的剪力以模拟主梁弯曲刚度、拉索间振动耦合等影响作用，通过微梁段间的弯矩平衡和有限差分法，修正了不同结构体系斜拉桥面内竖向整体动力学模型的运动方程，结合特征值法给出了斜拉桥面内竖向模态频率及振型计算方法.通过对比参考文献中动力学模型算法案例分析结果与某斜拉桥的实测值，进一步验证了本文关于斜拉桥面内竖向运动参数体系建模方法的适用性和正确性.本文计算方法无须建立大量细化的有限元模型，运用MATLAB、Excel 等软件按编码流程即可准确、快速地估算斜拉桥面内竖向模态频率及振型，简化了计算过程，便于工程应用.

1 多点弹性支撑梁的离散模型

约定下标“B”和“C”分别表示梁和索；下标“i”和“j”分别表示索和梁序号（i∈[1，I]，j∈[1，J]）.为便于区别有索区梁段与无索区梁段，对于Ci#拉索锚固处的有索区梁段表示为Bji#梁段.建立如图1（a）所示的多索-主梁模型.由于桥面质量远大于拉索质量，本文忽略了拉索振动对主梁振动的影响，将Ci#拉索竖向视为Bji#梁段的弹性支承，水平向的索力投影等效为Bji#梁段的轴向荷载HBCi，定义He为支座水平力，如图1（b）所示.为了精确模拟和求解具有分布质量、荷载和多点弹性支撑作用下斜拉桥主梁的动力行为，不考虑主梁的纵向运动，本文将主梁进一步简化为J个间距d相同、彼此铰接、带有I个竖向弹力支承和I个不同轴力的集中质量参数体系，主梁的质量和荷载均被视为作用在这些集中质量点之上，独立梁段间彼此通过理想铰连接，如图1（c）所示.

图1 斜拉桥主梁离散模型的简化流程Fig.1 The reduction process of main beam of cable-stayed bridges

图1（b）中，Ci#拉索对Bji#梁段竖向弹性支承系数kBCi为［7，10，23-24］：

式中：ECi、ACi分别为Ci#拉索的弹性模量及截面积；θCi表示Ci#拉索轴线与主梁大里程方向夹角；lCi表示 Ci#拉索上下端锚固点轴向距离.由于拉索振动中的索力增量远小于拉索初始索力，对于拉索索力的水平投影本文仅考虑初始索力［25-26］.因此，主梁上水平轴力HBCi和支座水平力HBe（下标e表示边界）满足：

式中：SCi表示Ci#拉索初始索力.图1（c）中，独立梁段彼此铰接，通过引入离散梁段的左右侧剪力以模拟具有分布质量的主梁竖向弯曲刚度、拉索间振动耦合作用在振动过程中的相互影响.梁段间的受力如图2（a）所示，梁段处的受力如图2（b）和（c）所示.

图2 微梁段受力示意图Fig.2 The force schematic of the micro-beam segment

图2 中，FB（j-1，j）、FB（j，j+1）表示梁段左右两侧的剪力；NB（j-1，j）、NB（j，j+1）表示梁段左右两侧的轴力-表示相邻梁段间的节段左右两侧受到的弯矩作用表示相邻梁段运动夹角；aBj表示梁段运动加速度.考虑系统初始为平衡状态，基于D’Alembert 原理可以得到Bj#梁段在竖向的动力平衡方程：

式中：VBj（t）表示Bj#梁段与时间相关的竖向振动位移变化因子，后文中简写为VBj；kT表示斜拉桥主塔对主梁面内竖向自由运动的刚度弹簧系数，需依据图纸和规范，以及不同结构体系下塔-梁处的边界条件进行取值.δ(j-ji)为狄拉克（Dirac）函数，由式（5）～式（6）定义：

值得注意的是，靠近边界的梁段需根据结构体系边界条件进行求解运算，如附表1所示.

附表1 不同主梁边界条件下的振动方程系数表达式Additional Table 1 Expressions of vibration equation coefficients under different boundary conditions of main beams

2 数学表达与验证

2.1 基于有限差分法的方程优化

假设质量体系分布较密，梁段间相对位移较小，则其振动的几何关系满足以下关系式：

图2（a）中，显然梁段间存在剪力与弯矩平衡：

由上式可得梁段间剪力、弯矩与振动位移间关系：

对于无索区梁段（即j≠ji），左右侧截面受到同一方向常轴力影响，可以得到：

对于有索区梁段（即j=ji），左右侧截面轴力突变量为拉力的水平投影：

整合式（7）～式（15）并代入式（4）后，可以得到Bj#梁段的振动方程：

整合B1#～BJ#梁段方程后可以得到斜拉桥面内竖向固有振动方程：

式中：Γ为等效剪力效应系数矩阵，表征了剪力效应对主梁运动的影响，与斜拉桥边界条件有关；Ξ为等效轴力效应系数矩阵，表征了轴力效应对主梁运动的影响，与拉索数量和锚固位置有关，矩阵具体形式如附录1 所示皆为相同形式的J维列向量.为避免赘述，在此仅展示形式：

式（20）中：ψBj单位与频率一致，是Bj#梁段的局部模态频率方程，表征了Bj#梁段参与整体模态的模态频率，由κDj、κCi、κSj、κT构成.κDj表示索力水平投影，即主梁轴力对Bj#梁段局部模态频率的影响；κCi表示拉索竖向弹性支承作用的影响；κSj表示主梁轴力的影响，系数λj与斜拉桥边界条件有关，具体形式如附表1所示；κT表征了不同结构体系主塔的影响.

2.2 基于特征值法的结构模态分析

参数体系的自由振动方程组——式（17）实际上是一个J维的齐次方程组，方程有解的前提是系数矩阵行列式为0.因此，不计结构阻尼，根据式（17）构造结构特征矩阵形式如下：

式中：MBj表示质量对角矩阵；KBj表示主梁的等效刚度矩阵.根据简谐振动理论，构造频率表达式：

式中：eig 表示求解矩阵特征值与特征向量.采用以上计算公式和流程求解斜拉桥的面内竖向固有振动模态参数，在确定质量矩阵并根据相应结构边界条件选择和编辑相应轴力与剪力系数矩阵后，仅需借助MATLAB、Excel 等工具就能简单地计算和求解固有振动频率和振型，无须进行大量细化的有限元建模分析.本文依托MATLAB编制了运行算法程序，其流程如图3所示.

图3 斜拉桥面内竖向固有振动模态参数计算流程Fig.3 The computation solution process of a cable-stayed bridge’s in-plane vertical natural vibration modal properties

2.3 参考文献算例验证

以2020年Cong 团队的双索-梁结构理论解析研究结果为对比对象［14，28］，该研究通过简化斜拉桥为多索-梁动力学模型，考虑结构构件间的几何非线性边界条件，基于传递矩阵法研究了多索斜拉结构的固有振动特性与面内外振动响应.其研究模型如图4所示.

图4 文献［14，28］双索-梁结构示意图Fig.4 Schematic of a double-cable-beam structure in the references ［14，28］

本文代入了文献［14，28］参数，选取主梁划分节段参数d=1 m，即J=300，I=2，运用MATLAB 软件根据图3 计算流程编写计算程序，讨论和对比多索斜拉结构的面内固有振动频率及以主梁为主要振型的前五阶振动模态，如表1所示.

表1 两索-梁结构面内竖向固有振动频率Tab.1 In-pane vertical natural vibration frequencies of the two-cable-beam structure Hz

根据表1，本文解析法得到的频率数值平均绝对误差仅为0.2‰，精确度超过了原文解析方法绝对误差1.4%.分别按照本文有限元数值模拟与解析法求解该结构以主梁为主的前五阶固有振型，如图5所示.

图5 两索-梁结构的竖向振型Fig.5 Vertical modal shapes of the two cable-beam structure

图5 显示两种方法得到的结构振型一致性良好，进一步验证了本文计算方法的有效性和准确性.

3 斜拉桥固有振动特性及影响性分析

3.1 斜拉桥模态分析

以我国西北地区某混凝土斜拉桥为对象开展固有振动特性及其影响参数分析.该桥全长166.8 m（39 m+88.8 m+39 m），采用三跨双台、双塔、双索面对称布置，墩塔处固结，为半漂浮支承体系，其立面图如图6所示.钢筋混凝土主梁由节段预制双箱梁和预制行车道板组合形成，箱梁高1.2 m，桥面净宽8.5 m.在两箱梁间锚固板处设横系梁一道，纵向长约0.22 m，为方便引用，汇总主梁各截面参数设置如表2 所示；全桥现有48 根斜拉索，从小里程边跨至大里程边跨方向以C1#～C24#对单索面拉索依次编号，参数如表3所示（仅示出一侧，另一侧参数与之相近，斜拉索弹性模量经恩斯特公式修正后取200 GPa）.

表2 主梁参数设置表Tab.2 Parameters of the beam

表3 斜拉索参数设置表Tab.3 Parameters of stay cables

图6 桥梁立面图（单位：cm）Fig.6 Bridge elevation drawing（unit：mm）

根据桥梁结构形式，采用商业有限元软件对该桥进行动力特性分析，得到该结构前三阶自振频率、振型特征，如图7所示.

图7 有限元法得到的斜拉桥前三阶固有振动模态参数Fig.7 The first-three order natural modes of the cable-stayed bridge by the finite element model

为进一步对比和验证本文解析公式的正确性，采用实桥实测、有限元分析、本文解析三种方法对该桥进行自振特性分析.设置参考算例参数详情如下所示：

1#算例（Referred Case 1#，RC1），实体结构有限元法：根据实际桥梁结构形式，采用商业有限元软件对该桥进行动力特性分析，识别该结构自振频率、振型特征等.

2#算例（Referred Case 2#，RC2），现场实测法：在桥梁边跨0.4L（L表示跨径）截面及中跨跨中截面布设加速度传感器，应用脉动激励法进行桥梁结构的振动试验，识别大桥前3 阶整体模态的动力特性参数，采用DHSAS 频谱分析及模态分析软件对其进行快速傅里叶变换得到相应的功率谱图，再对其作进一步的频谱分析可得到桥梁结构的自振频率、阻尼比.现场动载试验布置如图8所示.

图8 现场动载试验Fig.8 On-site dynamic load test

3#算例（Referred Case 3#，RC3），简化模型解析法：根据本文动力简化模型及振动方程，塔梁连接处按照图纸取kT=5.3 × 109N/m，大小里程结合墩采用简支边界条件，位于主梁同一截面的双索考虑为动力弹簧的并联关系，基于式（17）采取微梁段长度d1=0.1 m，运用MATLAB 软件根据图3 计算流程编写计算程序，对该结构的面内竖向模态频率和振型进行了计算和分析.

汇总以上工况下得到的该桥以主梁面内竖向为主要振型的前三阶自振频率，如表4所示.

表4 斜拉桥面内竖向固有振动频率Tab.4 In-plane vertical natural vibration frequencies of the cable-stayed bridge Hz

表4 中，若以实桥测得的模态参数为标准，有限元法得到的结果的绝对误差平均值为3.3%，而本文解析法得到的结果的绝对误差平均值为2.7%.此外，相较于其他两种结果，实测频率值整体偏小，这是由于该桥建设时间较长，结构刚度在通行运营中有所下降.汇总RC1 和RC3 工况下得到的该斜拉桥前三阶竖向固有振型，如图9所示.

图9 RC1、RC3工况下的斜拉桥前三阶面内竖向振型Fig.9 The first three order modes of cable-stayed bridge under the conditions of RC1 and RC3

图9 中两种工况下结构面内竖向前三阶振型一致性良好，上述情况进一步说明了本文方法计算斜拉桥面内竖向固有振动模态参数的精确性和适用性.

3.2 梁轴力对斜拉桥竖向模态的影响

为研究索力水平投影为梁提供的轴向力对斜拉桥面内竖向振动模态参数的影响，引入μa表示轴向力对Bj#梁段动平衡方程的放大系数，由式（28）定义：

按照本文解析法研究成果，基于RC3参数设置，图10展示了轴力对斜拉桥固有振型的影响，图11展示了轴力对斜拉桥固有振动频率的影响.

图10 轴力对斜拉桥竖向固有振型的影响Fig.10 Influence of the axial force on the vertical natural modes of the cable-stayed bridge

图11 梁轴力对斜拉桥面内竖向频率的影响Fig.11 Influence trends of the beam’s axial force on the in-plane vertical frequencies of the cable-stayed bridge

图12 μa与结构前五阶频率关系Fig.12 Relation between the first-five order in-plane vertical frequencies with μa

图12 显示，μa<28 时，V1～V5 随轴力增大而呈现线性缓慢下降.当μa达到30 左右时，V1 发生频率跃迁现象，此时尽管V2频率数值低于V1，但V1的振型是结构的一阶振型.当μa达到34.4时，V1～V3值产生共轭对称解，其频域信号值幅值相同而相位不同，而V2 在μa接近38.7 时急速下降为0，此时斜拉桥一阶固有振动频率为0，结构整体失稳.表明从轴力角度考虑，该桥目前暂无整体失稳风险，只有当轴力达到现有轴力38 倍以后会出现整体失稳现象.此外，随着结构轴力的增大，结构基频将发生跃迁，结构体系内易产生非线性内共振，需进行更深入的研究.

3.3 断索对斜拉桥面内竖向模态的影响

从式（20）与式（22）中可发现，斜拉桥拉索为主梁提供了轴力与弹性支撑作用（kBCi），为结构提供了有效刚度.因此，结构的固有振动模态频率与每一根拉索息息相关.以此实体结构为背景，模拟单拉索出现断裂的极端情况对结构固有振动频率及振型的影响，图13显示了该工况下低阶频率变化规律.

图13 不同位置的拉索断裂对结构V1～V5频率的影响Fig.13 Influence on the V1～V5-order modal frequencies when cables of different positions broke down

拉索的索力水平投影对主梁产生的轴压力降低了整体刚度，而其竖向弹性支承作用则增大了整体刚度，两者作用下使得断索对结构基频影响较小，如图13 所示.此外，左右侧索面断索后影响效应变化规律基本一致，其中中跨拉索断裂对V1、V2和V5影响较大，而边跨拉索对V3 和V4 影响较大，对比图10中结构前五阶面内竖向振型，表明断索后的影响效应与该拉索对应的主梁质点的振型参与系数相关.拉索发生损伤甚至断裂会改变结构整体频率，因此，当此情况发生时需进一步考虑结构因局部-整体模态耦合发生非线性共振的问题.

4 结论

1）本文方法考虑了主梁截面的变刚度、变轴力作用，建立的斜拉桥整体模态动力学模型更贴近工程实际结构.通过代入已有文献中理论模型算例参数并对比其理论计算结果，本文方法计算结果平均误差率为0.2‰；代入某斜拉桥参数并对比其面内竖向固有振动频率实测值，本文方法计算结果误差率为2.7%，进一步说明本文解析方法具有较高的精确度和较好的适用性.

2）斜拉桥低阶面内竖向固有振动频率随着轴力增加而降低，振型基本无变化，当轴力增加到一定值后，结构的低阶面内竖向整体模态频率将发生频率跃迁、分叉现象.

3）斜拉桥发生断索对整体结构固有振动频率影响较小，断索对各阶频率值影响效应与拉索锚固处主梁的对应阶次振型参与系数相关.

4）本文方法为快速计算斜拉桥整体模态频率以避免斜拉桥整体-局部模态耦合而发生非线性内共振问题提供了有效帮助.下一步将考虑拉索、主塔与主梁间几何非线性边界条件，通过建立更加细化的斜拉桥整体动力学模型开展斜拉桥的相关非线性共振研究.

附录1

为简化表达，定义Pp和Pp，q表达式为：

Ξ为推导得到的子项系数矩阵，如式（a-3）所示，其与拉索数量、锚固位置有关.不同主梁边界条件下的振动方程系数表达式见附表1.

附录2

系数矩阵的具体形式与边界条件相关，本文简支-简支表达式如式（b-1）所示.

系数矩阵的具体形式与拉索数量、锚固位置、边界条件相关，本文简支-简支表达式如式（b-2）所示.