薄平板断面气动自激力振幅效应数值模拟

刘志文，张瑞林，陈政清，2

（1.湖南大学风工程与桥梁工程湖南省重点实验室，湖南长沙 410082；2.湖南大学土木工程学院，湖南长沙 410082）

1940年，旧塔科马桥在19 m/s左右的风速下发生大幅振动，扭转振幅大约为35°，并在长时间振动下因吊杆断裂导致垮塌[1].此后，桥梁结构气动稳定性问题受到广泛关注，Scanlan等[2]将航空领域的颤振导数概念加以推广，建立了适用于桥梁主梁断面的线性颤振理论，为现代大跨度桥梁抗风设计提供指导.随着桥梁跨度的进一步增加，大跨桥梁主梁在强风作用下可能会出现大幅振动现象，传统的Scanlan线性颤振理论难以适用于桥梁主梁大幅振动情况.国内外许多学者开始关注振幅对桥梁主梁断面的颤振性能影响.

桥梁主梁断面气动力在大振幅条件下存在高次谐波分量.陈政清等[3]通过强迫振动风洞试验研究发现钝体桥梁主梁断面存在明显的高次谐波分量.Diana等[4-5]通过墨西拿海峡大桥主梁断面风洞试验研究，发现气动力存在高次谐波分量和气动力的迟滞效应.王骑[6]通过风洞试验研究发现大振幅下薄翼断面与流线型箱梁断面气动力存在显著的高次谐波分量.唐煜[7]采用计算流体动力学方法（Computational Fluid Dynamics，CFD）进行了薄平板断面与流线型箱梁断面强迫振动研究，结果表明：大振幅下气动力存在高次谐波分量，并且振幅越大，高次谐波分量越明显.熊龙等[8]、Lin等[9]分别进行了流线型箱梁断面与B/H=5矩形断面（B、H分别为矩形断面的宽度与高度）强迫振动风洞试验研究，结果表明：较大扭转振幅下流线型箱梁断面与B/H=5矩形断面气动力存在显著的高次谐波分量.Gao等[10]、Zhang等[11]分别在双边肋主梁断面与流线型箱梁断面极限环颤振中发现气动力存在二阶与三阶分量.

颤振导数是表征桥梁结构断面气动力的主要参数，为此部分学者进行了典型断面颤振导数的振幅效应研究.Noda等[12]通过风洞试验研究了宽高比B/H=13与150的薄矩形断面在不同振幅下的颤振导数，研究表明：扭转振幅对薄矩形断面颤振导数影响较大.熊龙等[8]、Lin等[9]分别进行了流线型箱梁断面与B/H=5矩形断面大振幅的强迫振动风洞试验研究，结果表明：扭转振幅对流线型箱梁断面与B/H=5矩形断面颤振导数影响较大，其中对的影响最为明显.Zhang等[11]采用CFD的方法进行了流线型箱梁断面的强迫振动研究，结果表明：箱梁断面颤振导数具有显著的振幅相关性，而对竖向扭转模态耦合因素不敏感.

综上所述，国内外许多学者开展了主梁断面气动力的振幅效应的试验与数值模拟研究，已有研究表明：大振幅下气动力存在高阶分量，而且随振幅的增大颤振导数会发生改变.但目前对桥梁断面气动力幅值与迟滞相位随振幅变化以及非线性气动力作用下颤振位移响应演变规律的研究较少.

为此，本文选取薄平板断面作为研究对象[3，7，13]，采用CFD方法进行薄平板断面气动自激力的振幅效应研究.首先采用强迫振动数值模拟方法，比较了不同振幅下的薄平板断面颤振导数与气动力迟滞相位，分析了薄平板断面气动力的频谱特性，然后采用自由振动数值模拟方法，分析了薄平板断面的颤振响应演变规律.此外，为验证本文CFD数值模拟结果的可靠性，将颤振导数与颤振临界风速结果与理想平板理论解以及已有文献试验结果进行了对比.

1 数值计算方法

1.1 流体力学控制方程

目前，在桥梁结构风工程领域，基于雷诺平均法（RANS）的数值模拟方法被广泛采用[14-15]，其原理是将瞬态方程中的场变量分解为时均值和脉动值：

1.2 求解与计算

本文采用大型商业软件ANSYS Fluent求解流场，采用剪应力输运湍流模型（SST k-w）以使流体控制方程封闭，基于有限体积法求解控制方程.空间离散采用二阶迎风格式，时间离散采用二阶隐式时间积分.基于SIMPLEC算法处理压力-速度耦合.

