桥梁断面颤振导数的分状态多频强迫振动识别

王林凯, 刘志文， 陈政清

(1. 安徽省交通规划设计研究总院股份有限公司，合肥 230000；2. 湖南大学 风工程与桥梁工程湖南省重点实验室，长沙 410082)

桥梁断面颤振导数是研究桥梁气动性能的重要参数，可通过风洞试验或计算流体力学数值模拟获得。风洞试验识别桥梁颤振导数可以分为自由振动法与强迫振动法两类。自由振动法试验测试相对简单，Scanlan等[1]最早提出分阶段的颤振导数识别方法，即先通过竖弯和扭转单自由度振动试验识别出直接颤振导数，再用耦合振动试验识别出耦合项颤振导数。Poulsen等[2]通过风洞试验用自由振动方法识别了大带东桥梁断面的颤振导数。之后，国内外部分学者对自由振动识别方法进行了进一步研究与与应用[3-8]。强迫振动法则需要能够驱使桥梁断面以某一频率和振幅做简谐运动的专门装置。早在1952年，Halfman[9]采用强迫振动法直接测量了机翼非定常气动力，采用不同振动形式识别了不同攻角的颤振导数。Ukeguchi[10]首次将强迫振动测力法用于桥梁断面颤振导数的识别。Li[11]采用强迫振动法在水洞中识别了桥梁断面的颤振导数。陈政清等[12]率先在国内实现了用强迫振动法识别桥梁断面颤振导数，并且研究了颤振导数识别的时域和频域方法。郭震山[13]、牛华伟[14]在此基础上进行了三自由度桥梁断面的颤振导数识别与研究。相比自由振动识别方法，强迫振动法所需装置复杂，但是具有程序稳健性高，数据重复性好，可控制折算风速范围广等优点而在近年来的工程以及研究应用中被广泛地采用。

随着计算流体力学的发展，桥梁断面气动导数的强迫振动数值模拟识别方法逐渐受到关注。Larsen[15]借助离散涡方法模拟了大带东桥主梁截面的绕流场，并用分状态强迫振动法识别了该桥主梁断面的颤振导数。后来该方法得到了广泛的研究与应用[16-19]。崔益华等[20]使用了耦合强迫振动数值模拟识别法，识别了桥梁断面的颤振导数，并计算相应得桥梁颤振临界风速，计算结果与试验吻合较好，验证了该方法的可行性。

综上所述，目前关于桥梁断面气动导数试验识别采用耦合强迫振动法或自由振动识别法，数值模拟识别则以分状态单频强迫振动识别法为主。随着计算流体动力学的发展，从提高桥梁主梁断面颤振导数识别精度和效率的角度考虑，进行桥梁断面颤振导数快速识别方法的研究仍具有十分重要的价值和意义。本文在现有研究成果的基础上，分别进行了薄平板断面和流线型主梁断面颤振导数弯扭耦合两自由度强迫振动识别和分状态多频强迫振动识别数值模拟研究。

1 数值方法

1.1 流体力学控制方程

直角坐标系下，黏性不可压缩流体基于雷诺平均(RANS)的连续性方程和N-S方程分别为：

(1)

(2)

式中：i,j=1，ρ为空气密度，ρ=1.225 kg/m3，μ为动力黏性系数，μ=1.789 4×10-5kg/m·s。

1.2 湍流模型

(3)

涡黏模型包括零方程模型，一方程模型以及两方程模型。常用的两方程RANS涡黏模型有标准k-ε湍流，标准k-ω模型以及SSTk-ω模型等。

在标准k-ε模型和k-ω模型的基础上，Menter[23]提出了SST(Shear Stress Transport)k-ω模型，湍流方程表示为：

(4)

(5)

式中：μeff, k与μeff, ω表示湍流项的等效黏度;Sk与Sω表示湍动能和比耗散率的源项。

该模型将k-ε模型中的ε方程写成ω的形式，并且将其与标准k-ω模型进行线性组合，组合系数则表示为关于与壁面之间距离的函数。SSTk-ω模型湍流模型通过线性组合函数随与壁面之间距离的变化现实了从近壁面的k-ω模型向远离壁面的k-ε模型转变，从而结合了两种模型各自的特点，在黏性子层及远离壁面湍流充分发展区都具有很好的计算性能。

1.3 动网格技术

