0 引言

颤振是一种由于气动不稳定引起的自激发散振动，在大跨度桥梁中若发生颤振，可能导致结构的整体破坏，因此大跨度桥梁结构颤振问题的研究显得尤为重要。为了对桥梁结构进行颤振响应分析，首先要确定桥梁断面的颤振导数。准确识别桥梁的颤振导数，成为大跨度桥梁结构颤抖振分析^[1]中的重要环节。

对于颤振机理^[2]的研究，杨詠昕^[3]在Matsumoto^[4]的研究基础上，导出了二维三自由度耦合颤振分析方法，用该方法对颤振机理问题进行了研究，但该方法未实现系统各自由度振动方程之间的解耦。丁泉顺，朱乐东^[5]通过对系统竖向和扭转耦合振动方程的解耦，将系统振动参数的频率和阻尼比隐式地表达为桥梁断面颤振导数的函数形式，导出了二自由度耦合颤振分析方法。Scanlan^[6]等最早建立了桥梁颤振的多模态分析方法，并从能量观点对桥梁的颤振稳定性进行了很有价值的研究。刘高^[7]从结构与气流系统内部能量平衡的观点对系统的颤振进行研究，发展了一种全桥多模态颤振能量分析方法。刘祖军等^[8-9]应用激励—反馈的分步分析方法对耦合颤振方程进行解耦，并将耦合颤振微分方程改造为能量方程形式，从能量转换的角度对平板弯扭耦合颤振进行了分析。朱乐东、高广中^[10-13]基于能量等效原理，建立了非线性自激力模型的等效线性化形式，通过对比非线性自激振动响应的理论预测值和试验值，验证了非线性自激力模型的适用性。本文基于能量等效原理，依据弯扭耦合颤振自激力各项对系统振动的作用效应不同将其分为纯阻尼效应项H₁^*、A₂^*，纯刚度效应项A₃^*、H₄^*和既有刚度效应又有阻尼效应的双重效应项A₁^*、H₂^*、H₃^*、A₄^*这样分类不仅明确了颤振自激力各项的物理意义，而且为从功能角度对经典耦合颤振机理进行研究提供了基础。

对于颤振导数的识别，张斐针对桥梁断面颤振导数识别的具体特点，在加权整体最小二乘法^[14]和修正总体最小二乘迭代法^[15]的基础上，引入加权矩阵，推导出用于颤振导数识别的加权最小二乘迭代(WLS)法^[16]。此外，最小二乘复指数法^[17](LSCE)，基于特征系统实现算法^[18](ERA)的识别方法，等一系列以二维节段模型自由振动试验为基础的颤振导数识别方法得到提出和发展。上述颤振导数的识别方法大多通过颤振导数与系统模态参数之间的关系，将颤振导数的识别问题转化为系统模态参数的识别问题。本文从功能转换角度，提出了一种基于瞬时做功的颤振导数识别方法。

桥梁和风是两个相互作用的系统，风对桥梁系统的作用效应可以分为两部分：(1)风的阻尼作用耗散或增加桥梁系统的机械能，即阻尼效应。(2)风与桥梁系统交换机械能，并以弹性势能的形式储存，即刚度效应。本文基于能量等效原理，首先将颤振自激力各项进行了分类。在有限元分析软件ANSYS中建立理想平板模型，其颤振自激力通过MATRIX27单元实现。以理想平板的经典耦合颤振为例，依据功能转换原理，通过对具有阻尼效应自激力的做功时程和桥梁系统机械能时程的对比分析，验证了上述颤振自激力分类的合理性。然后将颤振微分方程通过积分运算转化为功能方程的形式，提出了一种基于瞬时做功的颤振导数识别方法，并以理想平板的颤振导数识别为例对该方法的可靠性进行了验证。

1 颤振自激力分类及其合理性验证 1.1 颤振自激力分类

按照Scanlan颤振自激力表达式，桥梁断面单位长度上受到的自激升力L_se(t)和自激扭矩M_se(t)可以表示为竖向位移、扭转位移及其一阶导数的函数L_se($h, \dot{h}, \alpha, \dot{\alpha}$)和M_se($h, \dot{h}, \alpha, \dot{\alpha}$)。经典弯扭耦合颤振的运动方程表达式为:

$ \begin{array}{c} I_{m} \ddot{\alpha}+2 I_{m} \xi_{\alpha} \omega_{\alpha} \dot{\alpha}+I_{m} \omega_{a}^{2} \alpha=\rho B^{3} \omega A_{1}^{*} \dot{h}+ \\ \rho B^{3} \omega^{2} A_{4}^{*} h+\rho B^{4} \omega^{2} A_{3}^{*} \alpha+\rho B^{4} \omega A_{2}^{*} \dot{\alpha} \end{array} $

(1)