数值模拟计算在DELL Precision 7920 Tower工作站进行，采用4个核并行计算，网格数量约为8.57万，计算效率约为1 h/千步.

1.3 强迫振动与自由振动数值模拟方法

采用大型流体力学分析软件ANSYS Fluent动网格技术，通过动网格驱动宏（DEFINE-CG-MOTION）实现结构断面运动.对于结构强迫振动，指定结构断面相应速度即可实现结构运动状态更新.对于结构自由振动，通过求解结构振动方程得到结构振动的位移响应，再通过动网格实现结构运动状态的更新.结构在升力和升力矩作用下的振动方程为：

式中：m为单位长度质量；ξh0为竖向运动阻尼比；wh0为竖向运动圆频率；I为单位长度惯性矩；ξa0为扭转运动阻尼比；wa0为扭转运动圆频率；Lse、Mse分别为升力与升力矩.结合Fluent动网格技术，通过用户定义函数（UDF）将Newmark-β 数值算法嵌入Fluent中以求解式（4）（5），计算薄平板断面自由振动响应.

1.4 薄平板断面及计算域

薄平板断面宽度B为0.45 m，高度H为0.001 m，宽高比B/H为450 ∶1.计算域设为40B×20B（B为薄平板断面宽度）.图1所示为薄平板断面形状以及计算域与边界条件.计算域与边界条件如下：计算域左侧为速度入口边界（Velocity inlet），计算域右侧为压力出口边界（Pressure outlet），计算域上下侧为对称边界（Symmetry），薄平板断面为固定壁面边界（Wall）.

图1 薄平板断面形状、计算域及边界条件（单位：mm）Fig.1 Thin plate section，computational domain and boundary conditions（unit：mm）

1.5 网格无关性与时间步无关性验证

以扭转振幅α0=12°、折算风速U/fB=17.78的强迫振动为例进行网格无关性与时间步无关性验证.不同网格划分方案如表1所示，图2（a）所示为气动升力与力矩时程计算结果，其中时间步Δt均为0.000 5 s.由图2（a）可知不同网格方案计算得到的气动升力与力矩时程曲线几乎完全一致.对网格2进行时间步无关性验证，时间步Δt分别为0.000 3 s、0.000 5 s与0.000 8 s，图2（b）所示为气动升力与力矩时程计算结果，由图2（b）可知不同时间步计算得到的气动升力与力矩时程曲线几乎完全一致.

表1 薄平板断面网格参数Tab.1 Grid parameters of the thin plate section

图2 气动力时程曲线（α0=12°、U/fB=17.78）Fig.2 Aerodynamic time history curve（α0=12°，U/fB=17.78）

为确保计算结果可靠，同时考虑计算效率，确定网格2作为计算模型，即网格数量为8.57万，首层网格高度为10-5B，时间步Δt为0.000 5 s，近壁面网格情况如图3所示.计算结果显示y+＜1，如图4所示.

图3 近壁面网格Fig.3 Near-wall grids

图4 薄平板断面y+分布Fig.4 y+distribution of thin plate section

2 振幅对气动自激力的影响

2.1 不同振幅下颤振导数

颤振导数是表征结构断面气动自激力的重要气动参数，其与结构断面运动状态线性组合表示气动力的线性部分.在颤振导数识别方面有强迫振动方法与自由振动方法，考虑到强迫振动方法识别颤振导数具有重复性好、折算风速范围广的优点[16]，本文采用分状态单自由度强迫振动方法识别颤振导数.

风速U分别为4 m/s、6 m/s、8 m/s、10 m/s、12 m/s、16 m/s和18 m/s；雷诺数Re为1.2×105～5.5×105；扭转运动振幅α0分别取为1°、2°、4°、6°、8°、10°和12°；竖向运动振幅h0分别取为0.02B、0.08B、0.16B和0.20B；振动频率f均为2.0 Hz.薄平板断面线性气动自激力表达式为：

当薄平板断面分别作扭转强迫振动与竖向强迫振动时，对应的扭转与竖向位移分别为：

式中：α0为扭转运动振幅；h0为竖向运动振幅.

根据不同振幅条件下薄平板断面气动自激力时程，采用最小二乘法进行薄平板断面颤振导数识别，不同振幅对应的薄平板断面的颤振导数随折算风速的变化曲线如图5所示.

