基于ANSYS的桥梁颤振时域分析及程序实现

刘修平， 刘焕举， 赵 越， 李光玲

(1.天水师范学院 土木工程学院，甘肃 天水 741000;2.河北工程大学 土木工程学院，河北 邯郸 056038;3.长安大学 公路学院，西安 710064)

颤振是一种发散性气动力自激振动，是桥梁风致振动当中最危险的振动形式。自Scanlan引入颤振导数描述钝体桥梁自激力以来，分离流颤振机理在桥梁颤振分析应用中日趋完善[1-3]。然而随着大跨缆索承重桥梁在跨越江海、深切峡谷的国家公路和铁路交通控制性工程广泛使用，结构非线性、气动非线性等非线性因素对桥梁稳定性的影响也日益突出，风致颤振稳定性成了大跨桥梁设计施工过程中需考虑的关键性问题。由于Scanlan线性自激力模型不能很好的描述考虑各类非线性因素在颤振中的影响，在桥梁抗风研究日益精细的情况下，自激力时域模型研究成为了风致颤振分析的热点[4-6]。

自激力的时域表达式有两种方法：其一是Scanlan提出的Wagner在航空中提出的经典阶跃函数表达法，其二是Buther & Lin提出的单位脉冲响应函数融入Roger有理函数的表达方法。张志田等[7]分别研究了两种自激力时域表达式的瞬态特性及极限特性，提出了通过限定有理函数表达式若干参数克服瞬态特性模拟失真的方法；Chen等[8]基于有理函数表达的自激气模型，提出了桥梁大跨度桥梁结构颤振分析的时域分析流程及特征值求解方法；郭增伟等[9-10]通过颤振导数的有理函数近似式得到自激力脉冲响应函数在桥梁颤振分析中的表现方式和求解方案，提出了颤振分析中初始响应失真及积分时间步长设定的解决方法；伍波等[11-12]发现相比经典阶跃函数颤振分析方法，在考虑静风效应的大跨度桥梁结构颤振特性分析中，有理函数自激力表达式模型可以有效减少计算风速下的迭代量，显著提高计算效率。由于桥梁颤振非线性问题分析和方程求解的复杂性，一直以来，桥梁颤振时域分析多是基于C/C++、Fortran开发的专用程序进行[13-17]。近些年，随着大型通用有限元分析程序在桥梁工程中具有广泛的应用，基于通用有限元软件ANSYS的颤振频率分析方法日趋成熟，直接采用ANSYS进行桥梁颤振时域分析的研究较少。华旭刚等[18]基于ANSYS中MATRIX27矩阵单元建立了桥梁结构全模态颤振分析模型，提出以调整风速和振动频率求解动力系统各阶复模态特性的颤振频域分析方法，曾宪武等[19]基于MATRIX27矩阵单元对自激力以单元气动阻尼矩阵和单元气动刚度矩阵提出了基于ANSYS的大跨桥梁抖振分析方法；为避免在使用过程中需假定振动频率的步骤，谢炼等[20]利用了ANSYS内置的重启动算法，迭代计算当前时步真实振动状态，提出了桥梁颤振时域计算的重启动法。

为避免MATRIX27矩阵单元输入振动频率的误差干扰，同时简化计算时步内的自激力反复迭代确定的重复过程。本文基于传统有理函数自激力表达式推导处理，使得自激力表达式中的“记忆效应”项仅与上一时刻的桥面运动状态相关联，继而重新形成有理函数自激力有表达式的质量项、刚度项、阻尼项和时间历程项，建立基于MATRIX27矩阵单元的自激力时域模拟方法，并通过经典平板颤振分析成功验证了所提方法的准确性。

1 桥梁断面自激力

1.1 自激力有理函数表达

(1)

p=K(ξ+i)，ξ为桥梁结构阻尼比一般较小，可以忽略；其他参数含义如表1所示。

表1 8个颤振导数的参与振动方向Tab.1 Directions of eight aerodynamic derivatives

Scanlan自激力列式和振动圆频率关联，不具备直接进行时域分析计算的条件。Lin提出可以完全在时域中描述风致振动中产生的耦合气动自激力的关系式，将之表达为

(4)

根据式(3)，有理函数推广用于近似桥梁断面自激力的气动矩阵Q，形式为

(5)

式中：C1，…Cl+3均与频率无关的气动力系数；C1为自激力的气动刚度项；C2为气动阻尼项；C3为气动质量项，其值一般较小，可以忽略；第4项以后都为自激力的记忆效应，与结构前一时刻和当前时刻的运动有关；λl为大于0的非定常项衰减常数。

1.2 气动力系数识别

颤振导数的拟合对象为高度非线性的方程组，为避免大量的复数运算过程，令上述方程的实部和虚部相等，以升力的计算为例，则可以得到如下的关系式

(6)

