基于两种湍流模型的桥梁颤振导数识别研究及比较＊

祝志文，夏 昌

（湖南大学风工程试验研究中心，湖南长沙 410082）

风洞试验是研究大跨度桥梁颤振稳定性的主要手段.风洞试验简单直观，但风洞试验涉及模型制作、复杂的仪器设备，并且其试验周期长，费用高，紊流风场的有效模拟比较困难.随着计算机技术的快速发展，采用数值计算的方法来识别颤振导数成为可能，Larsen的离散涡法［1］，Zhu等开发的有限体积法［2］，曹丰产等开发的有限元方法［3］，都取得了重要的成果.其中部分方法还是在均匀流场中开展的，并没有实现对湍流的模拟，虽然计算工作量减小，但网格等可能并没有达到要求.大型CFD商业软件FLUENT和CFX等以其通用性和便利的后处理功能受到数值风洞研究者的肯定，同时提供了多种湍流模型，因此它们在结构风工程领域的研究中越来越受到重视.

标准k－ε湍流模型是由Launder和Spalding［4］于1972年首先提出的，该湍流模型属于高Re数模型，在近壁区并不直接求解，而是通过壁面函数将近壁区的变量与湍流核心区变量联系起来求解，因此存在固有近边界模拟的缺陷，同时在此模型中假设湍流黏性系数μt是相同的，而在流线弯曲和应变率大的情况下，湍流是各向异性的，因此，此模型用于强旋流、弯曲壁面流动或弯曲流线流动时，会产生一定的失真［5］.Menter在1993年对Wilcox提出的标准两方程k－ω湍流模型进行修正后提出分区剪应力输送（Shear Stress Transport）模型［6］，即SST湍流模型.该模型实际上是从边界层内部的标准k－ω模型逐渐转变到边界层外部的高Re数的标准k－ε模型，在近壁区，它不再引入如壁面函数这样的经验公式，而是采用标准k－ω模型直接求解，能有效地使模型流线分离；在湍流核心区采用标准k－ε模型计算，同时，对湍流黏性系数μt进行修正，这样既提高了数值模拟的精度，又有效减小了计算工作量.

过去由于受到计算条件的限制，标准k－ε湍流模型得到广泛应用，同时人们根据需要提出了多种相关的修正模型，有力地推动了数值风洞的发展.近几年来，随着计算速度的大幅提高，研究者们才开始重视SST模型的使用，基于此模型利用数值计算的方法，成功地计算了桥梁断面的三分力系数，数值计算值与试验值吻合较好，体现了SST湍流模型的计算优势.采用数值方法识别颤振导数时，研究者［7］常用标准k－ε湍流模型，采用SST湍流模型的相关报道却很少见.

采用数值模拟桥梁断面强迫振动，需要采用动网格技术.研究者［7］提出采用断面外包矩形刚性网格的方法，以此提高数值模拟精度.但在矩形直角周围的动网格质量却不理想（疏密不均），由此可能引起在气动力时程曲线中，出现大量“毛刺”.本文采用外包椭圆形刚性网格的形式，既能保证刚性网格区的质量，又能兼顾动网格的质量.

本文CFD计算基于大型商业软件FLUENT平台开展，湍流模型分别采用FLUENT提供的基于雷诺平均（RANS）的标准k－ε模型和SST模型，分析两种湍流模型的使用条件，尝试采用SST湍流模型识别桥梁断面颤振导数，同时比较两湍流模型在识别颤振导数上的计算精度.

1 控制方程及湍流模型分析

基本控制方程为连续方程与雷诺平均N－S方程（RANS）：

式中：ρ，μ分别表示流体的密度和动力粘度；ui，uj代表某个方向上的平均流速，u'i为速度分量的脉动量，对于二维问题，i与j的取值范围为1，2；－项定义为Reynold应力［8］，这属于新的未知量，为了使方程组可以封闭，故引入湍流模型（雷诺平均N－S方程），以便求解.

