CN103218481A - 大跨度桥梁风致灾变全过程的模拟方法 - Google Patents

Abstract

本发明提供了一种大跨度桥梁风致灾变全过程模拟方法，该方法将三分力系数、有效瞬时风攻角和相对风速统一作用在主梁的气动力上，然后对其中的自激力成分进行非定常修正；将气动力荷载施加到有限元模型上进行动力时程计算，并在各级风速下进行扭转基频和竖弯基频的迭代；根据位移响应，判断桥梁结构在某级风速下的稳定状态。本发明实现了静风失稳、抖振和颤振的统一计算，克服了传统风振计算中单独计算的不足，提高了桥梁风振计算的精度。

Description

大跨度桥梁风致灾变全过程的模拟方法

技术领域

本发明属于桥梁抗风设计领域，尤其是一种考虑了各种风振形式耦合作用的大跨度桥梁风致灾变全过程的模拟方法。

背景技术

随着科学技术的进步和设计水平的提高，21世纪的世界桥梁工程已经进入了跨海连岛建设的新时期，现代桥梁不断向着长大化方向发展。桥梁跨径大幅度增长带来的主要问题是结构刚度的急剧下降，使得桥梁结构对风作用比较敏感，因此桥梁风致响应问题变得越来越重要。

桥梁风致灾变全过程模拟旨在用一个统一的数学模型解决风速从小到大的所有结构风振形式模拟问题，具体就是低风速时的抖振、极限风速时的静风失稳和颤振。静风失稳、颤振和抖振是同一事物的三个方面，它们是桥梁结构风振响应中人为剥离出来的不同形态。在静风力、自激力、抖振力甚至涡激力的共同作用下，桥梁结构产生一般性的风振响应。静风失稳和颤振是均匀流场中的特例，抖振是紊流场中的振动响应。传统的风致响应计算中，将静风失稳、颤振和抖振单独计算，忽略了它们之间的耦合作用，难以适应超大跨度桥梁风振分析的需要。

传统的抗风理论中，将来流风速分解为平均风速和脉动风速，将气动力分解为静风荷载、抖振力和自激力三部分，然后叠加在桥梁结构上。其中，静风荷载按常规静力三分力系数计算；抖振力通常采用Scanlan的准定常理论计算,并使用Davenport引入的气动导纳函数进行修正；自激力通常采用Scanlan提出的气动导数的线性表达式，并利用Lin提出的脉冲响应函数方法将自激力进行时域化表达。上述叠加气动力模型存在以下不足：在紊流风作用下，主梁静风力、抖振力和自激力三项的效应实际上是相互影响的，严格来说是不能截然分开再进行叠加的。

发明内容

发明目的：为克服传统桥梁风振计算方法的上述两方面不足，本发明提供了一种同时考虑静风失稳、随机抖振和颤振发散之间相互耦合作用的桥梁结构风致灾变的全过程模拟方法。

发明内容：本发明的大跨度桥梁风致灾变全过程模拟方法包括以下步骤：

第一步：建立桥梁有限元模型，确定桥梁结构在恒载作用下的初始状态；

第二步：进行无阻尼的结构模态分析，提取其一阶竖弯频率f_v0和一阶扭转频率f_t0；

第三步：给定初始平均风速U₁和风速增量ΔU；

第四步：设当前平均风速U_i，采用谐波合成法生成水平和垂直脉动风时程；假定系统的竖弯振动频率f_vi和扭转振动频率f_ti分别为上级风速下的竖弯振动频率f_vi-1和扭转振动频率f_ti-1，当i=1时即为桥梁结构的一阶竖弯频率f_v0和一阶扭转频率f_t0；

第五步：根据初始风攻角和各单元的扭转变形均值计算各单元的有效攻角；然后根据当前风速、试算振动频率和有效攻角二次插值计算各单元的颤振导数，即先在不同初始攻角下相应的折减风速处进行第一次插值计算，而后根据实际的有效攻角进行第二次插值计算；

第六步：计算统一气动力荷载，进行动力时程分析，具体如下列公式表示：

F_z(t)＝D(t)cos(α₀+ψ)-L(t)sin(α₀+ψ)    (1)

$F_y\left(t\right)=\frac{1}{2}\mathrm{\ρ}{V}_r^{2}B(C_{\mathrm{Dw}}\sin {\mathrm{\ψ}}_w+C_{\mathrm{Lw}}\cos {\mathrm{\ψ}}_w)+\frac{1}{2}\mathrm{\ρ}{V}_r^{2}B(h_1{H}_1\frac{\stackrel{\·}{h}}{U_r}+h_2{H}_2\frac{B\stackrel{\·}{\mathrm{\θ}}}{U_r}+h_3{H}_3\mathrm{\θ})---\left(2\right)$

$M_x\left(t\right)=\frac{1}{2}\mathrm{\ρ}{V}_r^2{B}^2{C}_{\mathrm{Mw}}+\frac{1}{2}\mathrm{\ρ}{V}_r^2{B}^2(a_1{A}_1\frac{\stackrel{\·}{h}}{U_r}+a_2{A}_2\frac{B\stackrel{\·}{\mathrm{\θ}}}{U_r}+a_3{A}_3\mathrm{\θ})---\left(3\right)$

其中，

$D\left(t\right)=\frac{1}{2}\mathrm{\ρ}{V}_r^2\left(t\right)B{C}_D\left[{\mathrm{\α}}_e\left(t\right)\right]---\left(4\right)$