$ \begin{array}{c} m_{\mathrm{h}} \ddot{h}+2 m_{\mathrm{h}} \xi_{h} \omega_{h} \dot{h}+m_{\mathrm{h}} \omega_{h}^{2} h=\rho B^{2} \omega H_{1}^{*} \dot{h}+ \\ \rho B^{2} \omega^{2} H_{4}^{*} h+\rho B^{3} \omega^{2} H_{3}^{*} \alpha+\rho B^{3} \omega H_{2}^{*} \dot{\alpha} \end{array} $

(2)

式中：ρ为空气密度，U为风速，B为模型宽度，折算频率K=ωB/U，H_i^*和A_i^*(i=1, 2, 3, 4)是K的无量纲系数，称为颤振导数。

依据能量等效原理，当流线型桥梁断面采用Scanlan颤振自激力模型计算颤振响应时，为保证振动响应计算结果的等效性，Scanlan颤振自激力和实际颤振自激力(L_se(t)、M_se(t))的气动阻尼系数和气动刚度系数在每一个瞬时均应相等，即满足:

$ \int_{t_{1}}^{t 2} L_{\mathrm{se}}(\tau) \dot{h}(\tau) \mathrm{d} \tau=\int_{t_{1}}^{t 2} L_{\mathrm{se}}(h, \dot{h}, \alpha, \dot{\alpha}) \dot{h}(\tau) \mathrm{d} \tau $

(3)

$ \int_{t_{1}}^{t_{2}} L_{\mathrm{se}}(\tau) h(\tau) \mathrm{d} \tau=\int_{t_{1}}^{t 2} L_{\mathrm{se}}(h, \dot{h}, \alpha, \dot{\alpha}) h(\tau) \mathrm{d} \tau $

(4)

$ \int_{t_{1}}^{t_{2}} M_{\mathrm{se}}(\tau) \dot{\alpha}(\tau) \mathrm{d} \tau=\int_{t_{1}}^{t_{2}} M_{\mathrm{se}}(h, \dot{h}, \alpha, \dot{\alpha}) \dot{\alpha}(\tau) \mathrm{d} \tau $

(5)

$ \int_{t_{1}}^{t_{2}} M_{\mathrm{se}}(\tau) \alpha(\tau) \mathrm{d} \tau=\int_{t_{1}}^{t_{2}} M_{\mathrm{se}}(h, \dot{h}, \alpha, \dot{\alpha}) \alpha(\tau) \mathrm{d} \tau $

(6)

式(3)、式(5)为瞬时气动阻尼等效性的约束条件，其积分项代表自激力对节段模型振动系统做功的瞬时积累值，若保证自激力做功的当前积累值等效，那么节段模型振动系统的瞬时机械能也是等效的，由此，可以保证计算出的瞬时振幅是相等的；式(4)、式(6)为瞬时气动刚度等效性的约束条件，其积分项代表自激力对节段模型振动系统所做的无功瞬时积累值，若保证自激力所做无功的当前积累值等效，则可以保证计算出的瞬时相位是一致的。在自激振动过程中，由于气动自激力耦合项的存在，扭转模态和竖弯模态的振动响应都包含两个不同的频率成分，竖弯振动和扭转振动的位移与速度函数可以近似表达为:

$ {h(t) = {a_{h1}}\cos {\varphi _{h1}} + {a_{h2}}\cos {\varphi _{h2}}} $

(7)

$ {\dot h(t) = {a_{h1}}{\omega _{h1}}\sin {\varphi _{h1}} + {a_{h2}}{\omega _{h2}}\sin {\varphi _{h2}}} $

(8)

$ {\alpha (t) = {a_{a1}}\cos {\varphi _{a1}} + {a_{a2}}\cos {\varphi _{a2}}} $

(9)

$ {\dot \alpha (t) = {a_{a1}}{\omega _{a1}}\sin {\varphi _{a1}} + {a_{a2}}{\omega _{a2}}\sin {\varphi _{a2}}} $

(10)

经典耦合颤振属于竖弯和扭转模态完全耦合的颤振，当颤振发生时竖弯和扭转两个自由度的振动频率基本相等，并具有一定的相位差。

将式(1、2)中的每一项自激力和式(7~10)分别带入式(3~6)以考察该项是否具有气动阻尼效应或刚度效应。为了便于说明，下面以A₂^*项在颤振临界状态时引起的自激扭矩在一个周期内作用效应为例具体说明甄别过程，将该项分别带入式(5)和式(6)并结合式(9、10)得到:

