引言

由于长期腐蚀、疲劳、自然灾害或人为事故等原因^[1]，悬索桥主缆可能出现一定程度的损伤.国内一些学者开展了带损伤桥梁的相关研究，但大多是针对悬索桥结构损伤识别与健康监测^[2-5]，或缆索的腐蚀试验与劣化机制分析等^[6]，并未揭示带损伤悬索桥梁的动力学特征.一旦缆索发生损伤，将很可能导致原本对称的结构体系两侧刚度不再一致，从而影响结构的动力特性.

一些国外学者通过建立悬索桥的连续体模型开展了此类动力学问题的相关研究.Malík等^[7]提出了一种原塔科马悬索桥的整桥模型，通过分析周期外力作用下的振动特性，揭示了该桥倒塌的重要原因是缆索的松脱.Materazzi和Ubertini^[8]使用连续介质公式研究了主缆损坏的悬索桥的垂直振动，并考虑主缆损伤的强度、位置和范围参数的影响，对固有频率和振型进行了敏感性分析.在此基础上，Ubertini^[9]建立了一根主缆上存在任意损伤的整桥模型，应用Hamilton原理得到控制方程，研究了非线性耦合垂直和扭转振动.Lepidi，Gattulli和Vestroni^[10]建立了一个受损主缆的弹性连续模型，将主缆损伤描述为轴向刚度的降低，结果表明损伤引起的截面对称性破坏，导致频率轨迹交叉转化为跃迁现象的出现.

相比于整桥模型，横截面模型更简单，可以更高效地开展各类大跨度桥梁的定性分析.Gattulli等^[11-12]提出了一种截面模型来研究悬索桥中梁和缆索之间的非线性动力作用.Plaut和Davis^[13]采用二自由度和四自由度截面模型来分析原塔科马海峡大桥的动态响应.最近，Hui等^[14]建立了模拟大跨度悬索桥截面的六自由度模型，研究了结构的非线性振动问题.这些研究都表明截面模型可以体现结构振动的基本特征.

增量谐波平衡法（IHB法）最初是由Lau和Cheung^[15-16]提出的非线性动力计算方法，它是谐波平衡法和Galerkin法的结合.最近，Hui ^[17]提出了拓展增量谐波平衡法，该方法无需对非线性刚度进行泰勒展开，从而避免了由此带来的计算误差.

为了定性分析主缆刚度损伤对悬索桥动力特征的影响，首先，以单侧主缆刚度的折减来表示损伤，基于文献^[17]中Hui所用的悬索桥模型建立了七自由度截面模型.其次，通过对系统特征方程进行模态分析，讨论了由于单侧主缆刚度折减对系统模态跃迁效应的影响.最后，利用拓展的增量谐波平衡法（EIHB法），对悬索桥模型在1∶1内共振条件下的非线性振动问题开展了数值分析.

1 悬索桥横截面模型

1.1 模型概述

本文以Hui^[17]建立的七自由度截面模型为计算模型，如图1所示.该模型由一个长方形刚体、2个质点和7根弹簧组成.长方形刚体代表主梁，具备两个方向的平动及扭转自由度.2个质点分别表示两侧缆索的质量，各自具备平面内水平和竖直方向的两个自由度.7根弹簧分为四类:S-1提供主缆的面内刚度; S-2提供吊杆的面内刚度; S-3、S-4和S-5提供主梁的竖直、水平和扭转刚度.该模型可用于研究悬索桥主梁和缆索的振动情况，并可考虑缆索和主梁之间的作用，及主缆损伤的影响.

H是主缆在竖直方向上的投影; D是主梁宽度的一半; l_i（i=1，2）是弹簧S-i（i= 1，2）的长度; m是主缆质量; M是主梁质量; J是主梁惯性矩; k_1j （j=1，2）是弹簧S-1两侧的刚度; k_i（i=2，3，4，5）是弹簧S-i（i= 2，3，4，5）的刚度; u_j（j=1，2）是m的水平自由度; v_j（j=1，2）是m的竖直自由度; x是M的水平自由度; y是M的竖直自由度; θ是M的转动自由度; F_cs是静态平衡下S-1的内力; F_hs是静态平衡下S-2的内力; F_cd是弹簧S-1的动内力; F_hd是弹簧S-2的动内力.