由图5可知，薄平板断面在小振幅下（α0=1°、h0=0.02B）的颤振导数均与理想平板颤振导数理论解[17]吻合较好，表明本文强迫振动数值模拟结果具有良好的精度.

图5 不同振幅下薄平板断面颤振导数随折算风速的变化Fig.5 Flutter derivatives of the thin plate section at different amplitudes vs.reduced wind velocity

由图5（a）～（d）可知，扭转振幅对薄平板断面的颤振导数影响较大，其中对颤振导数的影响最为明显.当扭转振幅α0=6°，且折算风速较大时，颤振导数与小振幅下（α0=1°）颤振导数开始出现差异；在高折算风速下，随着振幅增大（α0=10°、12°）颤振导数均出现了由负转正的变化，这与文献[12]薄矩形板（B/H=150）的风洞试验结果一致.当扭转振幅为α0=6°、8°、10°、12°时，随折算风速的增大颤振导数与小振幅下（α0=1°）颤振导数的差别逐渐增大.扭转振幅对颤振导数的影响较小，但在高折算风速下也有一定影响.

由图5（e）～（h）可知，竖向振幅对薄平板断面的颤振导数影响总体较小.

2.2 气动自激力的幅值与迟滞相位

薄平板断面的线性气动自激力表示为正弦函数形式，即：

式中：L0、M0为气动力幅值；φL、φM为气动力相对于位移的迟滞相位.比较式（6）～（9）与式（10）（11），可以将颤振导数用气动力幅值与气动力迟滞相位表示，如式（12）～（19）所示.由式（12）～（19）可知，颤振导数的改变可由L0/α0、M0/α0、L0/h0、M0/h0与φL、φM表征.

根据薄平板断面强迫振动模拟所得的气动力，分别对式（10）（11）进行最小二乘拟合，可得到薄平板断面的气动力幅值与迟滞相位，具体结果如图6所示.

由图6（a）～（f）可知，对于扭转运动，L0/α0、M0/α0随振幅变化较小.升力与力矩的迟滞相位的正弦值sin φLα、sin φΜα随振幅变化较大，其中sin φMα在折算风速U/fB为20.0时，其值在1°振幅时为0.28，在12°振幅时为-0.26，这也引起了颤振导数在高折算风速下，随着振幅的增大发生了由负转正的变化.升力与力矩的迟滞相位余弦值cos φLα、cos φΜα随振幅变化不明显，且其绝对值大多大于0.94，所以迟滞相位并没有引起颤振导数发生较大变化.

图6 气动力幅值与迟滞相位正弦值与余弦值Fig.6 Amplitude and hysteresis phase sinusoidal and cosine values of aerodynamic forces

由图6（g）～（l）可知，对于竖向运动，L0/h0、M0/h0随振幅变化较小.升力与力矩的迟滞相位的正弦值sin φLh、sin φΜh变化亦较小，且其绝对值大多大于0.96，所以颤振导数随振幅变化较小.升力与力矩的迟滞相位的余弦值cos φLh、cos φΜh在低折算风速下略有变化，高折算风速下变化不明显，这引起颤振导数在低折算风速下有较小变化.

2.3 气动自激力频谱特性

采用快速傅里叶变换方法（Fast Fourier Transform，FFT）对薄平板断面单自由度强迫振动得到的气动自激力进行频谱分析，图7为折算风速U/fB为17.78的薄平板断面气动力频谱图.

由图7（a）（b）可知，当薄平板断面的扭转振幅为12°时，升力与力矩均存在明显高阶分量，主要为三阶、五阶等奇数倍频分量；由图7（c）（d）可知，竖向振幅为0.20B时，气动力同样存在高阶分量，主要为三阶分量，与扭转运动相比，竖向运动引起的薄平板断面的气动力高阶分量则相对较小.

图7 气动力频谱特性（U/fB=17.78）Fig.7 Spectrum characteristics of aerodynamic forces（U/fB=17.78）

为进一步分析扭转振幅对薄平板断面气动力分量的影响，定义气动力分量幅值比Ri：

式中：Fi为i阶气动升力或力矩幅值；f为频率；Ff为频率f对应的气动力幅值.由式（20）计算不同扭转振幅下薄平板断面气动力一阶、三阶与五阶分量的幅值比，结果如图8所示.

图8 不同扭转振幅下气动力分量幅值比（U/fB=17.78）Fig.8 Amplitude ratio of aerodynamic components at different torsional amplitude（U/fB=17.78）