标准k－ε模型是高Re数的湍流模型，它只能求解湍流核心区的流动，直接来求解桥梁断面近壁区的流动是不正确的.因此，使用此湍流模型时需要引入壁面函数，对湍流核心区的流动使用k－ε模型求解，对近壁区不进行求解，而是直接采用一组经验公式，将近壁区的变量与湍流核心区的求解变量联系起来求解.采用标准k－ε模型时，只要把第一层网格节点布置在湍流充分发展的区域就可以了，并不需要在桥梁断面近壁区加密网格，因此减少了网格数量，节省了数值计算的时间.当计算机计算速度有限时，此模型显示出明显的优势.为了合理地使用壁面函数，需要对断面的Yplus值进行控制，一般要求满足Yplus≥11.63［5］.

SST湍流模型是在k－ω模型基础上发展而来的，且融合了k－ε模型的优点.k－ω湍流模型的优势是在低Re数下的近壁计算，k－ε模型适合湍流核心区的计算.SST模型克服了标准k－ω湍流模型对自由流参数变化比较敏感的缺点，在近壁区采用k－ω湍流模型，在远离壁面的流场中采用k－ε湍流模型.这样充分利用了k－ε湍流模型对逆压梯度流动具有较高的模拟精度和k－ε湍流模型对湍流自由流参数不敏感的优点.采用SST模型时，需要对近壁区进行数值计算，而不是采用类似的壁面函数，因此需要加密近壁区的网格，同时要合理控制第一层网格高度，大体上使得Yplus≤6［8］.

2 气动自激力描述

气动自激力的Scanlan表达式为：

式中：h（t），α（t）分别为竖弯与扭转位移；Κ＝ωB／U为无量纲频率；?（t），（t）分别为竖弯速度与扭转速度；和（i＝1，2，3，4）即为桥梁断面颤振导数，是桥梁断面外形和折算风速Ur的函数（Ur＝U／fB）.

本文通过数值模拟的方法分别提取每个折算风速下的气动自激力，然后采用最小二乘法识别相应折算风速下的8个颤振导数.

3 算 例

3.1 研究对象及计算域网格划分

本文选取丹麦大带东桥为研究对象，该桥为主跨1 624m的三跨连续钢箱梁悬索桥，其主梁是上下游侧带风嘴的扁平箱梁断面（如图1所示）.箱梁截面全宽B＝31m，高宽比为7.05∶1.CFD计算的模型采用与风洞试验相同的几何缩尺比1∶80，数值模拟时不考虑栏杆和防撞墙等附属物.计算域采用二维圆形区域，左侧半圆弧为来流进口，到模型中心的距离为30B，右侧半圆弧为来流出口，到模型中心的距离为30B.

图1 大带东桥主梁断面（单位：m）Fig.1 Girder section of the Great Belt Bridge（unit：m）

由于CFD模拟时刚性桥梁断面在每一时间步上运动，因而在每一时间步上需要重新对计算域网格进行划分.为确保桥梁最大运动位移处有较好的网格质量，不至出现畸变网格甚至负体积网格，本文将计算域进行分区划分，分成3个大小近似合理的区，如图2所示，并采用不同的网格进行剖分.围绕桥梁的称为刚性网格区.桥梁断面运动时，该区域网格与桥梁断面刚性固定，并在每一时间步上与桥梁断面同步运动.该区域外边界为椭圆，通过对该椭圆适当分段，便于对该域进行四边形结构网格划分，以便能获得较好的正交网格.计算域绝大部分区域称静止网格区，该区域外边界是计算域外边界，内边界为离开刚性网格区外椭圆一定距离，且包围刚性网格区的圆形.静止网格区采用四边形单元剖分，从内到外采用合适的网格放大比例.静止网格区和刚性网格区在整个CFD模拟过程中一直使用计算开始时的网格系统，不再重新划分网格.在静止网格区和刚性网格区之间为动网格区，动网格区采用三角形单元剖分.在每一时间步上，该区域根据桥梁断面的运动位置并由设定的网格系统质量要求重新进行网格划分.紧靠桥梁断面的区域流场变量变化剧烈，特别是断面迎风侧和断面法向，因此网格划分需要能适应流场变量的变化程度，并沿各个方向采用适度的网格放大率，实现与动网格区域网格的平顺过渡，如图3所示.