1.2 模型参数

针对主缆损伤导致截面刚度不对称状态下的悬索桥动力问题，以弹簧S-1的刚度折减系数表示单侧主缆损伤的程度，设右侧主缆受损，其与左侧完好主缆的刚度之比可表示为:

建模中认为$l_1\gg D\gg l_2$和$k_2\gg k_{11}\left(k_{12}\right)$，悬索桥缆索的质量m远小于主梁的质量M，m=βM，主梁惯性矩J=γMD².根据文献^[17]，除主缆刚度k₁₁和k₁₂以外，其他参数的取值为:

$\begin{array}{c}D=1,M=0.015,l_1=6,l_2=0.2,\beta =0.05,\\ \gamma =0.5,k_2=100,k_3=0.605,k_4=40,k_5=0.175\end{array}$

1.3 振动方程

上述截面模型中共有7个自由度，分别用u_j（j=1，2）和v_j（j=1，2）表示两侧主缆水平和竖直方向的自由度，x、y和θ表示主梁的三个自由度.则无阻尼自由振动方程为:

$J\ddot{\theta }+k_4\theta -\left(F_{hs}+F_{hd1}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\phi }_{h1}\right)D\mathrm{c}\mathrm{o}\mathrm{s}\theta +$

式中，F_cs和F_hs分别是弹簧S-1和S-2在重力作用下静态平衡的内力，F_cdj和F_hdj（j=1，2）分别是弹簧S-1和S-2的动内力，j表示弹簧的位置（j=1代表主梁的左侧，j=2为右侧）.

式中，$l_i^j$和$l_i^0$是弹簧S-i（i=1，2）的实际长度和初始长度.

式（2）中，φ_cj和φ_hj （j=1，2）分别为模型振动过程中弹簧S-1和S-2与竖直方向的夹角，可以用系统自由度（y，x，θ，u₁，v₁，u₂，v₂）表示，限于篇幅，在文后附录中给出.

2 模态分析

为了求系统特征方程，即方程（2）所对应的线性化振动方程，将方程（2）泰勒展开，仅保留第一阶级数，可得:

式中，M、K和q分别是质量矩阵、刚度矩阵和解向量.

$\mathit{q}=\left[\begin{array}{c}yx\theta u_1{v}_1{u}_2{v}_2\end{array}\right]$

$\mathit{M}=\left[\begin{array}{c}M\\ M\\ J\\ m\\ m\\ m\\ m\end{array}\right]$

其中，

$\begin{array}{c}{k}_{12}={\eta }_k{k}_{11},{\stackrel{-}{k}}_2=k_2-\frac{F_{hs}}{l_2},\\ {\stackrel{-}{k}}_{11}=k_{11}-\frac{F_{cs}}{l_1},{\stackrel{-}{k}}_{12}=k_{12}-\frac{F_{cs}}{l_1},\\ b_{11}=\frac{F_{cs}}{k_{11}{l}_1-F_{cs}},b_{21}=1+\frac{F_{cs}}{k_{11}{l}_1-F_{cs}},\\ b_{12}=\frac{F_{cs}}{k_{12}{l}_1-F_{cs}},b_{22}=1+\frac{F_{cs}}{k_{12}{l}_1-F_{cs}},\\ a_1=1+\frac{F_{hs}}{k_2{l}_2-F_{hs}},a_2=\frac{F_{hs}}{k_2{l}_2-F_{hs}},\\ a_3=D+\frac{F_{hs}D}{k_2{l}_2-F_{hs}}.\end{array}$

系统的特征方程如下所示:

式中，λ_k和Ф_k分别为特征值和特征向量.

通过特征值分析^[18]，可以得到系统的固有频率与振型.图2（a）和图2（b）为系统的前两阶振型，分别与桥面主梁的垂直振动和扭转振动有关.

图3给出了前两阶固有频率随主缆刚度的变化，从图中可看出，在k₁₁=0.2位置附近，两条频率曲线非常接近，使得悬索桥模型的前两阶模态之间存在着发生1∶1内共振的条件.

