0 引言

涡激振动是桥梁结构一种典型的自激、限幅非线性振动现象，会引起结构疲劳，并影响行车舒适性。随着我国大跨度桥梁建设的越来越多，研究与防止涡激振动已成为我国大跨度桥梁抗风设计的关键问题之一。

近年来随着计算流体动力学的发展，结构涡激振动数值模拟研究受到越来越多学者的关注。Nomura采用基于任意拉格朗日-欧拉法计算了雷诺数分别为Re=100，400时圆柱和H形柱体的涡激振动响应^[1-2]。曹丰产等采用动网格方法考虑圆柱体与流体的耦合，采用二阶投影法和多重网格法求解N-S方程，柱体振动方程采用Newmark-β方法进行求解，实现了圆柱涡激振动响应的数值模拟^[3]。徐枫、欧进萍采用动网格技术和滑移网格技术对二维方柱进行非定常绕流与涡激振动的数值模拟^[4]。Pan等采用SST k-w湍流模型对二维圆柱进行涡振响应数值模拟。研究表明，二维数值结果不能真实反映漩涡脱落的三维效应^[5]。方平治^[6]、李永乐^[7]均采用松耦合方法对二维方柱的涡振进行了数值模拟。刘志文等采用SST k-w湍流模型和动网格技术对宽高比为4的矩形断面涡激振动响应进行了计算，数值模拟结果与试验结果吻合较好^[8]。Al-Jamal H等采用大涡模拟对雷诺数为Re=8 000时的二维圆柱进行了涡振数值模拟^[9]。Sarwar等采用大涡模拟方法分别对宽高比为4的矩形断面和箱梁断面进行了涡激振动三维数值模拟研究^[10]。Zhao等采用PG-FEM求解三维N-S方程的方式实现弹性支承三维圆柱涡激振动数值模拟，并与二维数值模拟结果进行了比较^[11]。龚慧星等通过三维LES模拟与二维RANS数值模拟对宽高比5的矩形断面的涡激振动机理进行了研究^[12]。

随着计算流体动力学的发展，开源软件的应用也逐渐增多，其中OpenFOAM(Open Field Operation and Manipulation)是计算流体动力学领域影响较大的开源软件之一，国内外许多学者采用该开源软件进行结构气动性能数值模拟研究。Robertson基于OpenFOAM对钝体断面静止绕流进行了数值模拟，并且将计算设置和结果与Fluent进行了比较，结果表明OpenFOAM计算精度尚可，但是其对求解设置相对较难，而且对网格质量要求更高^[13]。Lysenko采用OpenFOAM并行计算对圆柱断面绕流问题进行了深入研究，计算结果与试验结果吻合较好^[14]。Abe采用OpenFOAM耦合延迟分离涡湍流模型求解流体运动方程，采用外接有限元求解器求解刚体运动方程，研究了圆柱断面的涡激振动问题^[15]。Tawekal采用OpenFOAM研究了不同雷诺数下圆柱断面静止绕流和涡激振动响应数值模拟^[16]。国内部分学者采用OpenFoam进行了圆柱、方柱和桥梁断面的绕流研究^[17-19]。

综上所述，目前结构涡激振动的数值模拟以圆柱、方柱为主，部分学者也进行了矩形断面和箱梁断面涡激振动数值模拟研究。实际桥梁工程中，矩形断面的应用则较多，如大跨度桥梁主梁断面、桥塔断面等均是以矩形断面为原型发展起来的。因此，关于矩形断面的气动性能的研究是钝体空气动力学研究的热点问题，如近年来国际风工程协会发起了宽高比为5的矩形断面梁气动性能标准问题研究，简称BARC(A Benchmark on the Aerodynamics of a Rectangular 5:1 Cylinder)^[20]。本研究采用开源流体力学软件OpenFOAM对宽高比为5的矩形断面进行静止绕流和涡激振动响应数值模拟研究，为今后采用OpenFOAM进行桥梁主梁断面气动性能数值模拟研究奠定基础。