图2 计算域分区Fig.2 Computational domain decomposition

值得注意的是，为了分别满足本文所采用两湍流模型的使用要求，需要划分两套网格，如图3所示.这两套网格计算域划分形式完全一样，只有刚性网格区的第一层网格高度和网格数量不一样.因为需要通过控制第一层网格的高度来调节断面的Yplus值，本文中图3（a）第一层网格高度取为2mm（约为0.005B），而图3（b）第一层网格高度取为0.25mm（约为0.000 6B）；两套网格的网格数量分布见表1.两套网格在数量上比较，只有刚性网格区的网格数量不同，这是由于使用SST模型时需要在近壁区加密网格.两套网格中，网格数量都集中在动网格区，这是因为经过反复试算表明：适当增加动网格区的网格数量，可以有效消除气动力系数时程曲线中的“毛刺”.

表1 计算域各区网格数量Tab.1 Meshes in each region for models

计算区域的网格划分在Gambit中实现，然后分别导入到Fluent中进行数值计算.根据上面对两湍流模型使用要求的分析，分别控制Yplus值.由图4可知，断面Yplus值符合两湍流模型的使用要求（见第1部分分析）.

图3 桥梁断面近壁网格及局部Fig.3 Girds arrangement around bridge section and close view

图4 桥梁断面Yplus值分布Fig.4 Yplus distribution around bridge section

3.2 计算条件

3.2.1 桥梁断面运动模式

采用单自由度单频等幅正弦位移激励桥梁断面运动.对纯竖弯运动，有扭转自由度α（t）＝0，桥梁断面竖弯运动位移是：

对于纯扭转运动，竖弯运动位移h（t）＝0，桥梁断面扭转运动位移是：

式中：fh，fα分别表示竖弯运动频率和扭转运动频率；h0，α0分别表示竖弯运动和扭转运动幅值，本文统一h0和α0的取值，h0取为0.025B，α0取为3°，小幅振动以满足线性小扰动假设.

3.2.2 边界条件和其他相关参数

计算域左侧进口为模拟大气边界层速度和紊流度的速度入口，紊流度为5%；计算域右侧为出流边界条件，对应沿出口边界法向速度梯度为零.识别不同折算风速下的颤振导数，本文通过改变强迫振动频率fh和fα来实现无量风速Vr＝U／fB的改变，这样保证了不同折算风速下的Re数相同，采用的非稳态计算时间步长为0.005s.

本文采用了FLUENT软件中的动网格技术，主要参数设置有：采用Smoothing和Remeshing两种动网格更新方法；网格光滑更新迭代次数设为200，弹性常数因子和边界节点松弛都设为0.6；局部网格重划分中网格最大畸变控制为0.4，网格尺寸重划分迭代次数设为100.

4 速度场比较

为了更直观地说明两湍流模型的计算特点和比较两湍流模型的计算精度，本文给出桥梁断面近壁区在nT时刻处的速度矢量图，T为模型强迫振动周期，nT时刻表示模型回到平衡位置的时刻.在模型强迫振动过程中，没有看到明显的漩涡脱落；同时对气动升力进行频谱分析，可以看到在主频率中，只有模型进行周期运动的频率.

如图5（b）所示，使用SST湍流模型时，在桥梁断面上下缘流线分离处成功地捕捉了涡，且在近壁区处，速度大小沿断面法线方向呈现梯度变化，图5（a）却无法看到这些现象.这可以解释为：标准k－ε模型在近壁区采用了壁面函数，在模拟近壁区绕流、流线分离和涡的形成时，出现了一定程度的失真.SST湍流模型在计算近壁区流动时，采用了标准kω模型，由于它属于一种低Re模型，因此有效地提高了近壁区的计算精度.可见，在模拟桥梁断面绕流和流线分离时，SST湍流模型的计算优势明显.

图5 近壁区速度矢量分布前缘及局部放大Fig.5 Velocity vector around bridge section and close view

5 气动导数识别结果