$L\left(t\right)=\frac{1}{2}\mathrm{\ρ}{V}_r^2\left(t\right)B{C}_L\left[{\mathrm{\α}}_e\left(t\right)\right]---\left(5\right)$

α_e(t)＝α₀+ψ(t)+θ(t)    (6)

$V_r\left(t\right)={\{U_r^2\left(t\right)+{[w\left(t\right)-\stackrel{\·}{h}-r\stackrel{\·}{\mathrm{\θ}}]}^2\}}^{1/2}---\left(7\right)$

$\mathrm{\ψ}\left(t\right)=[w\left(t\right)-\stackrel{\·}{h}-r\stackrel{\·}{\mathrm{\θ}}]/U_r\left(t\right)---\left(8\right)$

$U_r\left(t\right)=U+u\left(t\right)-\stackrel{\·}{p}---\left(9\right)$

$h_1=\frac{2}{(C_{L0}-K_{D0})\sin {\mathrm{\α}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\α}}_0}K{H}_1^{*}---\left(10\right)$

$h_2=\frac{2B}{[(C_{L0}-K_{D0})\sin {\mathrm{\α}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\α}}_0]r}K{H}_2^{*}---\left(11\right)$

$h_3=\frac{2}{K_{D0}\sin {\mathrm{\α}}_0+K_{L0}\cos {\mathrm{\α}}_0}{K}^2{H}_3^{*}---\left(12\right)$

$a_1=\frac{-2}{K_{M0}}K{A}_1^{*},a_2=\frac{-2B}{K_{M0}r}K{A}_2^{*},a_3=\frac{2}{K_{M0}}{K}^2{A}_3^{*}---\left(13\right)$

$C_{D0}=C_D\left({\mathrm{\α}}_0\right),C_{L0}=C_L\left({\mathrm{\α}}_0\right),K_{D0}=\frac{d{C}_D}{\mathrm{d\α}}{|}_{\mathrm{\α}={\mathrm{\α}}_0},K_{L0}={\frac{d{C}_L}{\mathrm{d\α}}|}_{\mathrm{\α}={\mathrm{\α}}_0},K_{M0}={\frac{d{C}_M}{\mathrm{d\α}}|}_{\mathrm{\α}={\mathrm{\α}}_0}---\left(14\right)$

ψ_w＝α₀+w/U_r    (15)

$C_{\mathrm{Dw}}=C_D\left({\mathrm{\ψ}}_w\right),C_{\mathrm{Lw}}=C_L\left({\mathrm{\ψ}}_w\right),C_{\mathrm{Mw}}=C_M\left({\mathrm{\ψ}}_w\right),K_{\mathrm{Dw}}={\frac{d{C}_D}{\mathrm{d\α}}|}_{\mathrm{\α}={\mathrm{\ψ}}_w},K_{\mathrm{Lw}}={\frac{d{C}_L}{\mathrm{d\α}}|}_{\mathrm{\α}={\mathrm{\ψ}}_w},K_{\mathrm{Mw}}={\frac{d{C}_M}{\mathrm{d\α}}|}_{\mathrm{\α}={\mathrm{\ψ}}_w}---\left(16\right)$