1 计算方法 1.1 控制方程

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

(1)

(2)

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

1.2 流固耦合计算方法

将矩形断面竖向涡激振动简化为弹性支承单自由度振子模型，如图 1所示。对应的结构振动控制方程为：

(3)

式中，分别为t时刻时结构竖向振动位移、速度和加速度；ξ_h为弹性支承结构竖向振动系统阻尼比；ω_h为弹性支承结构竖向振动系统圆频率；m为单位长度结构质量；F_L(t)为时刻t对应的单位长度结构所受到的升力。

采用OpenFOAM中的动网格求解更新器sixDoFRigidMotion进行动网格计算与更新，选择Newmark-β法进行结构涡激振动方程求解，获得结构振动位移、速度等，然后通过拉普拉斯扩散法来改变网格中的点坐标场从而实现网格更新，即：

(4)

式中，d为单元中心点的位移场；r为网格中心与运动边界的距离；γ(r)为关于r的函数，一般采取线性插值的方式取值。

2 静止绕流数值模拟 2.1 计算域及网格

取矩形断面宽为B=600 mm，高为H=120 mm。选取计算域为30B×20B的矩形计算域进行计算，如图 2计算域所示。计算域边界条件为：左侧、下侧为入口条件(Inlet)，右侧和上侧为出口条件(Outlet)，矩形断面为无滑移壁面条件(Wall)。

采用分块思路进行网格划分，所有网格均设为四边形网格，网格增长比率控制在1.1以内，首层网格厚度参考文献[8]拟为0.001B。整体与局部网格如图 3所示。

利用OpenFOAM中网格质量检验工具checkMesh进行检测，该网格的网格数量为33 611，网格最大畸变角为34.57°，网格最大畸变率为0.61。网格畸变率和近壁面渐变率符合Bakker的建议值，即四边形网格最大畸变不得超过0.85，近壁面渐变率不得超过1.2^[21]。

2.2 边界条件与离散格式设置

入口边界(Inlet)为速度入口，风速绝对值取5.4 m/s, 其分量值根据不同的角度进行换算。湍流强度取为0.5%，湍流黏性比取10%，雷诺数为4.4×10⁴(以断面高度为参考尺寸)。OpenFOAM中入口处湍流特性值按Saxena的建议^[22]计算：

(5)

(6)

(7)

式中，k为湍动能；U为速度值大小；I为湍流强度；w为湍流比耗散率；β为湍流黏性比；ε为湍流耗散率；ν为空气动力黏性系数；ν=m/r=1.46×10^-5 m²/s，边界处各变量的边界条件设置如表 1所示。

表中，U为速度；p为压强；nut为μ，即动力黏性系数；k为湍动能；omega和epsilon分别为湍流比耗散率(ω)和耗散率(ε)。zeroGradient为诺依曼边界条件，fixedValue为狄里克雷边界条件，kqRWall Function，omegaWallFunction，epsilon WallFuntcion表示在nut选择壁面函数时对应的各自的壁面函数方程，calculated表示OpenFOAM程序会根据其他变量进行内部计算。

值得注意的是nutUSpaldingWallFunction壁面函数采用从黏性底层连续变化到对数层y⁺-u⁺多项式，将壁面处网格落于黏性层、缓冲层或者对数层都是有效的，因此当速度变化时可不用修改网格尺寸而满足y⁺的要求，并且整体效应(如三分力)的计算结果与首层网格高度是无关的。

相比于Ansys Fluent，OpenFOAM中的离散格式设置更为苛刻，这在文献[13]中已有研究。为此，参考该文献进行离散格式设置如表 2所示。

表中，phi表示控制方程通用变量，对连续方程而言取值为1，对动量方程取值为速度U。default表示所有的该项格式使用同一种默认离散格式，除对流项离散格式需要特殊指定外，其他均为统一设置。上述中div(nuEff*dev2(T(grad(U))))表示湍流模型方程中的某一湍流项。上述格式设置不唯一，但是二维RANS对流项离散格式至少需要具备二阶精度。