$ \begin{aligned} \left.D_{a}\right|_{A_{2}^{*}} &=\rho B^{4} \omega A_{2}^{*} \int_{0}^{2 \pi} \dot{\alpha} \cdot \dot{\alpha} \mathrm{d} \varphi \\ &=\rho B^{4} \omega A_{2}^{*} \int_{0}^{2 \pi}\left(a_{a} \omega \sin \varphi\right)^{2} \mathrm{d} \varphi \\ &=\frac{\rho B^{4}}{2} \omega^{3} a_{a}^{2} A_{2}^{*} \int_{0}^{2 \pi}(1-\cos 2 \varphi) \mathrm{d} \varphi \\ &=\frac{\rho B^{4}}{2} \omega^{3} a_{a}^{2} A_{2}^{*} \varphi \\ & \neq 0 \end{aligned} $

(11)

$ \begin{aligned} \left.K_{a}\right|_{A_{2}^{*}} &=\rho B^{4} \omega A_{2}^{*} \int_{0}^{2 \pi} \dot{\alpha} \cdot a_{a} \mathrm{d} \varphi \\ &=\rho B^{4} \omega A_{2}^{*} \int_{0}^{2 \pi} a_{a} \omega \sin \varphi \cdot a_{a} \cos \varphi \mathrm{d} \varphi \\ &=\frac{\rho B^{4}}{2} \omega^{2} a_{a}^{2} A_{2}^{*} \int_{0}^{2 \pi} \sin 2 \varphi \mathrm{d} \varphi \\ &=0 \end{aligned} $

(12)

由此可以得出A₂^*项在一个周期内对桥梁系统做功，起到消耗或增加机械能的作用，具有阻尼效应；该项在一个周期内刚度效应的贡献值为零，不具有传递和储存能量的作用。由此可知，A₂^*项为纯阻尼效应项。通过上述甄别过程，依据每一项作用效应的不同。将颤振自激力分为3类，即纯阻尼效应项H₁^*、A₂^*，纯刚度效应项A₃^*、H₄^*和既有刚度效应又有阻尼效应的双重效应项A₁^*、H₂^*、H₃^*、A₄^*。

1.2 颤振自激力分类的合理性验证

为了检验颤振自激力分类的合理性，下面分别求出颤振自激力各效应项的做功时程和无功时程。扭转和竖弯模态气动阻尼效应项做功时程表达式如下:

$ \begin{aligned} W_{\alpha, c}=& \int_{0}^{t} \rho B^{4} \omega A_{1}^{*} \dot{h} \cdot \dot{\alpha} \mathrm{d} \tau+\int_{0}^{t} \rho B^{4} \omega A_{2}^{*} \dot{\alpha} \cdot \dot{\alpha} \mathrm{d} \tau \\ &+\int_{0}^{t} \rho B^{3} \omega^{2} A_{4}^{*} h \cdot \dot{\alpha} \mathrm{d} \tau \end{aligned} $

(13)

$ \begin{aligned} W_{h, c}=& \int_{0}^{t} \rho B^{2} \omega H_{1}^{*} \dot{h} \cdot \dot{h} \mathrm{d} \tau+\int_{0}^{t} \rho B^{3} \omega H_{2}^{*} \dot{\alpha} \cdot \dot{h} \mathrm{d} \tau \\ &+\int_{0}^{t} \rho B^{3} \omega^{2} H_{3}^{*} \alpha \cdot \dot{h} \mathrm{d} \tau \end{aligned} $

(14)

扭转和竖弯模态气动刚度效应项无功时程表达式为:

$ \begin{aligned} W_{a, k}=& \int_{0}^{t} \rho B^{3} \omega A_{1}^{*} \dot{h} \cdot \alpha \mathrm{d} \tau+\int_{0}^{t} \rho B^{4} \omega^{2} A_{3}^{*} \alpha \cdot \alpha \mathrm{d} \tau \\ &+\int_{0}^{t} \rho B^{3} \omega^{2} A_{4}^{*} h \cdot \alpha \mathrm{d} \tau \end{aligned} $

(15)

$ \begin{aligned} W_{h, k}=& \int_{0}^{t} \rho B^{3} \omega H_{2}^{*} \dot{\alpha} \cdot h \mathrm{d} \tau+\int_{0}^{t} \rho B^{3} \omega^{2} H_{3}^{*} \alpha \cdot h \mathrm{d} \tau \\ &+\int_{0}^{t} \rho B^{2} \omega^{2} H_{4}^{*} h \cdot h \mathrm{d} \tau \end{aligned} $

(16)

颤振自激力阻尼效应项做功与系统机械能存在转换关系，因此为了对自激力分类的合理性进行验证，还需求出系统的机械能。扭转和竖弯模态机械能表达式如下:

$ \left\{\begin{array}{l} E_{a}=\frac{1}{2}\left(I_{m} \dot{\alpha}^{2}+K_{a} \alpha^{2}\right) \\ E_{h}=\frac{1}{2}\left(m \dot{h}^{2}+K_{h} h^{2}\right) \end{array}\right. $

(17)

本文在ANSYS中建立理想平板简支梁模型，并以模型的扭转模态为例具体说明验证过程。

该理想平板简支梁长L=300 m，宽B=40 m，两端扭转自由度固定。平板断面竖向和横向弯曲刚度为EI_z=2.1×10⁶ MPa·m⁴, EI_y=1.8×10⁷ MPa·m⁴，扭转刚度GI_t=4.1×10⁵ MPa·m⁴。每延米长度质量m=20 000 kg/m，质量惯矩I_m=4.5×10⁶ kg·m²/m，空气密度ρ=1.25 kg/m。结构模态阻尼比均假设为零。其中桥面主梁采用BEAM4单元模拟，质量惯矩采用MASS21单元模拟，自激力采用MATRIX27单元模拟。模型共有120个节点，237个单元。

为获得理想平板简支梁的振动响应，给予主梁节点单位瞬时速度的初始激励，基于ANSYS瞬态动力学分析求得理想平板简支梁跨中节点在风速U₀=142 m/s时扭转模态的速度和位移时程。依据阻尼效应的物理意义，桥梁系统机械能的改变只与自激力阻尼效应项的做功有关。通过理想平板颤振导数理论解和式(13)求得扭转模态气动阻尼效应项做功时程曲线，并通过式(17)求得扭转模态机械能时程曲线。从图 1中可以看出扭转模态阻尼效应项做功时程与机械能时程吻合良好，从而验证了自激力分类的合理性。

2 颤振驱动机理研究 2.1 颤振自激力的做功时程和无功时程

根据1.1节对颤振自激力的分类，可以看出扭转模态中具有阻尼效应的项为A₁^*、A₂^*、A₄^*，具有刚度效应的项为A₁^*、A₃^*、A₄^*，基于能量等效原理，通过积分运算可得扭转模态各效应项的做功时程和无功时程为:

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{l}} {A_1^c = \int_{{t_1}}^{t2} \rho {B^4}\omega A_1^*\dot h \cdot \dot \alpha {\rm{d}}\tau }\\ {A_2^c = \int_{{t_1}}^{t2} \rho {B^4}\omega A_2^*\dot \alpha \cdot \dot \alpha {\rm{d}}\tau }\\ {A_4^c = \int_{{t_1}}^{t2} \rho {B^3}{\omega ^2}A_4^*h \cdot \dot \alpha {\rm{d}}\tau } \end{array}} \right.\\ \\ \left\{ {\begin{array}{*{20}{l}} {A_1^k = \int_{{t_1}}^{{t_2}} \rho {B^3}\omega A_1^*\dot h \cdot \alpha {\rm{d}}\tau }\\ {A_3^k = \int_{{t_1}}^{t2} \rho {B^4}{\omega ^2}A_3^*\alpha \cdot \alpha {\rm{d}}\tau }\\ {A_4^k = \int_{{t_1}}^{{t_2}} \rho {B^3}{\omega ^2}A_4^*h \cdot \alpha {\rm{d}}\tau } \end{array}} \right. \end{array} $

(18)

竖弯模态中具有阻尼效应的项为H₁^*、H₂^*、H₃^*，具有刚度效应的项为H₂^*、H₃^*、H₄^*，竖弯模态各效应性的做功时程和无功时程为:

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{l}} {H_1^c = \int_0^t \rho {B^2}\omega H_1^*\dot h \cdot \dot h{\rm{d}}\tau }\\ {H_2^c = \int_0^t \rho {B^3}\omega H_2^*\dot \alpha \cdot \dot h{\rm{d}}\tau }\\ {H_3^c = \int_0^t \rho {B^3}{\omega ^2}H_3^*\alpha \cdot \dot h{\rm{d}}\tau } \end{array}} \right.\\ \\ \left\{ {\begin{array}{*{20}{l}} {H_2^k = \int_0^t \rho {B^3}\omega H_2^*\dot \alpha \cdot h{\rm{d}}\tau }\\ {H_3^k = \int_0^t \rho {B^3}{\omega ^2}H_3^*\alpha \cdot h{\rm{d}}\tau }\\ {H_4^k = \int_0^t \rho {B^2}{\omega ^2}H_4^*h \cdot h{\rm{d}}\tau } \end{array}} \right. \end{array} $

(19)

纯阻尼效应项H₁^*、A₂^*的无功时程和纯刚度效应项A₃^*、H₄^*的做功时程为:

$ \begin{array}{l} \left\{ {\begin{array}{*{20}{l}} {H_1^k = \int_0^t \rho {B^2}\omega H_1^*\dot h \cdot h{\rm{d}}\tau }\\ {A_2^k = \int_0^t \rho {B^3}\omega H_2^*\dot \alpha \cdot \alpha {\rm{d}}\tau } \end{array}} \right.\\ \\ \left\{ {\begin{array}{*{20}{l}} {A_3^c = \int_0^t \rho {B^4}{\omega ^2}A_3^*\alpha \cdot \dot \alpha {\rm{d}}\tau }\\ {H_4^c = \int_0^t \rho {B^2}{\omega ^2}H_4^*h \cdot \dot h{\rm{d}}\tau } \end{array}} \right. \end{array} $

(20)

2.2 从功能关系角度研究颤振驱动机理

以理想平板简支梁为例来研究桥梁颤振驱动机理，将理想平板简支梁跨中节点的速度和位移响应时程以及理想平板颤振导数理论解带入式(18~20)，求得各效应项在颤振前、颤振临界状态和颤振后做功时程和无功时程，如图 2和图 3所示。

从图 2中可以看出纯刚度效应项A₃^*、H₄^*在扭转和竖弯模态的整个做功时程内做功为零，说明了纯刚度效应项不会增加或耗散系统的机械能；纯阻尼效应项H₁^*、A₂^*分别在竖弯和扭转模态的整个做功时程内做负功，耗散系统的机械能，具有维持系统稳定的作用；双重效应项A₁^*、H₃^*在竖弯和扭转模态的整个做功时程中做正功，随着风速的增加该项所做的增加系统机械能的功将逐渐大于H₁^*、A₂^*所做的减少机械能的功，并最终导致颤振发生。从图 3中可以看出纯阻尼效应项H₁^*、A₂^*在竖弯和扭转模态的整个无功时程中做功为零，说明了纯阻尼效应项不会改变系统相位；纯刚度项A₃^*在扭转模态的整个无功时程中做正功，双重效应项H₃^*在竖弯模态的无功时程中做负功，这种作用效果使得扭转和竖弯模态具有一定的相位差，为颤振的发生提供了条件。

3 基于瞬时做功的颤振导数识别方法

Scanlan颤振自激力模型中包含了纯刚度效应项、纯阻尼效应项和既有刚度效应又有阻尼效应的双重效应项，这些项对自激振动响应的影响不尽相同。依据能量等效原理，气动阻尼项只耗散(或增加)系统的机械能，而气动刚度项只影响振动的相位，不会影响机械能的大小。因而，若从基于自激扭矩和自激升力的做功时程来识别气动阻尼参数，可以提高气动参数的识别精度。相似地，若基于自激扭矩和自激升力的无功时程来识别气动刚度参数，也可以提高气动刚度参数的识别精度。与传统识别方法相比，基于瞬时做功的颤振导数识别方法结算结果更加稳定。

3.1 将颤振微分方程转化为功能方程

本文在建立颤振导数识别方法时，采用自激力做功的形式，通过颤振微分方程两边同时积分将运动方程转化为功能方程的形式，分别求出了自激扭矩和自激升力的做功时程和无功时程，由此建立颤振导数与系统振动响应之间的关系，然后通过加权最小二乘法求得颤振导数。通过积分运算得到扭转颤振自激力做功时程，如式(21)所示:

$ \begin{aligned} W_{a, c}(t) &=\int_{t_{1}}^{t} M_{\mathrm{se}}(\tau) \dot{\alpha}(\tau) \mathrm{d} \tau \\ &=\int_{t_{1}}^{t}\left[I_{m} \ddot{\alpha}+2 I_{m} \xi_{\alpha} \omega_{\alpha} \dot{\alpha}+I_{m} \omega_{a}^{2} \alpha\right] \cdot \dot{\alpha} \mathrm{d} \tau \\ &=\int_{t_{1}}^{t} I_{m} \dot{\alpha} \mathrm{d} \dot{\alpha}+\int_{t_{1}}^{t} k_{\alpha} \alpha \mathrm{d} \alpha+\int_{t_{1}}^{t} c_{\alpha} \dot{\alpha} \dot{\alpha} \mathrm{d} \tau \\ &=\int_{t_{1}}^{t} \mathrm{d}\left[\frac{1}{2} I_{m} \dot{\alpha}^{2}+\frac{1}{2} k_{\alpha} \alpha^{2}\right]+W_{\alpha}^{c}(t) \\ &=E_{a}(t)-E_{\alpha}\left(t_{1}\right)+W_{\alpha}^{c}(t) \end{aligned} $

(21)

竖弯颤振自激力做功时程，如式(22)所示:

$ \begin{aligned} W_{h, c}(t) &=\int_{t_{1}}^{t} L_{\mathrm{se}}(\tau) \cdot \dot{h}(\tau) \mathrm{d} \tau \\ &=\int_{t_{1}}^{t}\left[m_{h} \ddot{h}+2 m_{h} \xi_{h} \omega_{h} \dot{h}+m_{h} \omega_{h}^{2} h\right] \cdot \dot{h} \mathrm{d} \tau \\ &=\int_{t_{1}}^{t} m \dot{h} \mathrm{d} \dot{h}+\int_{t_{1}}^{t} k_{h} h \mathrm{d} h+\int_{t_{1}}^{t} c_{h} \dot{h} \dot{h} \mathrm{d} \tau \\ &=\int_{t_{1}}^{t} \mathrm{d}\left[\frac{1}{2} m \dot{h}^{2}+\frac{1}{2} k_{h} h^{2}\right]+\int_{t_{1}}^{t} c_{h} \dot{h} \dot{h} \mathrm{d} \tau \\ &=E_{h}(t)-E_{h}\left(t_{1}\right)+W_{h}^{c}(t) \end{aligned} $

(22)

式中:E_α(t)、E_h(t)分别表示扭转和竖弯模态的瞬时机械能，W_α^c(t)、W_h^c(t)分别表示系统扭转和竖弯模态阻尼力所做的功。式(21、22)表明颤振自激力瞬时做功时程等于机械能变化量与系统阻尼力做功变化量之和。类似地，通过积分得到扭转颤振自激力无功时程，如式(23)所示:

$ \begin{aligned} W_{a, k} &=\int_{t_{1}}^{t} M_{\mathrm{se}}(\tau) \cdot \alpha(\tau) \mathrm{d} \tau \\ &=\int_{t_{1}}^{t}\left[I_{m} \ddot{\alpha}+2 I_{m} \xi_{a} \omega_{a} \dot{\alpha}+I_{m} \omega_{a}^{2} \alpha\right] \cdot \alpha \mathrm{d} \tau \\ &=\int_{t_{1}}^{t} I_{m} \ddot{\alpha} \alpha \mathrm{d} \tau+\int_{t_{1}}^{t} I_{m} \omega_{a}^{2} \alpha^{2} \mathrm{d} \tau+2 \int_{t_{1}}^{t} I_{m} \xi_{a} \omega_{a} \dot{\alpha} \alpha \mathrm{d} \tau \\ &=-\int_{t_{1}}^{t} I_{m} \omega_{t}^{2}(t) \alpha^{2} \mathrm{d} \tau+\int_{t_{1}}^{t} I_{m} \omega_{\alpha}^{2}(t) \alpha^{2} \mathrm{d} \tau+0 \\ &=I_{m} \int_{t_{0}}^{t}\left[\omega_{a}^{2}(t)-\omega_{t}^{2}(t)\right] \alpha^{2} \mathrm{d} \tau \end{aligned} $

(23)

竖弯颤振自激力无功时程，如式(24)所示:

$ \begin{aligned} W_{h, k} &=\int_{t_{1}}^{t} L_{\mathrm{se}}(\tau) \cdot h(\tau) \mathrm{d} \tau \\ &=\int_{t_{1}}^{t}\left[m \ddot{h}+2 m \xi_{h} \omega_{h} \dot{h}+m \omega_{h}^{2} h\right] \cdot h \mathrm{d} \tau \\ &=\int_{t_{1}}^{t} m \ddot{h} h \mathrm{d} \tau+\int_{t_{1}}^{t} m \omega_{h}^{2} h^{2} \mathrm{d} \tau+2 \int_{t_{1}}^{t} c_{h} \dot{h} h \mathrm{d} \tau \\ &=-\int_{t_{1}}^{t} m \omega_{t}^{2}(t) h^{2} \mathrm{d} \tau+\int_{t_{1}}^{t} m \omega_{h}^{2}(t) h^{2} \mathrm{d} \tau+0 \\ &=m \int_{t_{1}}^{t}\left[\omega_{h}^{2}(t)-\omega_{t}^{2}(t)\right] h^{2} \mathrm{d} \tau \end{aligned} $

(24)

由于竖弯和扭转模态的位移和速度近似正交，式中系统阻尼力的无功为零。式(23、24)表明颤振自激力的无功时程等于刚度恢复力的无功时程之差。

3.2 加权最小二乘法识别颤振导数

本节以扭转模态为例具体说明颤振导数的识别过程。由式(13)、式(15)并结合式(21)、式(23)可得到颤振微分方程两边同时积分后的功能方程，如式(25)所示:

$ \left\{ {\begin{array}{*{20}{c}} {{E_a}(t) - {E_\alpha }\left( {{t_1}} \right) + W_\alpha ^c(t) = }\\ {\int_{{t_1}}^t \rho {B^4}\omega A_1^*\dot h \cdot \dot \alpha {\rm{d}}\tau + }\\ {\int_{{t_1}}^t \rho {B^4}\omega A_2^*\dot \alpha \cdot \dot \alpha {\rm{d}}\tau + }\\ {\int_{{t_1}}^t \rho {B^3}{\omega ^2}A_4^*h \cdot \dot \alpha {\rm{d}}\tau }\\ {{I_m}\int_{{t_1}}^t {\left[ {\omega _t^2(t) - \omega _a^2(t)} \right]} {\alpha ^2}{\rm{d}}\tau = }\\ {\int_{{t_1}}^t \rho {B^3}\omega A_1^*\dot h \cdot \alpha {\rm{d}}\tau + }\\ {\int_{{t_1}}^t \rho {B^4}{\omega ^2}A_3^*\alpha \cdot \alpha {\rm{d}}\tau + }\\ {\int_{{t_1}}^t \rho {B^3}{\omega ^2}A_4^*h \cdot \alpha {\rm{d}}\tau } \end{array}} \right. $

(25)

将上式写为矩阵形式:

$ \left\{\begin{array}{l} W_{a \cdot c}=W_{m \cdot c} A_{a \cdot c} \\ W_{a, k}=W_{m, k} A_{a, k} \end{array}\right. $

(26)

其中

$ \begin{array}{l} W_{a . c}= \\ {\left[\begin{array}{c} E\left(t_{2}\right)-E_{h}\left(t_{1}\right)+W_{a}^{c}\left(t_{2}\right)-W_{a}^{c}\left(t_{1}\right) \\ E\left(t_{3}\right)-E_{h}\left(t_{2}\right)+W_{a}^{c}\left(t_{3}\right)-W_{a}^{c}\left(t_{2}\right) \\ \vdots \\ E\left(t_{n}\right)-E_{h}\left(t_{n-1}\right)+W_{a}^{c}\left(t_{n}\right)-W_{a}^{c}\left(t_{n-1}\right) \end{array}\right]} \end{array} $

(27)

$ W_{a, k}=\left[\begin{array}{c} I_{m} \int_{t_{1}}^{t 2}\left[\omega_{t}^{2}(t)-\omega_{a}^{2}(t)\right] \alpha^{2} \mathrm{d} \tau \\ I_{m} \int_{t_{2}}^{t_{3}}\left[\omega_{t}^{2}(t)-\omega_{\alpha}^{2}(t)\right] \alpha^{2} \mathrm{d} \tau \\ \vdots \\ I_{m} \int_{t_{n-1}}^{t_{n}}\left[\omega_{t}^{2}(t)-\omega_{a}^{2}(t)\right] \alpha^{2} \mathrm{d} \tau \end{array}\right] $

(28)

$ \begin{array}{l} {W_{m,c}} = \\ \rho \omega {B^3}\left[ \begin{array}{l} \begin{array}{*{20}{c}} {\int_{{t_1}}^{t2} B \dot h\dot \alpha {\rm{d}}\tau }&{\int_{{t_1}}^{t2} B \dot \alpha \dot \alpha {\rm{d}}\tau }&{\int_{{t_1}}^{t2} \omega h\dot \alpha {\rm{d}}\tau }\\ {\int_{t2}^{t3} B \dot h\dot \alpha {\rm{d}}\tau }&{\int_{t2}^{t3} B \dot \alpha \dot \alpha {\rm{d}}\tau }&{\int_{t2}^{t3} \omega h\dot \alpha {\rm{d}}\tau }\\ {}& \vdots &{} \end{array}\\ \int_{{t_{n - 1}}}^{{t_n}} B \dot h\dot \alpha {\rm{d}}\tau \;\;\;\int_{{t_{n - 1}}}^{{t_n}} B \dot \alpha \dot \alpha {\rm{d}}\tau \;\;\;\;\int_{{t_{n - 1}}}^{{t_n}} {\omega h\dot \alpha {\rm{d}}\tau } \; \end{array} \right] \end{array} $

(29)

$ \begin{array}{l} {W_{m,k}} = \\ \rho \omega {B^3}\left[ {\begin{array}{*{20}{l}} {\int_{{t_1}}^{t2} {\dot h} \alpha {\rm{d}}\tau }&{\int_{{t_1}}^{t2} B \omega {\alpha ^2}{\rm{d}}\tau \;\;\;\int_{{t_1}}^{t2} {\omega h\alpha {\rm{d}}\tau } }\\ {\int_{{t_2}}^{{t_3}} {\dot h} \alpha {\rm{d}}\tau }&{\int_{{t_2}}^{t3} B \omega {\alpha ^2}{\rm{d}}\tau \;\;\;{\kern 1pt} \int_{{t_1}}^{t2} {B\omega } h{\alpha ^2}{\rm{d}}\tau }\\ {}&{\;\;\;\;\;\;\; \vdots }\\ {\int_{{t_{n - 1}}}^{{t_1}} {\dot h} \alpha {\rm{d}}\tau }&{\int_{{t_{n - 1}}}^{{t_1}} B \omega {\alpha ^2}{\rm{d}}\tau \;\;\;{\kern 1pt} \int_{{t_{n - 1}}}^{{t_n}} B \omega {\alpha ^2}{\rm{d}}\tau } \end{array}\;} \right] \end{array} $

(30)

$ \begin{array}{l} A_{\alpha, c}=\left[\begin{array}{lll} A_{1}^{*} & A_{2}^{*} & A_{4}^{*} \end{array}\right]^{\mathrm{T}} \\ A_{\alpha, k}=\left[\begin{array}{lll} A_{1}^{*} & A_{3}^{*} & A_{4}^{*} \end{array}\right]^{\mathrm{T}} \end{array} $

(31)

其中n为颤振过程的采样点数，W_{α, c}和W_{α, k}表示通过振动响应计算得到的扭转自激力的做功矩阵和无功矩阵，W_{m, c}和W_{m, k}为扭转自激力各效应项瞬时的做功矩阵和无功矩阵，A_{α, c}为气动阻尼效应项矩阵，A_{α, k}为气动刚度效应项矩阵。

在式(26)中，颤振导数估计的误差分别为:

$ \left\{\begin{array}{l} e_{\alpha, c}=W_{\alpha, c}-W_{m, c} A_{\alpha, c} \\ e_{\alpha, k}=W_{\alpha, k}-W_{m, k} A_{\alpha, k} \end{array}\right. $

(32)

误差函数可采用如下加权形式：

$ \left\{\begin{aligned} J_{a \cdot c} &=w_{g}^{2} e_{a, c}^{\mathrm{T}} e_{a, c} \\ &=w_{g}^{2}\left(W_{a, c}-W_{m, c} A_{a \cdot c}\right)^{\mathrm{T}}\left(W_{a \cdot c}-W_{m \cdot c} A_{a \cdot c}\right) \\ J_{a \cdot k} &=w_{s}^{2} e_{a \cdot k}^{\mathrm{T}} e_{a \cdot k} \\ &=w_{s}^{2}\left(W_{a, k}-W_{m, k} A_{a, k}\right)^{\mathrm{T}}\left(W_{a, k}-W_{m, k} A_{a, k}\right) \end{aligned}\right. $

(33)

其中w_g²和w_s²分别为做功时程和无功时程的权重，该加权矩阵建议取自激力在做功和无功时程所做的功。为了使总体误差函数取得最小值，需满足

$ \left\{\begin{array}{l} \frac{\partial J_{a, c}}{\partial A_{a, c}}=2 w_{g}^{2}\left(-W_{m, c}^{\mathrm{T}} W_{a, c}+W_{m, c}^{\mathrm{T}} W_{m, c} A_{a, c}\right)=0 \\ \frac{\partial J_{a, k}}{\partial A_{a, k}}=2 w_{s}^{2}\left(-W_{m, k}^{\mathrm{T}} W_{a, k}+W_{m, k}^{\mathrm{T}} W_{m, k} A_{\alpha, k}\right)=0 \end{array}\right. $

(34)

由此可得颤振导数的最小二乘估计：

$ \left\{\begin{array}{l} A_{a, c}=\left(W_{m, c}^{\mathrm{T}} w_{g}^{2} W_{m, c}\right)^{-1} W_{m, c}^{\mathrm{T}} x w_{g}^{2} W_{a, c} \\ A_{a, k}=\left(W_{m, k}^{\mathrm{T}} w_{s}^{2} W_{m, k}\right)^{-1} W_{m, k}^{\mathrm{T}} w_{s}^{2} W_{a, k} \end{array}\right. $

(35)

通过上述建立的基于自激力瞬时做功的颤振导数识别方法求得理想平板与扭转振动有关的颤振导数如图 4所示。通过对比发现，本文解与理论解吻合良好，证明了上述颤振导数识别方法的可靠性。

4 结论

(1) 基于能量等效原理甄别出颤振自激力中只起到阻尼效应的纯阻尼项H₁^*、A₂^*；只起到刚度效应的纯刚度项A₃^*、H₄^*和既有刚度效应又有阻尼效应的双重效应项A₁^*、H₂^*、H₃^*、A₄^*。通过阻尼效应项做功时程与系统机械能时程的对比分析发现这种按作用效应对颤振自激力进行分类的方法是合理的。

(2) 纯阻尼项H₁^*、A₂^*在做功时程内做负功，是竖弯和扭转振动主要耗能因素；双重效应项H₃^*、A₁^*在做功时程内做正功，是竖弯和扭转振动主要的增能因素。随着风速的增大，颤振自激力的增能能力将逐渐大于耗能能力，并导致颤振的发生。

(3) 以理想平板为例，将基于瞬时做功的颤振导数识别方法的结果与颤振导数理论解对比发现，该方法所识别的颤振导数具有较好的可靠性，可以考虑将该方法应用到二维节段模型自由振动试验的颤振导数识别中。