H₁=(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w          (17)

H₂=[(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w]r/B     (18)

H₃=K_Dwsinψ_w+K_Lwcosψ_w                        (19)

A₁=－K_Mw,A₂=－K_Mwr/B,A₃=K_Mw                   (20)

在上述表达式中，F_z(t)、F_y(t)和M_x(t)分别为主梁单位长度上任一时刻在总体坐标系上所受的气动阻力、升力和升力矩；D(t)和L(t)分别为主梁单位长度上任一时刻在相对风轴坐标系上的气动阻力和升力；ρ为空气密度；B为主梁宽度；U为主梁高度处的平均风速，u(t)、w(t)分别为顺风向和竖向脉动风速，U_r(t)为来流风轴合成速度，V_r(t)为相对风速；α₀为初始风攻角，α_e(t)为瞬时风攻角，ψ(t)为来流风轴坐标系内的瞬时风攻角分量，ψ_w为总体坐标系内竖向脉动风速产生的瞬时风攻角分量；p、h和θ分别为主梁振动的水平、竖向和扭转位移，、

和

分别为主梁竖向、水平和扭转振动速度；C_D、C_L和C_M是主梁风轴静力三分力系数（参考长度均为主梁宽度B），分别是阻力系数、升力系数和升力矩系数，它们是风攻角的函数；C_D0和C_L0分别为初始风攻角α₀对应的阻力系数和升力系数；K_D0、K_L0和K_M0分别为初始风攻角α₀对应的阻力系数导数、升力系数导数和升力矩系数导数；C_Dw、C_Lw和C_Mw分别ψ_w对应的阻力系数、升力系数和升力矩系数；K_Dw、K_Lw和K_Mw分别ψ_w对应的阻力系数导数、升力系数导数和升力矩系数导数；r为主梁特征长度；K=ωB/U，为折减频率，ω为振动圆频率；

和

是颤振导数，均为折减频率K的函数，与桥梁断面的几何构形和来流有关；h_i和a_i(i=1,2,3)是和折减频率相关的修正系数；H_i和A_i(i=1,2,3)是准气动导数。

第七步：根据各级风速的位移响应结果评价抖振性能和判断静风失稳：在较低风速下，通过均值、根方差和功率谱来评价结构的抖振性能；以位移响应均值的出现峰值（或者说是突然增大）来判断静风失稳。

其中，r是反映气动力作用点与截面形心之间距离的系数，一般取B/4。

所述第六步进一步包括：对竖弯和扭转位移响应进行频谱分析，确定竖弯振动卓越频率f_v和扭转振动卓越频率f_t。将f_v和f_t与试算值f_vi和f_ti分别进行比较，如果|(f_v-f_vi)/f_vi|<0.05且|(f_t–f_ti)/f_ti|<0.05，说明频率迭代收敛。根据位移响应判断是否出现了颤振发散，随着风速的增长，如果振动形式逐渐从随机振动过渡到谐波发散振动，振幅逐渐增大,相应振动的阻尼将逐渐减小到0时就是颤振。如果出现了颤振发散，则结束计算；否则，增加风速，进入第四步；如果|(f_v-f_vi)/f_vi|>0.05或|(f_t–f_ti)/f_ti|>0.05，则计算各单元的扭转变形均值，并令f_vi=f_v，f_ti=f_t，重复第五步和第六步。

有益效果：在统一气动力模型中，综合考虑了静风荷载、抖振力和自激力，不需要分别单独列式再对各自效应进行叠加，而且计入了脉动风、结构位移和结构速度等的高阶项，这更符合气动力非线性的特点；由于该模型实现了静风失稳、抖振和颤振的统一计算，考虑了它们之间的耦合作用，提高了桥梁风振计算的精度。其次，该模型通过有效攻角和相对风速的概念来考虑自激作用，能更直观地体现自激力的产生机理。最后，相对于叠加气动力模型，该气动力模型仍保留了列式简单直观、物理意义明确的优点。

附图说明

图1为本发明相对风速和有效瞬时风攻角示意图，其中，附图标记①为总体坐标系，②为来流风轴坐标系，③为相对风轴坐标系；

图2为本发明风致灾变过程模拟程序流程图。

具体实施方式：

下面通过附图和实施例，对本发明技术方案进行详细说明，但是本发明的保护范围不局限于这些实施例。

如图1所示，本发明首先介绍统一气动力模型，将桥梁结构振动速度等效成来流风速，即将主梁的侧向振动速度等效成来流顺风向的风速，将主梁的竖向振动速度等效成来流垂直向的风速；顺风向风速和垂直向风速合成为相对风速，其与运动中主梁的夹角称为有效瞬时风攻角，包含三部分：初始风攻角α₀、风速的垂直成分产生的攻角ψ，以及主梁的扭转角θ；然后利用三分力系数、有效瞬时风攻角和相对风速统一表达作用在主梁上的气动力。具体而言，基于准定常假定，主梁单位长度上任一时刻所受的气动阻力、升力和升力矩在总体坐标系上可分别表示如下：

F_z(t)＝D(t)cos(α₀+ψ)-L(t)sin(α₀+ψ)    (1a)

F_y(t)＝D(t)sin(α₀+ψ)+L(t)cos(α₀+ψ)    (1b)

M_x(t)＝M(t)                               (1c)

式中，α₀为初始风攻角；ψ为来流风轴坐标系内的瞬时风攻角分量；D(t)、L(t)和M(t)分别为主梁单位长度上任一时刻在相对风轴坐标系上的气动阻力和升力，它们可表示为：

$D\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B{C}_D\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(2a\right)$

$L\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B{C}_L\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(2b\right)$

$M\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B^2{C}_M\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(2c\right)$

式中，ρ为空气密度；B为主梁宽度；C_D、C_L和C_M是主梁风轴静力三分力系数（参考长度均为主梁宽度B），分别是阻力系数、升力系数和升力矩系数，它们是风攻角的函数；α_e(t)和V_r(t)分别为瞬时风攻角和相对风速，它们可以表示为：

α_e(t)＝α₀+ψ(t)+θ(t)    (3)

$V_r\left(t\right)={\{U_r^2\left(t\right)+{[w\left(t\right)-\stackrel{\&CenterDot;}{h}-r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}]}^2\}}^{1/2}---\left(4\right)$

其中，

$\mathrm{\&psi;}\left(t\right)=[w\left(t\right)-\stackrel{\&CenterDot;}{h}-r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}]/U_r\left(t\right)---\left(5\right)$

$U_r\left(t\right)=U+u\left(t\right)-\stackrel{\&CenterDot;}{p}---\left(6\right)$

式中，α₀为初始风攻角；ψ(t)为来流风轴坐标系内的瞬时风攻角分量；θ为任一时刻主梁振动的扭转位移，

、

和

分别为主梁竖向、水平和扭转振动速度；r为主梁特征长度，反映气动力作用点与截面形心之间距离的系数，取B/4；U为主梁高度处的平均风速，u(t)、w(t)分别为顺风向和竖向脉动风速，U_r(t)为来流风轴合成速度。

由于主梁振动所引起的攻角改变相对较小，所以式(2)中的阻力D、升力L和升力矩M可以展开成泰勒级数并截取线性项如下：