2.3 静止绕流计算

在划分好网格后，首先采用OpenFOAM中势流求解器(potientialFoam)进行流场初始化，然后采用不可压缩流稳态求解器(simpleFoam)进行流场稳态求解至残差曲线平稳波动，最后将计算平稳后的稳态计算结果作为瞬态求解器(pisoFoam)初始条件进行瞬态计算。整个计算的流程图如图 4所示。

计算过程中，输出矩形断面的三分力系数以进行研究。矩形断面三分力系数和斯托罗哈数定义如下：

(8)

(9)

式中，C_D，C_L，C_M分别为矩形断面的阻力、升力及升力矩系数；F_D，F_L，M分别为矩形断面的阻力、升力及升力矩，均以风轴为参考坐标系；r为空气密度；B为矩形断面宽度；D为矩形断面高度；f_s为涡脱频率；V为来流风速；St为矩形断面斯托罗哈数。

需要注意的是，采用OpenFOAM中的pisoFoam求解器计算时，最大库朗特数(Courant number，Co)必须小于1。某个网格的库郎特数表示如下：

(10)

式中，|U|为某个网格内的速度模量；Δt为计算时间步；Δx为速度方向上的网格长度。

整个计算过程中，当监控器(pyFoam)监控的各个变量的残差曲线出现平稳波动的时候即认为稳态计算平稳，当最大库朗特数(maxCo)曲线出现平稳波动的时候，即认为瞬态计算已经平稳。

限于篇幅，本研究仅给出了0°攻角下SST k-ω湍流模型计算得出的三分力时程曲线，如图 5所示。表 3给出了对应的计算结果与已有文献以及试验结果的比较, 从中可以看出本研究计算结果精度较高。

不同风攻角(α=0°, 2°, 4°, 6°)下分别采用3种RANS湍流模型(k-ε，k-ω，SST k-ω)进行计算的计算结果如表 4~表 7和图 6所示。

从表 4~表 7及图 6中可以看出，对于阻力系数及斯托罗哈数St，3种湍流模型的计算结果总体较为接近，且与试验结果较为吻合，其中SST k-ω湍流模型的计算结果与试验结果吻合相对最好；对于升力系数和升力矩系数，3种湍流模型的计算结果都与试验结果存在一定的差异。产生这种差异的主要原因是：风洞试验中0°风攻角下矩形断面的升力系数较小，其试验测试值相对误差较大；随着风攻角的增加，矩形断面周围流场的分离与再附现象明显，数值模拟结果与试验结果存在一定的差异；因此采用数值模拟方法进行大攻角下矩形断面气动性能的数值模拟仍存在一定的不足，这与文献[23]的数值模拟研究结论是一致的。综合以上计算结果可以发现，总体而言SST k-ω湍流模型的计算精度相对较好，故采用该湍流模型进行为此选择SST k-ω湍流模型进行矩形断面0°攻角下竖向涡激振动响应计算。

3 矩形断面涡激振动数值模拟

为了便于比较，矩形断面涡激振动数值模拟时根据文献[25]来确定矩形断面相关参数，如表 8所示。

将瞬态计算稳定后的结果作为初始条件，采用pimpleDyMFoam结合sixDoFRigidMotion动网格更新器进行矩形断面涡激振动响应计算。其计算流程如图 7所示。

在上述0°攻角下采用SST k-ω湍流模型计算的静止绕流稳态场的基础上添加位移场(point Displacement)，并且在边界处中指定位移场的边界条件：壁面处为程序内部计算，进口和出口边界均设置为零。

动网格区域设置时，为防止网格变形过大而影响计算结果，将网格变形区域设置在内半径为0.5 m，外半径为1 m的圆环内，其他区域均为静止区域。动网格运动动力参数按表 7中数据选取，运动方程求解器选项Newmar k-β方法。

整个计算周期较长，且需确保每一个状态下求解均达到稳定状态。在进行稳态、瞬态计算时，时间步取0.000 1 s，并且设置每一步求解的外迭代