对于通量φ，在任一控制体V内，其边界运动时，N-S方程的积分形式表达式可以表示为：

(6)

式中：us为网格运动的速度。

为实现结构运动，采用“刚性网格区域+动网格区域”的方案，以薄平板为例的计算区域分区示意图，如图2所示。

图1 薄平板网格区域划分示意图(mm)Fig.1 Schematic diagram of plate grid division(mm)

薄平板附近包围一圈刚性运动区域网格，该部分网格随薄平板一致做刚性运动。在刚性运动区域网格外设置一层三角形动网格区域，该区域网格会在运动过程中重新划分网格。最外层为静止区域网格，该区域内网格在计算过程中不会发生任何变化。计算域左侧边界为速度入口，计算域右侧边界为压力出口，参考压力设为零；计算域上、下侧边界为对称边界条件，模型表面为无滑移边界。

动网格更新方法选择“弹性光顺+局部网格重构”的方式，既保证网格更新效率又可保证网格更新质量。动网格驱动宏选择DEFINE_CG_MOTION，以指定刚性网格运动区域的相应运动速度，进而实现网格位移的更新。

2 主梁断面气动导数识别方法

2.1 气动自激力描述

将Scanlan线性气动自激力表达式表述为基于耦合强迫振动的形式：

(7)

(8)

(9)

式中：fh与fα分别表示强迫振动的竖向振动频率和扭转振动频率。

2.2 弯扭耦合强迫振动气动导数识别方法

文献[20]中的弯扭耦合强迫振动识别方法为进行两次弯扭同频率的强迫振动识别，第一次弯扭位移同相位，第二次强迫振动弯扭位移设置为180°的相位差。为提高识别效率，本文选择文献[14]中的将弯扭频率设为不一致的识别方法。假设结构做耦合运动形式位移表示为：

(10)

通过CFD计算获取的气动自激力L与M的时程曲线，按照式(4)与式(5)的形式进行最小二乘拟合获取相应的颤振导数值。通过不同的弯扭频率进行错位组合识别出感兴趣折算风速下的颤振导数。

以CFD耦合强迫振动获取的升力系数和升力矩系数为数据来源，其三分力系数定义如下：

(11)

采用Matlab语言时域识别步骤可以表示为：

步骤1 赋予初值a=[0.1 0.1 0.1 0.1]；

步骤3 定义相关系数矩阵量为：

步骤4 定义目标拟合函数：

Fun=nlinfit(a,X)(a1N1X1+a2N2X2+a3N3X3+a4N4X4)

步骤5 进行非线性拟合：

H=nlinfit(X,CL,Fun,a)

A=nlinfit(X,CM,Fun,a)

2.3 分状态多频强迫振动气动导数识别方法

在Scanlan线性自激力可叠加性的理论框架下，推出一种新的识别方式。假设模型分别进行单自由度多频率强迫振动，即

(12)

每一振动项对颤振导数的影响系数取各自的振动频率。当结构做竖向分状态多频强迫振动时，对应的升力和升力矩系数分别为：

(13)

(14)

当结构做扭转分状态多频强迫振动时，对应的升力和升力矩系数分别为：

(15)

(16)

按式(12)分别驱动断面做竖向单自由度多频振动与扭转单自由度多频振动，通过两次计算即可获得主梁断面的气动力时程曲线，进而可识别出不同折算风速对应的桥梁断面气动导数。

3 典型断面气动导数识别

3.1 断面参数

分别选取薄平板[21]和大带东桥主梁断面进行断面气动导数识别研究。薄平板断面高D=20 mm，宽度B=450 mm，其外形如图2(a)所示。大带东桥主梁断面高D=4.4 m，宽B=31 m，高宽比为7.05，断面上有3%的横坡。该断面的剪切中心(S.C)位于距离桥面0.465倍截面高度处，其外形如图2(b)所示。

(a) 薄平板尺寸(mm)

(b) 大带东桥梁断面尺寸(m)图2 模型断面参数Fig.2 The Section parameters of models

3.2 网格划分与求解设置