该组方程式对应7个未知量，拟合过程需要参数共享，同时也要避免因初始值的影响陷入局部最优的问题。对于二维弯扭耦合颤振导数拟合共需识别28个气动力系数，因此，优化迭代算法的选择至关重要。

以理想平板颤振导数拟合为例，在最大迭代次数2 000、收敛误差为1×10-10的设定下，分别采用常用的麦夸特法、拟牛顿法、遗传算法、模拟退火算法、拟牛顿+通用全局优化算法求解上述方程组(2)式，求解结果和过程参数对比如表2所示。

表2 不同算法求解结果比较Tab.2 Comparison of the results of different algorithms

由表2可知，在颤振导数拟合的实际计算中，求解此类方程组，麦夸特法、拟牛顿法和拟牛顿+通用全局优化算法都能快速的得到收敛解，而遗传算法和模拟退火算法的计算效率较低，达到限定迭代步数后仍未收敛，同时对非定常项衰减常数λ计算值出入较大，甚至出现负值，显然，考虑卷积积分Φsel的存在，λ应是大于0的实数。因此，在气动力系数识别过程中应优先采用拟牛顿+通用全局优化算法。

2 基于MATRIX27的自激力时域模拟

2.1 MATRIX27矩阵单元的组集

(7)

(8)

(9)

MATRIX27矩阵单元无几何外形特性，其运动学响应用刚度、阻尼或者质量系数来指定，该矩阵单元连接2个节点，每个节点有6个自由度。在通过式(8)、(9)集成式(7)中的各矩阵后，将矩阵内的元素以实常数的方式赋予MATRIX27单元，并将其添加至结构单元的节点处，如图3所示，时间历程项作为外荷载形式输入，从而形成时域自激力桥梁颤振分析有限元模型。

2.2 非线性项处理方法

时域自激力计算的难点在于非线性项的荷载受当前的时步ti和上一时步ti-1节点运动状态共同影响，而商业有限元软件ANSYS的瞬态分析是基于隐式算法，无法推算出当前时步的荷载。虽然可以利用重启动实现了当前时步荷载的计算，但重启动每个计算时步需迭代确定，在通用有限元程序中计算步数相对较多，耗时较长，对于大跨桥梁风致抖振分析而言效率并不高。因此，真正意义地实现基于MATRIX27矩阵单元的ANSYS时域颤振计算的前提就是将有理函数自激力表达式重新改写，具体如下：

(1) 将非线性项重新改写，以单位长度上扭转运动对扭转的影响为例。

(10)

(11)

引入，

那么，

(12)

利用Φsel(ti)和Φsel(ti-1)的关系可以得到

(13)

(2) 形成自激力有理函数表达的质量项、刚度项、阻尼项和时间历程项。

在非线性项改写的基础上，那么，式(10)和式(11)可以表达为

(14)

(15)

令

(16)

2.3 时程分析实现过程

桥梁颤振时域分析过程中，首先要颤振分析有限元模型施加一个初始激励，由于自激振动系统的特征量，如频率、衰减率等与初始激励无关，在颤振分析时可对结构施加初始激励，激励形式可为集中力或者位移等，在ANSYS中的具体实现过程为：

(1) 假定计算风速。

(2) 施加任意初始荷载激励。

(3) 计算第ti时刻的位移、速度和加速度向量。

(4) 计算ti+1时刻的荷载向量。

(5) 计算ti+1时刻的位移、速度和加速度向量。

(6) 依次重复步骤(4)、(5)，直到设置时间步完成。

(7) 计算此风速下结构振动时程响应的对数衰减率r，若r>0增加风速返回步骤(2)，若r<0减小风速返回步骤(2)，若r=0，则为颤振临界风速，程序结束。

在步骤(7)中的颤振临界风速确定时，可采用折半法查找，首先确定颤振临界风速大致区间[Ul,Uh]，计算风速取Um=(Ul+Uh)/2，再进行迭代计算。同时，由于ANSYS中的APDL语言编程的限制，节点振动时程响应对数衰减率的计算可由MATLAB中的峰值读数法函数直接进行，通过增加最小间隔条件，可以过滤掉单周期内的干扰极值，更准确的判断振动是否衰减，在程序中迭代框架具体见图1。

图1 颤振时域分析实现框架Fig.1 The implementation framework of time-domain flutter analysis

3 算例验证

3.1 理想平板断面的简支梁

为验证所提时域自激力计算的准确性，以理想平板的颤振分析为例，并将计算结果和理论解对比分析。该简支梁长度为300 m，桥面宽度40 m，材料弹性模量取2.1×1012Pa，横向抗弯惯性取8.57 m4，竖向抗弯惯矩取1.0×104m4，扭转惯矩取0.386 m4，质量取2×104kg/m，质量惯矩取4.5×106kg·m，空气密度取1.225 kg/m3，简支梁前10阶的自振特性如表3所示，理想平板颤振导数如图2所示。