图3（a）η_k = 1.0时，两条固有频率曲线相交叉，k₁₁ ＜ 0.1964时，第一阶振型为扭转模态，第二阶振型为竖弯模态; k₁₁＞0.1965时，前两阶振型则相反.然而，η_k = 0.95时，频率曲线发生模态跃迁，图3（b）中振型图的变化展示了模态跃迁的过程:k₁₁＜0.16时，较为明显地，第一阶振型为扭转模态，第二阶振型为竖弯模态; k₁₁＞0.24时，前两阶振型则相反; 而在内共振条件（前两阶固有频率非常接近）下，振型表现出模态耦合的现象，耦合状态下的振型均同时具备竖弯与扭转振动的特征，且主缆刚度k₁₁在[0.16，0.24]的范围内，固有频率之间越接近，模态耦合越明显.将这两个振型根据频率的高低分别命名为mode-1和mode-2（频率较高者）.

为确定主缆损伤对固有频率的影响，进行参数η_k的敏感性分析.考虑四种不同的η_k:1.0、0.98、0.95和0.92.如图4所示，仅η_k = 1.0时，主缆未损伤，前两阶频率曲线存在交叉点; 而η_k ＜ 1.0，右侧主缆损伤，第二阶频率曲线通过交叉点，且η_k越小，损伤程度越大，前两阶频率曲线的间距增大，模态耦合现象将在主缆刚度k₁₁更大的区间范围内出现.

3 EIHB算法

有阻尼系统在简谐激励下的振动方程为:

式中，F为外部谐波激励的幅值向量，ω表示激励频率; F_R（q）为系统的恢复力向量; M和C分别为系统的质量矩阵和阻尼矩阵，采用结构计算中较为常用的瑞利阻尼，$\mathit{C}={\alpha }_0\mathit{M}+{\alpha }_1\mathit{K}，{\alpha }_0=0.005，{\alpha }_1=0.001$; 当η_k = 1.0时，竖弯和扭转模态的阻尼比分别为0.31%和0.32%.

激励频率为ω₀时，设方程（6）的位移响应为q₀.引入位移和激励频率的增量形式后，近似解为:

将式（7）代入方程（6）时，可以得到:

式中，τ = ωt，撇号（′）表示对变量τ的微分; R是误差向量; $\frac{{\mathit{F}}_{\mathit{R}}\left({\mathit{q}}_0+\mathrm{\Delta }\mathit{q}\right)-{\mathit{F}}_{\mathit{R}}\left({\mathit{q}}_0\right)}{\mathrm{\Delta }\mathit{q}}$为刚度，不再局限于多项式形式的非线性刚度.

位移响应设为傅里叶级数形式:

式中，$\mathit{S}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(C_{\mathrm{s}}，C_{\mathrm{s}}，C_{\mathrm{s}}，C_{\mathrm{s}}，C_{\mathrm{s}}，C_{\mathrm{s}}，C\mathrm{s}\right)，\mathit{A}={\left[A_1^{\mathrm{T}}，A_2^{\mathrm{T}}，\cdots ，A_7^{\mathrm{T}}\right]}^{\mathrm{T}}$，其中A_i是对应于第i自由度的振幅向量，C_s为:

经过对结果的试算和验证，式中，令m=4，激励频率等于基频的分量为一次谐波项，激励频率为基频的2倍、3倍和4倍的分量分别为二次、三次和四次谐波项，常数项为静态项.

将式（10）代入方程（8）中，并在τ为0到2π的一个周期内应用Galerkin平均过程，可得代数方程:

其中，