利用Gambit绘制相应的网格如图3～4所示。薄平板边界处使用边界层4∶2网格能够有效的减少网格数量，并且与薄平板两端以结构化网格形式拓扑出去使得网格边长增大导致相邻网格长度比增大的现象相平衡，三角形网格数为23 840个，四边形网格数为13 208个。大带东桥主梁断面采用外包椭圆加圆的方式向外拓扑，四边形网格数量为39 941个，三角形网格数量为65 774，壁面处首层网格高度为1×10-5B。大带东桥梁断面网格文件导入Fluent中时需要进行80∶1缩尺，断面宽度即为0.387 5 m。

图3 薄平板网格示意图Fig.3 Schematic diagram of flat plate grid

图4 大带东桥梁断面网格图Fig.4 Grid of the main girder section of the Great Belt East Bridge

采用商业流体力学计算软件Ansys Fluent进行求解计算。计算选择速度-压力解耦的SIMPLEC算法，湍流模型选择SSTk-ω模型，压力方程采用二阶格式离散，动量方程、湍动能方程及比耗散率方程均采用二阶迎风格式离散，计算迭代残差为10-6，给予初始均匀流场湍流特征，湍流强度为0.5%，湍流黏性比为2。薄平板流场相对简单，湍流模型对计算结果的影响较小[24]。大带东桥主梁断面流场较为复杂，在强迫网格做简谐运动前，先对大带东桥梁断面做静止绕流计算以验证网格与求解设置的可靠性。大带东桥主梁断面静止绕流计算中，来流速度取12 m/s，计算时间步长取0.000 1 s，计算至流场充分稳定。St数定义如下：

(17)

式中：fs为涡脱频率；D为断面迎风面高度。计算结果如表1和图5所示。由表1可知，本文的大带东桥主梁断面网格以及求解设置满足精度要求，壁面处Y+分布符合湍流模型的要求。

图5 大带东桥主梁断面静止绕流计算结果Fig.5 The calculation results of flow of Great Belt East Bridge

3.3 弯扭耦合两自由度强迫振动气动导数识别

不同于文献[14]中试验过程中通过改变风速来改变折算风速，本文通过改变运动频率来改变折算风速，计算工况如表2和表3所示，各个工况中的频率的设置尽量使得折算风速为一整数值，为此薄平板算例中的来流速度值取10倍的薄平板宽度，即4.5 m/s，大带东桥梁断面的来流速度值取24倍的梁宽，即9.3 m/s。计算时，模型竖向振动的振幅为0.02B，扭转振幅为2°，根据采样定理，计算时间步长建议不大于模型驱动周期的0.02倍，本次统一取0.000 5 s，计算时间建议不小于各运动成分的10个周期时长。

薄平板颤振导数的识别结果如图6所示，大带东桥主梁断面颤振导数识别结果如图7所示。

表2 薄平板颤振导数计算工况Tab.2 Calculation conditions of flutter derivatives of the thin plate

表3 大带东桥梁断面颤振导数计算工况Tab.3 Calculation conditions of the flutter derivatives of the section of the Great Belt East Bridge

图6 薄平板颤振导数识别结果Fig.6 Flutter derivative identification results

从图6中可以看出，薄平板各个颤振导数均与平板理论解吻合较好，验证了本文的耦合强迫振动识别方法的正确性与可靠性。从图7中可以看出，大带东桥主梁断面颤振导数结果与相关文献数值模拟的结果较为接近，并且与Poulsen等的试验结果吻合较好。需要说明的是，祝志文等采用的是分状态单自由度强迫振动法识别的颤振导数。

图7 大带东桥梁断面颤振导数计算结果Fig.7 Calculation results of the flutter derivatives of girder section of the Great Belt East Bridge

3.4 气动自激力可叠加性CFD验证

(a) 升力时程曲线

(b) 力矩时程曲线图8 气动自激力对比图Fig.8 Comparison of self-excited aerodynamic forces

从图8中可以看出，由CFD获取的升力、力矩时程曲线与用耦合振动法识别的颤振导数拟合表达式计算的时程曲线几乎完全吻合，从侧面为耦合强迫振动法提供了理论支持，验证了该桥梁断面的气动自激力的可叠加性。同时也证明了CFD气动力计算的准确性并不因为模型结构的运动复杂性提高而改，这与文献[24]的结论是一致的。

3.5 分状态多频强迫振动气动导数识别

在线性气动自激力可叠加性的基础上，按表4和表5中工况进行分状态多频强迫振动识别薄平板和大带东桥主梁断面的颤振导数。