(a) H*

(b) A*图2 理想平板的颤振导数Fig.2 Aerodynamic derivatives of ideal plate

表3 简支梁前10阶自振特性Tab.3 The first ten order frequencies of simply supported beam bridge

在MATRIX27矩阵单元的颤振分析模型建立过程中，首先将简支梁划分为30个单元，将质量等效到节点上，梁端节点约束扭转自由度，引入29个气动刚度矩阵、29个气动阻尼矩阵和29个气动质量矩阵，颤振分析有限元模型如图3所示，模型共包括单元210个，节点122个，基于拟牛顿+通用全局优化算法对理想平板颤振导数拟合后得到的气动力系数如表4所示。

图3 简支梁颤振分析有限元模型Fig.3 Flutter analysis FE model of simply supported beam bridge

表4 理想平板气动力系数Tab.4 Aerodynamic coefficients of ideal plate

从表4中可见：φlα和φmh的气动质量项足够小，可以忽略不计，而φlh和φmα的拟合值并非如此，因此，在进行气动导数拟合时，不考虑气动质量项可能使得拟合结果出现偏差，可能影响颤振临界风速的计算结果。

简支梁在颤振发生及前后时的跨中竖向位移和扭转位移的响应时程如图4所示，由图可见：随着风速的增加，简支梁跨中位置振动时程曲线由衰减到逐步发散，在颤振临界风速处，简支梁维持等幅振动。

(a) 竖向位移时程

(b) 扭转位移时程图4 颤振前后简支梁跨中位移响应Fig.4 Displacement responses of mid-span of simply supported beam bridge before and after flutter

表5为本文方法与其他方法计算结果的比较，由表可知：本文所提方法的计算结果与理论解吻合良好，误差不足1%。

若不重新改写非线性Φsel，直接由Newmark-β进行时程求解上述自激力模型时，则会出现文献[9]中的结果，结构的颤振临界风速值受时间步影响较大，这是由于计算时步的真实荷载受当前和上一时刻桥梁运动状态共同的影响，直接递推得到的自激力在相位差和大小上与实际值存在一定的差距。为研究时间步长对本文所提方法的影响程度，分别计算不同的时间步长下的颤振临界风速值，计算结果如表6所示。由表6可知：对非线性Φsel处理后，时间步长的选取对颤振临界风速的计算结果依赖性不明显。

3.2 某大跨斜拉桥颤振临界风速计算

以某大跨度斜拉桥颤振分析为例，设计采用三塔斜拉桥方案，桥梁跨径具体布置为110 m+129 m+258 m+258 m+129 m+110 m，全长994 m，桥梁宽度为38.8 m，箱梁高为4.5 m，主梁标准断面如图5所示。

图5 桥梁标准断面(mm)Fig.5 Standard section of a bridge (mm)

图6为该主梁断面节段模型0°风攻角下的(几何缩尺比为1∶50)颤振导数试验及拟合结果(拟合值为无符号的曲线)，表7为气动力系数计算结果。

表7 桥梁断面平板气动力系数Tab.7 Aerodynamic coefficients of bridge section

图6 0°风攻角下的颤振导数和拟合结果Fig.6 Flutter derivatives and fitting results under 0° wind attack angle

(a) 竖向位移时程

(b) 扭转位移时程

0°风攻角下，试验颤振临界风速>17 m/s，换算实桥颤振临界风速>110.5 m/s，通过本文方法计算得到颤振临界风速Ucr=198.6 m/s，颤振临界频率fcr=0.930 8 Hz，远高于成桥状态颤振检验风速79.88 m/s，说明该桥的颤振稳定性良好。图7为颤振临界风速条件下桥梁主梁断面竖向、扭转位移响应时程曲线，可以直观地判断结构正处于临界运动状态。

4 结 论

(1) 基于自激力有理函数表达式，推导了时域自激力递推计算过程，提出了一种基于ANSYS的桥梁颤振时域分析模型和求解方法，并采用MATLAB和ANSYS混合编程技术实现了颤振临界风速的程序求解。

(2) 在相同约束条件的限定下，相比其他几种优化算法，拟牛顿+通用全局优化算法在气动力系数识别过程中迭代步数更少，收敛更快，拟合效果更好。

(3) 以经典平板颤振为例，采用本文所提方法计算得到的颤振临界风速和颤振频率与理论吻合度很高，并且本文所提方法对计算时间步长的选取不具备依赖性，某大跨桥梁颤振临界风速计算说明所提方法在工程应用中具备实用性和可操作性。

(4) 伴随着桥梁抗风精细化研究的进程，结合通用有限元程序在处理非线性问题上优势，本文所提方法为复杂的非线性颤振问题研究分析提供了参考。