步为20步。为了提高计算效率，采用6核并行运算。计算结束后，通过编程处理计算过程中生成地脚本文件，提取出运动位移、速度与加速度信息。

涡激振动计算风速分别取V_∞=2.54，2.94，3.20，3.53，3.89，4.06，4.20和4.32 m/s，对应的折减风速V_red=V_∞/f_hD分别为6.88，7.86，8.70，9.56，10.54，11，11.40，11.81。限于篇幅，本研究仅给出风速为V_∞=2.54，3.20以及3.89 m/s时(对应折算风速分别为V_red=V_∞/f_hD=6.88，8.70，10.54) 矩形断面升力、无量纲振动位移响应时程和升力系数FFT变换曲线，如图 8~图 10所示。

从图 8~图 10中可以看出，当折算风速为V_red=6.88时，矩形断面涡脱频率为f_s=5.050 Hz，小于结构自振频率f_h=6.152 Hz，未发生“锁定”现象，对应的矩形断面竖向无量纲振幅约为y/D=6.18×10^-5。当折算风速为V_red=8.70时，矩形断面涡脱频率为f_s=6.176 Hz，与结构的自振频率f_h=6.152 Hz接近，即发生“锁定”现象，对应的矩形断面竖向无量纲振幅约为y/D=0.069。当折算风速为V_red=10.54时，矩形断面涡脱频率为f_s=7.470 Hz，大于结构自振频率f_h=6.152 Hz，未发生锁定现象，但出现“拍”振现象，矩形断面最大无量纲振幅约为y/D=0.004。

值得注意的是，折算风速为V_red=8.70时，升力系数的频率成份随着矩形断面振动幅值的变化而变化，如图 11、图 12所示。图 11表明，当结构振幅逐渐增大时，升力系数的频率成份单一；图 12表明，当结构振动幅值达到最大并稳定后，升力系数在最大值和最小值附近均出现了一个高阶分量成份，从其幅值谱中可以看出该分量的卓越频率为18.552 Hz，约为结构涡脱卓越频率的3倍，可能是导致结构振幅不再进一步增大的一个因素，也可能是由于数值计算模型所导致的，具体机理尚需进一步探讨。

图 13所示为矩形断面涡激振动响应计算结果与试验结果的比较。从图 13中可以看出，矩形断面涡激振动响应数值模拟结果总体趋势与试验结果吻合较好，但矩形断面涡振锁定区范围比文献[25]数值模拟结果和试验结果略有提前，最大振幅较文献[25]数值模拟结果和试验结果有一定的偏大。该文献的数值模拟通过在Fluent中植入自定义函数(User Defined Function，UDF)的方式实现。需要说明的是，采用OpenFOAM3.0开源软件进行静止断面绕流的计算结果总体与试验结果较为接近；而矩形断面涡激振动响应数值模拟结果总体趋势与试验结果一致，尚存在一定的差异。这可能是试验模型的阻尼比随振幅增大有增大的现象，而数值模拟没有考虑这种现象，是导致数值模拟偏大的一个重要因素。

4 结论

本研究采用开源软件OpenFOAM对宽高比为5的矩形断面静止绕流和涡激振动数值模拟，得到如下主要结论：

(1) 采用3种湍流模型计算得到的矩形断面阻力系数、斯托罗哈数St总体较为接近，且与试验结果较为吻合；对于升力系数、升力矩系数，3种湍流模型的计算结果与试验结果均存在一定的差异，总体而言，SST k-ω湍流模型的计算结果与试验结果吻合相对最好。

(2) 基于OpenFOAM开源软件结合动网格技术对矩形断面涡激振动进行了数值模拟，结果表明：矩形断面涡激振动响应数值模拟结果总体趋势与试验结果吻合较好，但矩形断面涡振锁定区范围比试验结果略有提前，最大振幅较试验结果有一定地偏大。

总体而言，随着开源软件OpenFOAM不断地完善，采用OpenFOAM进行桥梁主梁断面气动性能与气弹效应数值模拟研究仍具有一定的可行性。