本文采用最小二乘法识别了基于两湍流模型下不同折算风速的颤振导数，并将Poulsen［9］的试验值同时列入以作参考，如图6和图7所示.由于该试验值是在均匀流中通过桥梁模型自由振动实现识别获得，而本文则考虑了紊流场，因而与本文结果没有严格数量上的可比性，但可作为定性参考.采用数值模拟的识别结果与试验值趋势大体一致（除外），气动导数，，，和的数值模拟值与试验值相差不大，符合实际工程要求，因此，基于CFD商业软件FLUENT，采用标准k－ε湍流模型和SST湍流模型来识别大带东桥的气动导数，都能取得较好的效果；本文采用SST湍流模型成功地预测了和曲线的趋势，相关文献中对这些颤振导数的识别却不理想，因此相比而言，SST湍流模型的数值模拟精度优于标准k－ε湍流模型.

图6 计算的大带东桥颤振导数（1）Fig.6 Simulated flutter derivatives of the Great Belt Bridge（1）

图7 计算的大带东桥颤振导数（2）Fig.7 Simulated flutter derivatives of the Great Belt Bridge（2）

6 结 论

本文针对基于标准k－ε湍流模型和SST湍流模型的大跨桥梁颤振导数数值识别及比较研究，得到下述结论：

1）本文提出的计算域分区网格划分方法，以及为保证内部刚性网格区网格正交性的椭圆外形边界，提供了桥梁断面绕流近流场模拟的较好计算网格，为标准k－ε湍流模型和SST湍流模型较准确地识别桥梁断面的大部分颤振导数提供了保障.

2）虽然标准k－ε湍流模型和SST湍流模型均能较准确地识别桥梁断面的大部分颤振导数，但在个别颤振导数识别上，标准k－ε湍流模型无法给出趋势性的结果，表明SST湍流模型比标准k－ε湍流模型具有明显的优势.

3）不同的湍流模型有不同Yplus要求，网格尺度和布置必须满足各个湍流模型的相应要求.

4）本文考虑了紊流场对颤振导数的影响，大部分导数数值模拟结果略大于试验值，但紊流场中的各个参数（包括紊流强度、积分尺度等）对各个颤振导数的具体影响，还必须结合相关的风洞试验和数值模拟进行精细化研究.

［1］ LARSEN A，WALTHER J H.Aeroelastic analysis of bridge girder sections based on discrete vortex simulations［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997（67／68）：253－265.

［2］ ZHU Zhi－wen，GU Ming，CHEN Zheng－qing.Wind tunnel and CFD study on identification of flutter derivatives of a long－span self－anchored suspension bridge［J］.Computer－Aided Civil and Infrastructure Engineering，2007，22：541－554.

［3］ 曹丰产，项海帆，陈艾荣.桥梁断面颤振导数和颤振临界风速的数值计算［J］.空气动力学学报，2000，18（1）：26－33.CAO Feng－chan，XIANG Hai－fan，CHEN Ai－rong.Numerical assessment of aerodynamic derivatives and critical wind speed of flutter of bridge decks［J］.ACTA Aerodynamic Sinica，2000，18（1）：26－33.（In Chinese）

［4］ LAUNDER B E，SPALDING D B.Lectures in mathematical models of turbulence［M］.London：Academic Press，1972.

［5］ VERSTEEG H K，MALALASEKERA W.An introduction to computational fluid dynamics：the finite volume method［M］.New York：Wiley，1995.

［6］ MENTER F R.Two－equation eddy－viscosity turbulence models for engineering applications［R］.New York：AIAA，1993.

［7］ HUANG Lin，LIAO Hai－li.Numerical simulation for aerodynamic derivatives of bridge deck［J］.Journal of Wind Engineering and Industrial Aerodynamics，2009，17（4）：719－729.

［8］ LIAW Kai－fan.Simulation of flow around bluff bodies and bridge deak sections using CFD［D］.Nottingham：University of Nottingham，School of Civil Engineering，2005.

［9］ POULSEN N K，DAMSGAARD A，REINHOLD T A.Determination of flutter derivatices for the great belt bridge［J］.Journal of Wind Engineering and Industrial Aerodynamics，1992，41（1／2／3）：153－164.