由图8可知，随着扭转振幅的增大，薄平板断面气动力一阶分量所占比例不断下降，当扭转振幅大于8°时，薄平板断面气动力出现较为明显的三阶与五阶分量.

3 颤振位移响应演变规律

3.1 颤振临界风速

采用自由振动方法进行薄平板断面颤振稳定性的直接计算[18-19]，薄平板断面自由振动计算特征参数采用文献[1]中的相关参数，如表2所示.

表2 薄平板断面自由振动计算特征参数Tab.2 Calculating dynamic parameters of free vibration of the thin plate section

颤振临界状态附近不同计算风速下薄平板断面位移时程曲线如图9所示.薄平板断面颤振临界风速与频率的计算结果及其与频域理论解如表3所示.由表3可知，CFD计算误差较小，表明本文自由振动数值模拟结果具有良好的精度.

图9 位移时程曲线（U=16.20、16.34、16.50 m/s）Fig.9 Displacement time history curve（U=16.20，16.34，16.50 m/s）

表3 薄平板断面颤振临界风速与频率Tab.3 Critical flutter wind velocity and freguency of the thin plate section

3.2 颤振位移响应演变规律

非线性气动力作用下，结构振动的频率和阻尼比不是一个常数，而是随振幅变化，可根据薄平板断面自由振位移响应时程计算得到：

式中：Q（ti）为第i个位移峰值；ti为第i个位移峰值对应的时间；ξ（t）为瞬时阻尼比；ω（t）瞬时圆频率；φhα（t）为竖向与扭转振动位移相位差.另外，瞬时阻尼比为系统阻尼与气动阻尼两者共同作用的结果.

图10为薄平板断面在风速U为17.0 m/s（U/fa0B=12.49）时自由振动位移响应时程曲线.根据位移响应时程由式（21）～（23）计算得到薄平板断面竖向与扭转振动对应的瞬时频率、瞬时阻尼比以及竖向与扭转振动位移相位差随时间变化曲线，如图11所示.

图10 位移时程曲线（U=17.0 m/s）Fig.10 Displacement time history curve（U=17.0 m/s）

图11 瞬时频率、瞬时阻尼比与竖向扭转位移相位差（U=17.0 m/s）Fig.11 Instantaneous frequency，instantaneous damping ratio and phase difference between the vertical displacement and torsional displacement（U=17.0 m/s）

由图11（a）可知，竖向与扭转瞬时频率在小振幅（α0＜1.44°）下不随振幅变化，且两者为同一数值.随着振幅的增大，两者均呈现了减小的趋势，在减小过程中竖向瞬时频率略大于扭转瞬时频率，换言之，是扭转瞬时频率下降略早于竖向瞬时频率.由于薄平板断面为弯扭耦合颤振，竖向与扭转频率变化趋势具有一致性，但在变化过程中两者频率不是同一个数值.值得说明的是瞬时频率存在一定跳跃性，这是由时间离散引起的数值误差.

由图11（b）可知，竖向与扭转瞬时阻尼比绝对值随着振幅增大而不断增大，说明颤振发散速率越来越大.这与颤振导数随振幅变化有关，随着扭转振幅的增大颤振导数逐渐变小，再变为负值，导致随着扭转振幅的增大颤振导数提供的气动正阻尼逐渐减小，再变为提供气动负阻尼.薄平板断面颤振为竖向与扭转耦合的振动形式，所以竖向与扭转瞬时阻尼比的绝对值均随振幅的增大而增大，且竖向与扭转瞬时阻尼比基本保持一致.

由图11（c）可知，扭转振幅α0＞1.44°时，竖向与扭转位移相位差随着振幅增大不断增大.竖向位移相位超前于扭转相位，又由于竖向瞬时频率略大于扭转瞬时频率，导致竖向与扭转位移相位差进一步增大.

4 结论

以薄平板为研究对象，采用计算流体动力学方法进行强迫振动与自由振动研究，分析了气动自激力的振幅效应，得到如下主要结论：

2）当扭转振幅大于8°时，薄平板断面气动力存在较为明显高次谐波分量，主要为三阶与五阶分量；竖向振幅引起的薄平板断面气动力高阶分量成分相对较小.

3）薄平板断面颤振发散过程中，竖向与扭转瞬时频率在小振幅（α0＜1.44°）下基本不随振幅变化，随着振幅进一步增大两者均呈现减小的趋势；竖向与扭转瞬时阻尼比绝对值随着振幅增大而不断增大；扭转振幅α0＞1.44°时，竖向与扭转位移相位差随着振幅增大不断增大.