$D\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(C_{\mathrm{Dw}}-K_{\mathrm{Dw}}\frac{\stackrel{\&CenterDot;}{h}}{U_r}-K_{\mathrm{Dw}}\frac{r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+K_{\mathrm{Dw}}\mathrm{\&theta;})---\left(7a\right)$

$L\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(C_{\mathrm{Lw}}-K_{\mathrm{Lw}}\frac{\stackrel{\&CenterDot;}{h}}{U_r}-K_{\mathrm{Lw}}\frac{r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+K_{\mathrm{Lw}}\mathrm{\&theta;})---\left(7b\right)$

$M\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2(C_{\mathrm{Mw}}-K_{\mathrm{Mw}}\frac{\stackrel{\&CenterDot;}{h}}{U_r}-K_{\mathrm{Mw}}\frac{r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+K_{\mathrm{Mw}}\mathrm{\&theta;})---\left(7c\right)$

$\sin ({\mathrm{\&alpha;}}_0+\mathrm{\&psi;})=\sin {\mathrm{\&psi;}}_w-(\frac{\stackrel{\&CenterDot;}{h}}{U_r}+\frac{r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r})\cos {\mathrm{\&psi;}}_w---\left(8a\right)$

$\cos ({\mathrm{\&alpha;}}_0+\mathrm{\&psi;})=\cos {\mathrm{\&psi;}}_w-(\frac{\stackrel{\&CenterDot;}{h}}{U_r}+\frac{r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r})\sin {\mathrm{\&psi;}}_w---\left(8b\right)$

其中，

ψ_w＝α₀+w/U_r,C_Dw＝C_D(ψ_w),C_Lw＝C_L(ψ_w),C_Mw＝C_M(ψ_w),

$K_{\mathrm{Dw}}={\frac{d{C}_D}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w},K_{\mathrm{Lw}}={\frac{d{C}_L}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w},K_{\mathrm{Mw}}={\frac{d{C}_M}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w}---\left(9\right)$

式中，ψ_w为总体坐标系内竖向脉动风速产生的瞬时风攻角分量；C_Dw、C_Lw和C_Mw分别ψ_w对应的阻力系数、升力系数和升力矩系数；K_Dw、K_Lw和K_Mw分别ψ_w对应的阻力系数导数、升力系数导数和升力矩系数导数；其他符号意义同前。

将式(7)和式(8)带入式(1)，并忽略高阶项，可得：

$F_y\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(C_{\mathrm{Dw}}\sin {\mathrm{\&psi;}}_w+C_{\mathrm{Lw}}\cos {\mathrm{\&psi;}}_w)+\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(H_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+H_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+H_3\mathrm{\&theta;})---\left(10a\right)$

$M_x\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2{C}_{\mathrm{Mw}}+\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2(A_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+A_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+A_3\mathrm{\&theta;})---\left(10b\right)$

其中，

H₁=(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w

H₂=[(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w]r/B

H₃=K_Dwsinψ_w+K_Lwcosψ_w

A₁=－K_Mw,A₂=－K_Mwr/B,A₃=K_Mw    (11)

式中，H_i和A_i(i=1,2,3)是准气动导数；其他符号意义同前。

由

、

和θ引起的自激效应实际上是与主梁振动折减频率有关的，因此式(10)需通过引入修正系数重新表达如下：

$F_y\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(C_{\mathrm{Dw}}\sin {\mathrm{\&psi;}}_w+C_{\mathrm{Lw}}\cos {\mathrm{\&psi;}}_w)+\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(h_1{H}_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+h_2{H}_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+h_3{H}_3\mathrm{\&theta;})---\left(12a\right)$

$M_x\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2{C}_{\mathrm{Mw}}+\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2(a_1{A}_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+a_2{A}_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+a_3{A}_3\mathrm{\&theta;})---\left(12b\right)$

式中，h_i和a_i(i=1,2,3)是和折减频率（或折减风速）相关的修正系数；其他符号意义同前。

令u=w=0，则V_r≈U_r=U，那么式(12)中的自激力可以表达为

$F_y^{\mathrm{se}}\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{U}^{2}B(h_1{H}_1^{\&prime;}\frac{\stackrel{\&CenterDot;}{h}}{U}+h_2{H}_2^{\&prime;}\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U}+h_3{H}_3^{\&prime;}\mathrm{\&theta;})---\left(13a\right)$

$M_x^{\mathrm{se}}\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{U}^2{B}^2(-a_1{K}_{M0}\frac{\stackrel{\&CenterDot;}{h}}{U}-a_2{K}_{M0}\frac{r}{B}\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U}+a_3{K}_{M0}\mathrm{\&theta;})---\left(13b\right)$

其中，

$H_1^{\&prime;}=(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0---\left(14\right)$

$H_2^{\&prime;}=[(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0]\frac{r}{B}---\left(15\right)$

$H_3^{\&prime;}=K_{D0}\sin {\mathrm{\&alpha;}}_0+K_{L0}\cos {\mathrm{\&alpha;}}_0---\left(16\right)$

$C_{D0}=C_D\left({\mathrm{\&alpha;}}_0\right),C_{L0}=C_L\left({\mathrm{\&alpha;}}_0\right),K_{D0}=\frac{d{C}_D}{\mathrm{d\&alpha;}}{|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0},K_{L0}={\frac{d{C}_L}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0},K_{M0}={\frac{d{C}_M}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0}---\left(17\right)$

式中，

和

分别为竖向自激力和扭转自激力；

、

和

为简化表达需要，无实际物理意义；C_D0和C_L0分别为初始风攻角α₀对应的阻力系数和升力系数；K_D0、K_L0和K_M0分别为初始风攻角α₀对应的阻力系数导数、升力系数导数和升力矩系数导数；其他符号意义同前。

该气动力模型有两点不足：(1)该模型认为风荷载沿桥面宽度是完全相关的，即要求脉动风的相关距离远大于主梁宽度尺寸，这对于脉动风中的高频分量是很难做到的；(2)该模型所考虑的自激作用与结构的频率无关，只有在较低的结构折算频率情况下才有效，即要求风横向吹过主梁的时间远小于结构的振动周期。上述两个缺点均源于对静力三分力系数的直接采用，而未对其作出修正。关于上述第一点的影响，可以通过在脉动风场模拟中引入气动导纳函数得到一定的补偿；针对第二点不足，可以通过引入由颤振导数计算获得的修正系数对气动力中的自激力成分进行修正。

根据Scanlan的建议，桥梁结构主梁单位长度上的气动自激力可以用六个颤振导数表示为

$F_y^{\mathrm{se}}\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{U}^2\left(2B\right)(K{H}_1^{*}\frac{\stackrel{\&CenterDot;}{h}}{U}+K{H}_2^{*}\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U}+K^2{H}_3^{*}\mathrm{\&theta;})---\left(18a\right)$

$M_x^{\mathrm{se}}=\frac{1}{2}\mathrm{\&rho;}{U}^2\left(2B^2\right)(K{A}_1^{*}\frac{\stackrel{\&CenterDot;}{h}}{U}+K{A}_2^{*}\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U}+K^2{A}_3^{*}\mathrm{\&theta;})---\left(18b\right)$

式中，K=ωB/U，为折减频率，ω为振动圆频率；

和

是颤振导数，均为折减频率K的函数，与桥梁断面的几何构形和来流有关，可通过桥梁断面的节段模型风洞试验或CFD计算获取；其他符号意义同前。

将式(18)分别与式(13)作对应部分的比较，可得

$h_1\frac{2}{(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0}K{H}_1^{*}---\left(19\right)$

$h_2=\frac{2B}{[(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0]r}K{H}_2^{*}---\left(20\right)$

$h_3=\frac{2}{K_{D0}\sin {\mathrm{\&alpha;}}_0+K_{L0}\cos {\mathrm{\&alpha;}}_0}{K}^2{H}_3^{*}---\left(21\right)$

$a_1=\frac{-2}{K_{M0}}K{A}_1^{*},a_2=\frac{-2B}{K_{M0}r}K{A}_2^{*},a_3=\frac{2}{K_{M0}}{K}^2{A}_3^{*}---\left(22\right)$

本发明在计算修正系数时进行关于振动频率的迭代，即先用零风速时的频率计算修正系数h_i和a_i(i=1,2,3)，计算（较少的时间步）后得出振动频率再算，直到前后两次的频率接近为止。事实上，为了加快收敛速度，可以利用前一级风速时的振动频率作为迭代初值。将h_i和a_i(i=1,2,3)的值代入式(12)，即可获得修正后的气动力F_y(t)和M_x(t)。至于气动力F_z(t)，可根据式(1a)、(2a)和(2b)直接算得。

下面介绍大跨度桥梁风致灾变全过程模拟方法的实现步骤。在风荷载作用下，桥梁运动微分方程可以表示为：

$M\stackrel{\&CenterDot;\&CenterDot;}{X}+C\stackrel{\&CenterDot;}{X}+\mathrm{KX}=F(X,\stackrel{\&CenterDot;}{X})---\left(23\right)$

式中，M、C和K分别为质量矩阵、阻尼矩阵和刚度矩阵；、

和X分别为加速度向量、速度向量和位移向量；

为用统一气动力模型表示的风荷载列阵，它除了与风速有关外，还是结构位移与速度的非线性函数。由于考虑几何非线性的影响，且作用在大跨度桥梁上的气动力与每时刻结构的响应有关，因此，风荷载作用下桥梁运动微分方程式是一个非线性方程，可以通过直接积分法进行求解。

如图2所示，本发明的大跨度桥梁风致灾变全过程模拟方法，包括以下步骤：

第一步：建立桥梁有限元模型，确定桥梁结构在恒载作用下的初始状态；

第二步：进行无阻尼的结构模态分析，提取其一阶竖弯频率f_v0和一阶扭转频率f_t0；

第三步：给定初始平均风速U₁和风速增量ΔU；

第四步：在当前平均风速U_i，采用谐波合成法生成水平和垂直脉动风时程，并读入内存；假定系统的竖弯振动频率f_vi和扭转振动频率f_ti分别为上级风速下的竖弯振动频率f_vi-1和扭转振动频率f_ti-1，当i=1时即为结构的一阶竖弯频率f_v0和一阶扭转频率f_t0；

第五步：根据初始风攻角和各单元的扭转变形均值计算各单元的有效攻角；然后根据当前风速、试算振动频率和有效攻角二次插值计算各单元的颤振导数，即先在不同初始攻角下相应的折减风速处进行第一次插值计算，而后根据实际的有效攻角进行第二次插值计算；这样可以有效地考虑气动导数沿展向的不均匀分布；

第六步：计算统一气动力荷载，进行动力时程分析；

F_z(t)＝D(t)cos(α₀+ψ)-L(t)sin(α₀+ψ)    (1)

$F_y\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(C_{\mathrm{Dw}}\sin {\mathrm{\&psi;}}_w+C_{\mathrm{Lw}}\cos {\mathrm{\&psi;}}_w)+\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(h_1{H}_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+h_2{H}_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+h_3{H}_3\mathrm{\&theta;})---\left(2\right)$

$M_x\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2{C}_{\mathrm{Mw}}+\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2(a_1{A}_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+a_2{A}_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+a_3{A}_3\mathrm{\&theta;})---\left(3\right)$

其中，

$D\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B{C}_D\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(4\right)$

$L\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B{C}_L\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(5\right)$

α_e(t)＝α₀+ψ(t)+θ(t)    (6)

$V_r\left(t\right)={\{U_r^2\left(t\right)+{[w\left(t\right)-\stackrel{\&CenterDot;}{h}-r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}]}^2\}}^{1/2}---\left(7\right)$

$\mathrm{\&psi;}\left(t\right)=[w\left(t\right)-\stackrel{\&CenterDot;}{h}-r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}]/U_r\left(t\right)---\left(8\right)$

$U_r\left(t\right)=U+u\left(t\right)-\stackrel{\&CenterDot;}{p}---\left(9\right)$

$h_1=\frac{2}{(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0}K{H}_1^{*}---\left(10\right)$

$h_2=\frac{2B}{[(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0]r}K{H}_2^{*}---\left(11\right)$

$h_3=\frac{2}{K_{D0}\sin {\mathrm{\&alpha;}}_0+K_{L0}\cos {\mathrm{\&alpha;}}_0}{K}^2{H}_3^{*}---\left(12\right)$

$a_1=\frac{-2}{K_{M0}}K{A}_1^{*},a_2=\frac{-2B}{K_{M0}r}K{A}_2^{*},a_3=\frac{2}{K_{M0}}{K}^2{A}_3^{*}---\left(13\right)$

$C_{D0}=C_D\left({\mathrm{\&alpha;}}_0\right),C_{L0}=C_L\left({\mathrm{\&alpha;}}_0\right),K_{D0}=\frac{d{C}_D}{\mathrm{d\&alpha;}}{|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0},K_{L0}={\frac{d{C}_L}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0},K_{M0}={\frac{d{C}_M}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0}---\left(14\right)$

ψ_w＝α₀+w/U_r    (15)

$C_{\mathrm{Dw}}=C_D\left({\mathrm{\&psi;}}_w\right),C_{\mathrm{Lw}}=C_L\left({\mathrm{\&psi;}}_w\right),C_{\mathrm{Mw}}=C_M\left({\mathrm{\&psi;}}_w\right),K_{\mathrm{Dw}}={\frac{d{C}_D}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w},K_{\mathrm{Lw}}={\frac{d{C}_L}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w},K_{\mathrm{Mw}}={\frac{d{C}_M}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w}---\left(16\right)$

H₁=(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w          (17)

H₂=[(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w]r/B     (18)

H₃=K_Dwsinψ_w+K_Lwcosψ_w                        (19)

A₁=－K_Mw,A₂=－K_Mwr/B,A₃=K_Mw                   (20)

在上述表达式中，F_z(t)、F_y(t)和M_x(t)分别为主梁单位长度上任一时刻在总体坐标系上所受的气动阻力、升力和升力矩；D(t)和L(t)分别为主梁单位长度上任一时刻在相对风轴坐标系上的气动阻力和升力；ρ为空气密度；B为主梁宽度；U为主梁高度处的平均风速，u(t)、w(t)分别为顺风向和竖向脉动风速，U_r(t)为来流风轴合成速度，V_r(t)为相对风速；α₀为初始风攻角，α_e(t)为瞬时风攻角，ψ(t)为来流风轴坐标系内的瞬时风攻角分量，ψ_w为总体坐标系内竖向脉动风速产生的瞬时风攻角分量；p、h和θ分别为主梁振动的水平、竖向和扭转位移，、

和分别为主梁竖向、水平和扭转振动速度；C_D、C_L和C_M是主梁风轴静力三分力系数（参考长度均为主梁宽度B），分别是阻力系数、升力系数和升力矩系数，它们是风攻角的函数；C_D0和C_L0分别为初始风攻角α₀对应的阻力系数和升力系数；K_D0、K_L0和K_M0分别为初始风攻角α₀对应的阻力系数导数、升力系数导数和升力矩系数导数；C_Dw、C_Lw和C_Mw分别ψ_w对应的阻力系数、升力系数和升力矩系数；K_Dw、K_Lw和K_Mw分别ψ_w对应的阻力系数导数、升力系数导数和升力矩系数导数；r为主梁特征长度，反映气动力作用点与截面形心之间距离的系数，取B/4；K=ωB/U，为折减频率，ω为振动圆频率；和

是颤振导数，均为折减频率K的函数，与桥梁断面的几何构形和来流有关；h_i和a_i(i=1,2,3)是和折减频率（或折减风速）相关的修正系数；H_i和A_i(i=1,2,3)是准气动导数。

然后对竖弯和扭转位移响应进行频谱分析，确定竖弯振动卓越频率f_v和扭转振动卓越频率f_t。将f_v和f_t与试算值f_vi和f_ti分别进行比较，如果|(f_v-f_vi)/f_vi|<0.05且|(f_t–f_ti)/f_ti|<0.05，说明频率迭代收敛。根据位移响应判断是否出现了颤振发散，随着风速的增长，如果振动形式逐渐从随机振动过渡到谐波发散振动，振幅逐渐增大,相应振动的阻尼将逐渐减小到0时就是颤振。如果出现了颤振发散，则结束计算；否则，增加风速，进入第四步。如果|(f_v-f_vi)/f_vi|>0.05或|(f_t–f_ti)/f_ti|>0.05，则计算各单元的扭转变形均值，并令f_vi=f_v，f_ti=f_t，重复第五步和第六步；需要说明的是，扭转位移均值对试算频率的选取是不敏感的，所以不需要寻求关于扭转位移均值的迭代收敛；

第七步：根据各级风速的位移响应结果评价抖振性能和判断静风失稳：在较低风速下通过均值、根方差和功率谱等统计指标来评价结构的抖振性能；以位移响应均值的突然增大来判断静风失稳。

根据各风速下的位移响应，可以判断桥梁结构在某级风速下的稳定状态：在较低风速下可通过均值、根方差和功率谱等统计指标来评价结构的抖振性能；随着风速的增长，如果振动形式逐渐从随机振动过渡到谐波发散振动，振幅逐渐增大,相应振动的阻尼将逐渐减小到0时就是颤振；以位移响应均值的突然增大来判断静风失稳。

另外，本发明除了解决上述问题外，还针对现有技术存在的其他问题进行了改进，包括：

首先，现有技术中自激力采用响应函数列式方法，表达为脉冲函数与桥面运动的卷积形式，能充分考虑气动导数与主梁振动折算频率有关的影响，但在自激力时域化过程中用到复杂的数学理论，且存在传递函数的取用以及其中涉及的非线性参数识别问题，不便于工程应用，而且待识别的参数较多，可能存在误差传递和累积。本发明通过引入修正系数对气动力中的自激力成分进行修正，并通过频率的迭代实现了自激力成分时域化，避免了参数识别问题。

其次，现有技术中，抖振力的列式略去了脉动风速平方项的影响，这在紊流强度较大时候会引起较大的误差。本发明计算气动力时采用的是相对风速的平方，相对风速包含了水平脉动风速和竖向脉动风速，因此计算所得的气动力包含了脉动风速高阶项。

总之，本发明的整个计算过程中实现了批处理式调用，采用了多种参数传递技术，使得计算更加自动化，使用起来更加方便。

Claims (2)

1.大跨度桥梁风致灾变全过程模拟方法，其特征在于它包括以下步骤：

第一步：建立桥梁有限元模型，确定桥梁结构在恒载作用下的初始状态；

第二步：进行无阻尼的结构模态分析，提取其一阶竖弯频率f_v0和一阶扭转频率f_t0；

第三步：给定初始平均风速U₁和风速增量ΔU；

第四步：设当前平均风速U_i，采用谐波合成法生成水平和垂直脉动风时程；假定系统的竖弯振动频率f_vi和扭转振动频率f_ti分别为上级风速下的竖弯振动频率f_vi-1和扭转振动频率f_ti-1，当i=1时即为桥梁结构的一阶竖弯频率f_v0和一阶扭转频率f_t0；

第五步：根据初始风攻角和各单元的扭转变形均值计算各单元的有效攻角；然后根据当前风速、试算振动频率和有效攻角二次插值计算各单元的颤振导数，先在不同初始攻角下相应的折减风速处进行第一次插值计算，而后根据实际的有效攻角进行第二次插值计算；

第六步：计算统一气动力荷载，进行动力时程分析，具体如下列公式表示：

F_z(t)＝D(t)cos(α₀+ψ)-L(t)sin(α₀+ψ)    (1)

$F_y\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(C_{\mathrm{Dw}}\sin {\mathrm{\&psi;}}_w+C_{\mathrm{Lw}}\cos {\mathrm{\&psi;}}_w)+\frac{1}{2}\mathrm{\&rho;}{V}_r^{2}B(h_1{H}_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+h_2{H}_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+h_3{H}_3\mathrm{\&theta;})---\left(2\right)$

$M_x\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2{C}_{\mathrm{Mw}}+\frac{1}{2}\mathrm{\&rho;}{V}_r^2{B}^2(a_1{A}_1\frac{\stackrel{\&CenterDot;}{h}}{U_r}+a_2{A}_2\frac{B\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}{U_r}+a_3{A}_3\mathrm{\&theta;})---\left(3\right)$ 其中，

$D\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B{C}_D\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(4\right)$

$L\left(t\right)=\frac{1}{2}\mathrm{\&rho;}{V}_r^2\left(t\right)B{C}_L\left[{\mathrm{\&alpha;}}_e\left(t\right)\right]---\left(5\right)$

α_e(t)＝α₀+ψ(t)+θ(t)    (6)

$V_r\left(t\right)={\{U_r^2\left(t\right)+{[w\left(t\right)-\stackrel{\&CenterDot;}{h}-r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}]}^2\}}^{1/2}---\left(7\right)$

$\mathrm{\&psi;}\left(t\right)=[w\left(t\right)-\stackrel{\&CenterDot;}{h}-r\stackrel{\&CenterDot;}{\mathrm{\&theta;}}]/U_r\left(t\right)---\left(8\right)$

$U_r\left(t\right)=U+u\left(t\right)-\stackrel{\&CenterDot;}{p}---\left(9\right)$

$h_1=\frac{2}{(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0}K{H}_1^{*}---\left(10\right)$

$h_2=\frac{2B}{[(C_{L0}-K_{D0})\sin {\mathrm{\&alpha;}}_0-(C_{D0}+K_{L0})\cos {\mathrm{\&alpha;}}_0]r}K{H}_2^{*}---\left(11\right)$

$h_3=\frac{2}{K_{D0}\sin {\mathrm{\&alpha;}}_0+K_{L0}\cos {\mathrm{\&alpha;}}_0}{K}^2{H}_3^{*}---\left(12\right)$

$a_1=\frac{-2}{K_{M0}}K{A}_1^{*},a_2=\frac{-2B}{K_{M0}r}K{A}_2^{*},a_3=\frac{2}{K_{M0}}{K}^2{A}_3^{*}---\left(13\right)$

$C_{D0}=C_D\left({\mathrm{\&alpha;}}_0\right),C_{L0}=C_L\left({\mathrm{\&alpha;}}_0\right),K_{D0}=\frac{d{C}_D}{\mathrm{d\&alpha;}}{|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0},K_{L0}={\frac{d{C}_L}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0},K_{M0}={\frac{d{C}_M}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&alpha;}}_0}---\left(14\right)$

ψ_w＝α₀+w/U_r    (15)

$C_{\mathrm{Dw}}=C_D\left({\mathrm{\&psi;}}_w\right),C_{\mathrm{Lw}}=C_L\left({\mathrm{\&psi;}}_w\right),C_{\mathrm{Mw}}=C_M\left({\mathrm{\&psi;}}_w\right),K_{\mathrm{Dw}}={\frac{d{C}_D}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w},K_{\mathrm{Lw}}={\frac{d{C}_L}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w},K_{\mathrm{Mw}}={\frac{d{C}_M}{\mathrm{d\&alpha;}}|}_{\mathrm{\&alpha;}={\mathrm{\&psi;}}_w}---\left(16\right)$

H₁=(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w         (17)

H₂=[(C_Lw－K_Dw)sinψ_w－(C_Dw+K_Lw)cosψ_w]r/B    (18)

H₃=K_Dwsinψ_w+K_Lwcosψ_w                       (19)

A₁=－K_Mw,A₂=－K_Mwr/B,A₃=K_Mw                  (20)

在上述表达式中，F_z(t)、F_y(t)和M_x(t)分别为主梁单位长度上任一时刻在总体坐标系上所受的气动阻力、升力和升力矩；D(t)和L(t)分别为主梁单位长度上任一时刻在相对风轴坐标系上的气动阻力和升力；ρ为空气密度；B为主梁宽度；U为主梁高度处的平均风速，u(t)、w(t)分别为顺风向和竖向脉动风速，U_r(t)为来流风轴合成速度，V_r(t)为相对风速；α₀为初始风攻角，α_e(t)为瞬时风攻角，ψ(t)为来流风轴坐标系内的瞬时风攻角分量，ψ_w为总体坐标系内竖向脉动风速产生的瞬时风攻角分量；p、h和θ分别为主梁振动的水平、竖向和扭转位移，

、

和

分别为主梁竖向、水平和扭转振动速度；C_D、C_L和C_M是主梁风轴静力三分力系数（参考长度均为主梁宽度B），分别是阻力系数、升力系数和升力矩系数，它们是风攻角的函数；C_D0和C_L0分别为初始风攻角α₀对应的阻力系数和升力系数；K_D0、K_L0和K_M0分别为初始风攻角α₀对应的阻力系数导数、升力系数导数和升力矩系数导数；C_Dw、C_Lw和C_Mw分别ψ_w对应的阻力系数、升力系数和升力矩系数；K_Dw、K_Lw和K_Mw分别ψ_w对应的阻力系数导数、升力系数导数和升力矩系数导数；r为主梁特征长度，反映气动力作用点与截面形心之间距离的系数；K=ωB/U，为折减频率，ω为振动圆频率；

和

是颤振导数，均为折减频率K的函数，与桥梁断面的几何构形和来流有关；h_i和a_i(i=1,2,3)是和折减频率（或折减风速）相关的修正系数；H_i和A_i(i=1,2,3)是准气动导数；

第七步：根据各级风速的位移响应结果评价抖振性能和判断静风失稳：在较低风速下，通过均值、根方差和功率谱来评价结构的抖振性能；以位移响应均值出现峰值来判断静风失稳。

2.根据权利要求1所述的大跨度桥梁风致灾变全过程模拟方法，其特征在于第六步进一步包括：对竖弯和扭转位移响应进行频谱分析，确定竖弯振动卓越频率f_v和扭转振动卓越频率f_t，并将f_v和f_t与试算值f_vi和f_ti分别进行比较，

如果|(f_v-f_vi)/f_vi|<0.05且|(f_t–f_ti)/f_ti|<0.05，说明频率迭代收敛；根据位移响应判断是否出现了颤振发散，随着风速的增长，如果振动形式逐渐从随机振动过渡到谐波发散振动，振幅逐渐增大,相应振动的阻尼将逐渐减小到0时就是颤振；如果出现了颤振发散，则结束计算；否则，增加风速，进入第四步；

如果|(f_v-f_vi)/f_vi|>0.05或|(f_t–f_ti)/f_ti|>0.05，则计算各单元的扭转变形均值，并令f_vi=f_v，f_ti=f_t，重复第五步和第六步。

Images