$\begin{array}{c}\tilde{M}={\int }_0^{2\pi }{\mathit{S}}^{\mathrm{T}}M{\mathit{S}}^{\mathrm{\text{'}}\mathrm{\text{'}}}d\tau ,\\ \tilde{C}={\int }_0^{2\pi }{\mathit{S}}^{\mathrm{T}}C{\mathit{S}}^{\mathrm{\text{'}}}d\tau \\ \tilde{K}={\int }_0^{2\pi }{\mathit{S}}^T\frac{\partial {\mathit{F}}_{\mathrm{R}}\left(\mathit{q}\right)}{\partial \mathit{q}}Sd\tau ,\\ \stackrel{-}{R}={\int }_0^{2\pi }{\mathit{S}}^{\mathrm{T}}Fcos\left(\tau \right)d\tau -\left({\omega }_0^2\tilde{M}+{\omega }_0\tilde{C}\right)A-\\ {\int }_0^{2\pi }{\mathit{S}}^{\mathrm{T}}{\mathit{F}}_{\mathrm{R}}\left({\mathit{q}}_0\right)d\tau ,\\ {\stackrel{-}{K}}_{\mathrm{m}\mathrm{c}}={\omega }_0^2\tilde{M}+{\omega }_0\tilde{C}+\tilde{K}\text{,}\\ {\stackrel{-}{R}}_{\mathrm{m}\mathrm{c}}=\left(2{\omega }_0\tilde{M}+\tilde{C}\right)A\\ \end{array}$

采用弧长法进行数值求解^[19]，在计算时应先确定初始估计值A和ω₀，指定一个增量ΔA，代入方程（12）求解得到Δω，然后判断R的值是否满足给定精度要求，不断循环直至得到新的解.重复迭代过程，得到系统在不同激励频率下的响应.

4 1∶1内共振

本文主要研究1∶1内共振条件下悬索桥的非线性振动，考虑两种不同η_k的工况:1）悬索桥单侧主缆刚度有损伤，η_k = 0.95，k₁₁=0.2，k₁₂=0.19，此时，系统的前两阶固有频率之比ω₁/ω₂=5.3075/5.3789=0.9867; 2）悬索桥主缆刚度未损伤，η_k= 1.0，k₁₁= k₁₂= 0.2，固有频率之比ω₁/ω₂=5.3730/5.3873=0.9973.

4.1 竖向简谐激励

在主梁的y自由度上施加简谐力，主梁竖向位移的幅-频响应曲线如图5所示.结果显示，一次谐波项占绝对主导地位，静态项，二次、三次和四次谐波项的幅值均比一次谐波项的幅值小得多，可以忽略不计.

图5中蓝色曲线为η_k=0.95工况下的系统响应，各频率分量均含有两个峰值，这是因为共振条件下有损伤的悬索桥发生模态跃迁，两个原本相互独立的振型相互耦合.图中左侧峰值为mode-1的响应幅值，右侧峰值为mode-2的幅值.η_k=1.0工况对应的黑色频响曲线仅一次谐波项有较大响应幅值，对其他频率分量的激励不敏感.

图5中线性化系统表示不考虑系统非线性刚度的情况（即方程（8）中刚度项替换为方程（4）中的线性刚度）.相较于线性化系统的结果，η_k=0.95时，mode-1的部分能量转移到mode-2的竖向振动中，系统发生了1∶1内共振; 而η_k=1.0时，两条曲线吻合一致，结果表明，未损伤的悬索桥即使固有频率非常接近也不能激发内共振现象.

图6给出了竖向激励作用下主梁扭转振动的幅-频响应曲线（一次谐波项）.η_k=1.0工况下没有扭转响应，而η_k=0.95工况下，有损伤的悬索桥前两阶模态mode-1和mode-2均有明显的扭转响应.

4.2 简谐力矩激励

在主梁的θ自由度上施加简谐力矩，系统响应计算结果如图7和图8所示.与竖向激励下的响应类似，扭转激励的一次谐波项为主要成分，限于篇幅，下文仅给出一次谐波项的结果.

简谐力矩激励引起的扭转振动如图7所示，频响曲线中可观察到“软化”效应，即峰值处的响应频率倾向于低频率方向.这是因为主梁扭转刚度的非线性随振幅增加逐渐增大，模型表现为“软弹簧”，发生共振时的频率呈降低的趋势.

主梁竖直振动的幅-频响应曲线如图8所示.η_k=1.0工况下，系统没有竖向响应，未能激发1∶1内共振现象; 而η_k=0.95工况下，悬索桥前两阶模态均有明显的竖向响应.相较于线性系统，非线性系统的内共振导致mode-1的扭转响应幅值明显增大，系统能量发生转移.

图5和图7的幅-频响应结果表明，单侧主缆刚度有损伤悬索桥的响应幅值较未损伤工况有所减小，但是未损伤的工况仅在外激励频率等于自振频率时有较大共振响应，而有损伤的悬索桥具有两个相邻的频响峰值，会对激励频率更为敏感，导致共振频带更宽.

4.3 响应对比

为验证EIHB法计算结果的准确性，采用龙格-库塔（Runge-Kutta，R-K）法比较系统的非线性动力响应.R-K法是一种常用的时域方法.限于篇幅，仅给出激励幅值a=0.001的两组结果进行对比，如图9所示.图9（a）为竖向激励频率ω=5.380 rad/s时位移y的相位图，图9（b）为扭转激励频率ω=5.308 rad/s时转角θ的相位图.EIHB法和R-K法的响应曲线均吻合一致，证明了所得数值结果准确可靠.

5 结论

本文基于Hui的悬索桥七自由度截面模型，考虑到主缆刚度损伤在单侧主缆上引入折减系数η_k，通过模态分析来研究模态跃迁过程，并采用拓展的增量谐波平衡法（EIHB法）对模型进行了动力分析，得到以下主要结论:

1）随主缆刚度变化，系统的前两阶固有频率曲线发生交叉或跃迁.主缆未受损，刚度折减系数η_k为1.0时，竖弯与扭转模态的频率曲线发生交叉，而单侧主缆受损状态下，频率曲线发生模态跃迁.

2）单侧主缆刚度有损伤的悬索桥较未损伤工况响应幅值明显减小，但未损伤的悬索桥仅在外激励频率等于自振频率时有较大共振响应，而有损伤的悬索桥有两个共振响应峰值，对激励频率更为敏感.

3）单侧主缆刚度有损伤的悬索桥发生了能量转移，而未损伤的悬索桥没有激发1∶1内共振.

附录

φ_cj和φ_hj （j=1，2）的正弦和余弦计算公式为:

$\begin{array}{c}{G}_1=(M/2+m)g,G_2=M/2g,\\ B_1=l_1^0{F}_x/\sqrt{F_x^2+G_1^2}+F_x/k_1,\end{array}$

$\begin{array}{c}{B}_2=l_2^0{F}_x/\sqrt{F_x^2+G_2^2}+F_x/k_2,\\ C_{11}=l_1+G_1/k_{11},C_{12}=l_2+G_1/k_{12},\\ C_2=l_2{G}_2+G_2/k_2,\end{array}$

$\begin{array}{c}{L}_{h1}=\sqrt{{\left(C_2+v_1-y-D\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\right)}^2+{\left(B_2-u_1+x+D(1-\mathrm{c}\mathrm{o}\mathrm{s}(\theta \left)\right)\right)}^2},\\ L_{h2}=\sqrt{{\left(C_2+v_2-y+D\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\right)}^2+{\left(B_2-u_2-x+D(1-\mathrm{c}\mathrm{o}\mathrm{s}(\theta \left)\right)\right)}^2},\end{array}$

$\begin{array}{c}cos\left({\phi }_{h1}\right)=\left(C_2+v_1-y-D\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\right)/L_{h1},\\ sin\left({\phi }_{h1}\right)=\left(B_2-u_1+x+D(1-\mathrm{c}\mathrm{o}\mathrm{s}(\theta \left)\right)\right)/L_{h1},\\ cos\left({\phi }_{h2}\right)=\left(C_2+v_2-y+D\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\right)/L_{h2},\\ sin\left({\phi }_{h2}\right)=\left(B_2-u_2-x+D(1-\mathrm{c}\mathrm{o}\mathrm{s}(\theta \left)\right)\right)/L_{h2},\\ L_{c1}=\sqrt{{\left(C_{11}-v_1\right)}^2+{\left(B_1+u_1\right)}^2},\end{array}$

$\begin{array}{c}{L}_{c2}=\sqrt{{\left(C_{12}-v_2\right)}^2+{\left(B_1+u_2\right)}^2},\\ cos\left({\phi }_{c1}\right)=\left(C_{11}-v_1\right)/L_{c1},\\ sin\left({\phi }_{c1}\right)=\left(B_1+u_1\right)/L_{c1},\\ cos\left({\phi }_{c2}\right)=\left(C_{12}-v_2\right)/L_{c2},\\ sin\left({\phi }_{c2}\right)=\left(B_1+u_2\right)/L_{c2}.\end{array}$

参考文献