表4 薄平板分状态多频强迫振动颤振导数识别工况Tab.4 Calculation conditions of the identification of flutter derivatives of the thin plate by step-by-step forced vibration with multi frequency

表5 大带东桥主梁断面分状态多频强迫振动颤振导数识别工况Tab.5 Calculation conditions of the identification of flutter derivatives of the girder section of Great Belt East Bridge by step-by-step forced vibration with multi frequency

驱动断面运动的振幅仍然按耦合识别法设置，即竖弯0.02B，扭转位2°，计算获取的气动力系数时程曲线以及气动力时程曲线FFT幅值谱如图9～16所示。

(a) 升力时程曲线

(b) 力矩时程曲线图9 薄平板竖向多频强迫振动自激力时程曲线Fig.9 The history curves of self-excited forces of thin flat plate by forced vertical vibration with multi frequency

(a)升力时程曲线FFT幅值谱

(b)力矩时程曲线FFT幅值谱图10 薄平板竖向多频强迫振动自激力FFT幅值谱Fig.10 The FFT amplitude spectrum of forced vertical vibration self-excited forces of thin flat plate with and multi frequency

(a) 升力时程曲线

(b) 力矩时程曲线图11 薄平板扭转多频强迫振动自激力时程曲线Fig.11 The history curves of self-excited forces of thin flat plate by forced torsional vibration with multi frequency

(a) 升力时程曲线FFT幅值谱

(b) 力矩时程曲线FFT幅值谱图12 薄平板扭转多频强迫振动自激力FFT幅值谱Fig.12 The FFT amplitude spectrum of forced torsional vibration self-excited forces of thin flat plate with multi frequency

(a) 升力时程曲线

(b) 力矩时程曲线图13 大带东桥主梁断面竖向多频强迫振动自激力时程曲线Fig.13 The history curves of self-excited forces of the Great Belt East Bridge girder section by forced vertical vibration with multi frequency

(a) 升力时程曲线FFT幅值谱

(b) 力矩时程曲线FFT幅值谱图14 大带东桥主梁断面竖向多频强迫振动自激力FFT幅值谱Fig.14 The FFT amplitude spectrum of forced vertical vibration self-excited forces of the Great Belt East Bridge girder section with multi frequency

(a) 升力时程曲线

(b) 力矩时程曲线图15 大带东桥主梁断面扭转分状态多频强迫振动自激力时程曲线Fig.15 The history curves of self-excited forces of Great Belt East Bridge girder section by forced torsional vibration with multi frequency

(a) 升力时程曲线FFT幅值谱

(b) 力矩时程曲线FFT幅值谱图16 大带东桥主梁断面扭转分状态多频强迫振动自激力FFT幅值谱Fig.16 The FFT amplitude spectrum of forced torsional vibration self-excited forces of the Great Belt East Bridge girder section with multi frequency

从图10、12、14、16中均可以看出，升力与力矩系数中，各个预设的振动频率能量均较为集中，并且包含了所有预设的扭转振动频率。

利用上述数据识别的薄平板颤振导数结果与Theodorson理论解进行对比，结果如图17所示；大带动桥主梁断面颤振导数结果与耦合振动法识别的结果进行对比，如图18所示。

图17 薄平板弯扭分状态识别法颤振导数识别结果Fig.17 Identification of flutter derivatives of a thin plate with a multi-frequency step-by-step identification method

图18 大带东桥断面弯扭分状态多频识别法颤振导数识别结果Fig.18 Identification of flutter derivatives of the Great Belt East Bridge main girder section with a multi-frequency step-by-step identification method

表6 大带东桥梁断面颤振导数不同计算方法识别效率对比Tab.6 Comparison of the efficiency of different methods for identifying flutter derivatives of the Great Belt East Bridge

4 结 论

采用计算流体动力学方法对桥梁断面气动导数识别方法进行了研究，取得了如下主要研究成果：

(1)分别针对薄平板和大带东桥主梁断面，采用弯扭耦合强迫振动方法进行了主梁断面气动导数识别，识别结果分别于西奥多森理论解和已有文献结果吻合较好，验证了本文弯扭耦合强迫振动识别方法的精度。

(2)在验证了气动自激力可叠加性的基础上，提出桥梁断面气动导数识别的分状态多频强迫振动识别法，计算结果表明该方法的精度与传统弯扭耦合识别方法精度相当，而计算效率得到较大的提升。