基于脉动风作用下的大跨度桥梁拉索风场数值模拟

孙亭亭,杨吉新,黎建华,周兴宇



基于脉动风作用下的大跨度桥梁拉索风场数值模拟

孙亭亭1,2,杨吉新2*,黎建华2,周兴宇3

1. 安庆职业技术学院 建筑工程系, 安徽 安庆 246003 2. 武汉理工大学 交通学院, 湖北武汉 430063 3. 中国市政工程中南设计研究总院有限公司, 湖北 武汉 430010

采用Kaimal脉动风速谱，利用MATLAB模拟脉动风的瞬时风速时程并导入计算流体力学软件中，选择RNG-模型，考虑在高雷诺数的作用下，对大跨度桥梁的拉索风场数值模拟。采用结构性网格建立流场模型，分别选取30 m和80 m高度处分析流场和风压力系数的变化；选取50 m高度处分析速度矢量；选取80 m高度处分析流速和迹线变化；选取拉索中心位置、拉索尾端分析拉索周围压力分布，并进行数值模拟。计算结果表明，对比文献中风洞试验的结果，流场中流体的分布特征符合类似研究的结论，可以得到合理精度的数值解；在拉索周围上、下表面形成涡街，涡流沿中轴线运动并呈现对称的形状；拉索表面压力分布和文献中风洞试验结果较吻合，并且符合工程实际。

Kaimal脉动风速谱; 三维风场; RNG-模型; 数值模拟

在大跨度桥梁中斜拉桥是常用的桥型，斜拉索是斜拉桥中重要的结构组成，斜拉索的风荷载在全桥风荷载中占得比重比较大，所以在风荷载作用下斜拉索周围流场及其响应的问题研究就显得尤为重要[1]。从以往的研究结果来看，大跨度桥梁风荷载作用主要采用风洞试验的方法来检验斜拉索和桥梁的抗风性能[2,3]。比如在苏通长江公路大桥的建设中，全桥有272根斜拉索，其中最长的斜拉索长度是577 m[4]。研究表明在横桥向风的作用下，斜拉索产生的位移和内力对于主梁位移来说占整个风荷载的60%~70%[5,6]。目前，关于斜拉索风荷载的计算依据主要有欧洲规范、日本规范和我国的JTG/T D60-01-2004《公路桥梁抗风设计规范》，同时这些规范也被作为风洞试验的依据[7]，本文在计算中的参数主要采用我国的JTG/T D60-01-2004《公路桥梁抗风设计规范》[8]。风洞试验对于研究斜拉索以及斜拉桥抗风具有重要意义，然而风洞试验造价昂贵，时间周期较长[9,10]。本文采用计算流体力学CFD数值风洞的方法来对斜拉索进行研究并结合相关文献的风洞试验结果进行分析[11]。

现在流体力学研究方法主要包括理论分析、数值计算、试验研究和现场实测四个方面[12]，这些研究方法针对不同的角度相互补充。理论分析研究主要是从各参数影响出发，为数值计算和实验研究提供有效的理论指导[13]；试验研究是验证客观科学规律的一种有效手段，用来验证理论分析和数值计算的正确性；现场实测是有效收集现场数据的一种方法，能够反映自然界中的真实情况，但是由于各种环境因素的影响，测量往往有困难，甚至是无法实施[14]；数值计算主要通过计算流体力学的理论来模拟实际的情况，包括实际环境中的复杂情况，从而补充理论及试验的空缺，更重要的是提供了廉价的模拟、设计和优化的工具，以及提供了分析三维复杂流动的工具[13,15]。从上世纪80年代以来，计算流体力学(Computational fluid dynamics CFD)发展异常迅速，结构风工程的研究除了继续采用风洞试验的传统方法以外，研究者开始采用计算流体力学的技术并使用数值模拟的方法研究结构在流场中的作用[16,17]。在数值风洞分析中高雷诺数下的流场变化无论对于理论研究还是工程实践都具有重要意义[18]，同时也是风工程中需要解决的关键问题之一[19,20]。本文利用CFD数值模拟的方法来研究斜拉桥拉索在高雷诺数作用下采用RNG-模型拉索表面压力分布和不同高度处流场的分布。

1 计算理论分析

1.1 脉动风速谱

在工程实际应用中，通常认为风速包含平均风和脉动风两部分，风荷载的模拟主要是针对脉动风而言[21]。脉动风有两个主要概率特性：脉动风的风速谱和相干函数[11]。本文在风荷载计算中主要采用Kaimal脉动风速谱[21]。

在桥梁风工程中，气流相对于桥梁而言是湍流充分发展的状态，而Kaimal脉动风速谱考虑了大气湍流运动中的湍流功率谱随高度的变化，所以采用Kaimal脉动风速谱能够更好的反应出风场的真实情况。

1.2 脉动风荷载的模拟

模型参数，设[()]={()]=[1(),2(),3(),…,u()]，u()为第点处的风速序列为时间序列，可以知道，根据自回归模型特点，第个变量不仅与本身前一个或几个时间阶段的风速u(-D)有关，而且与其他位置处同时段的风速u()和前一个或几个时段的风速u(-D)有关，此外还受到随机因素的影响[22]。写成数学表达式如下：

式中[(-D)]=[u(-D)]，是与时刻相关的个风速，是自回归模型的阶数。

采用Kaimal脉动风速谱，对所模拟出的风速时程离散值做数理统计，得到的模拟功率谱与目标功率谱吻合[23]。模型高度共分10个部分，选取第3部分和第8部分作为研究对象，通过MATLAB软件基于AR模型方法模拟拉索多变量互相关水平脉动风速时程[20]，对此斜拉索进行结构风场模拟，提取脉动风场湍流强度。将此拉索在竖直方向上分为10段，每10 m一段，斜拉索的10个点编号从下往上依次为1~10。

图1为3号节点和8号节点瞬时风速时程的模拟结果，模拟时间步长为0.1 s，共设置2048步。通过图1观察可以看出：脉动风为完全不规则的随机脉动，且同一时刻拉索不同位置的风速也不尽相同。脉动风速在-10 m/s~10 m/s之间浮动，节点瞬时风速在30 m/s~50 m/s，体现风速在时间及空间位置上的随机性。

图 1 3、8号节点瞬时风速时程

图 2 节点脉动风速模拟谱与目标谱

图2表示根据选取上述节点模拟结果，利用快速傅里叶变换(FFT)求得模拟风速和Kaimal脉动风速谱的对比。这里采用双对数坐标轴形式来表示，可见两者相当吻合[22,23]。

2 气动力和流场特性分析

2.1 相关参数分析

本算例以圆柱体作为计算模型，模型本身定义的坐标系为体轴坐标系，根据风向坐标定义的坐标系为风轴坐标系。采用RNG-模型进行数值模拟，RNG-模型来源于严格的统计技术。它和标准-模型很相似，但是在计算中有效改善了计算精度，在湍流漩涡中精度比标准-模型高，为湍流Prandtl数提供了一个解析公式[24]。这些特点使得RNG-模型比标准-模型在更广泛的流动中有更高的可信度和精度。所以本文选择RNG-模型对拉索的风场进行分析，使得数值计算的结果和实际工程更加符合。

由于缺乏风洞试验条件，为验证数值风洞模型计算准确性，采取该建模思想对同济大学土木实验室的试验模型进行建模，选取拉索倾角=90°，风速为35 m/s的工况，求解此工况下的拉索风阻系数，与物理风洞试验进行对比，对比结果如表1所示[25]。将计算结果与试验结果对比发现，采用此计算模型对流场部分进行模拟是可行的[25]。

表1 数值计算与试验结果对比

2.2 控制方程

计算风工程中，钝体绕流问题的控制方程是粘性不可压Navier-Stokes方程，基于雷诺平均的控制方程可写为[18-20]：

3 几何建模与网络划分

拉索尺寸高度100 m直径D是0.15 m，在建模时采用三维仿真实体建模，模型尺寸3.75 m×3 m×100 m。计算区域为一个长方体：上游来流为10 D，下游15 D，宽20 D，高度为100 m，其中拉索直径D为0.15 m，其中拉索的材料密度为8400 kg/m3（考虑防护材料质量），弹性模量为195 GPa。

在数值计算中网格的质量直接影响到计算结果的精确性和计算时间的长短，同时网格的质量也是作为工程应用的有效工具所面临的关键问题之一[14]。在计算敏感区(壁面附近、尾流区、外形曲率大的表面处)参数变化梯度大的区域，对这些区域需特别注意，选取较密网格。按照本文模型的特点，采用前处理软件Icem中的O-Block型划分，并对拉索周围流场进行加密，流场采用结构化网格，计算分析结构化网格比非结构化网格计算精度较高，共建立1210824个单元，经检查网格质量较好。网格划分如图3（a，b，c）所示。

图 3 网格划分示意图

a:总体模型，上顶面和下底面 Overall model, top and bottom surface；b:整体网格划分Overall mesh generation；c:拉索周围加密的网格区域 Encrypted mesh area around the cable

4 CFD数值结果讨论与分析

4.1 风压分析

以风压系数表示拉索周围的压力分布。

为相对压力（风压—标准大气压力），正值为压力，负值为吸力。

本文计算考虑在高雷诺数2.4×105下的风场的变化，计算采用RNG-模型。分别考虑在拉索30 m和80 m高度处建立切面，研究在各个不同高度处拉索表面风压的变化和流场的压力分布如图4和5所示：

图 4 30 m处拉索周围压力分布和30 m高度处流场中压力分布

Fig.4 Pressure distribution around the cable and the flow field at 30 m

图 5 80 m处拉索周围压力分布和80 m高度处流场中压力分布

图4中(a)显示在30 m高度时的拉索表面风压分布，区间在（-2.121~2.406），(b)表示30 m高度处流场压力的分布；图5中(a)显示在80 m高度时的拉索表面风压分布，区间在（-2.745~1.768），(b)表示80 m高度处流场压力的分布。

对比不同高度处的风压力分布可以看出随着拉索高度的不断变化在拉索迎风面压力不断减小，背风面压力不断增加（图6）。

从压力系数和压力云图中可以看出，拉索周围压力场分布的情况，拉索迎风面为正压区最大值1500 KN，两侧为负压区最大值2000 KN，拉索末端压力呈现-500 KN~0 KN的变化，这是由于流体绕物体流动产生边界层分离，形成涡旋，涡旋消耗来流能量，因而尾流中压力降低所引起的。

图 6 不同高度处压力系数分布

4.2 流场分析

对于圆柱体，由于空气的粘性效应，在表面会形成两个区域：一个是离表面很近的边界层区，另一个是主流区，其空气粘性摩擦力很小，基本可以忽略不计。当主流区内存在逆压梯度使气流加速或减速时，这种压力差也会影响边界层区流体的流动。

图 7 Z=50 m时速度矢量图

图 8 Z=80时流速分布图和迹线图

选取50 m高度处的速度矢量图进行分析，如图7所示，由速度矢量图可见，在拉索周围上、下表面脱落的涡沿尾流区的中轴线运动，并且在中轴线附近存在着一个上、下交界区，并保持着一定的惯性持续向流场下流运动，形成一个完整的下流涡街。

选取80 m处的流场切片显示索端部流场，如图8所示。在边界层中的流体将会受到更为强烈的阻滞，当能量消耗完后会被迫折回发生倒流，而主流能量相对较大会继续前进。于是越来越多的流体在边界区和主流区堆积，形成一个或多个旋涡。由流线图可见，因为在计算中设置的计算域关于x轴对称，所以在涡流的形状、拉索表面的压力分布、拉索端部和拉索尾部的压力分布也呈现对称的形状。

依照建模坐标选取=0、=0.1、=0.025和=0.068处的位置作为切面，包含了拉索的中心位置和尾端，通过数值计算分析各表面的压力分布，如图9所示。从计算结果中可以看出各表面的压力等值线与文献[26]风洞试验结果比较吻合，Fluent中RNG-模型对钝体结构迎风面的模拟结果较为理想。

图9 拉索周围压力分布

5 结论

本文主要通过数值模拟的方法分析大跨度桥梁拉索在脉动风作用下三维流场，得到以下结论：

（1）采用Kaimal脉动风速谱进行计算，考虑在平均风速的作用下，利用MATLAB模拟脉动风在参考点的瞬时风速时程，并采用快速傅里叶变换(FFT)模拟风速和Kaimal脉动风速谱的对比，使得模拟结果更符合工程实际；

（2）为了提高计算精度，对构建的三维流场进行结构性网格的划分，并采用RNG-湍流模型进行数值模拟，计算结果和风洞试验的结果吻合并能得到合理精度的数值解，说明选择RNG-湍流模型能够反应出工程实际的风场状态；

（3）对拉索进行分段研究，分别计算在拉索30 m和80 m高度处建立切面，比较各个高度的风压力系数。拉索迎风面正压力最大值1500 KN，两侧为负压力最大值2000 KN，拉索末端压力呈现-500 KN~0 KN的变化；

（4）通过数值模拟出的速度矢量图和流线分布图可以看出，在拉索周围上、下表面脱落的涡沿尾流区的中轴线运动。在气流经过拉索壁面时，在壁面摩擦力和压力梯度的共同作用下，使得拉索壁面的空气质点的运动受阻，质点的动量不断消耗，不能维持原来的运动，一些质点无法沿壁面向前运动，出现回流；

（5）选择拉索的中心位置（关于轴对称）进行研究，通过数值计算分析各表面的压力分布，从计算结果中可以看出涡流和压力呈现对称的分布，是因为拉索在竖向平面内周期性运动，拉索壁面上下分离的剪切层之间的相互作用形成涡街；

（6）对拉索的尾端以及下游区域进行研究，从计算结果中可以看出随着拉索高度的不断增加，拉索表面及周围的压力梯度不断增大，旋涡的尺寸也发生变化，边界层出现层流分离，其尾流逐渐变成紊流涡街。

[1] 张志田,张显雄,陈政清.桥梁气动力CFD模拟中湍流模型的应用现状[J].工程力学,2016,33(6):1-8

[2] 陶天友,王浩.基于Hermite插值的简化风场模拟[J].工程力学,2017,34(3):182-188

[3] 曾要争,喻梅.大跨度独塔连续刚构斜拉桥抗风性能研究[J].公路交通技术,2015(6):55-59

[4] 张挣鑫.斜拉桥拉索参数振动的半主动控制及MR阻尼器优化布置研究[D].长沙:中南大学,2013

[5] Zheng Z, Lei J. Application of the-Re θ, Transition Model to Simulations of the Flow Past a Circular Cylinder[J]. Flow, Turbulence and Combustion, 2016,97(2):401-426

[6] 吕东阳.大跨度屋盖结构表面风压的数值模拟研究[D].西安:长安大学,2015

[7] 张建,李天,杨庆山.雷诺数对大跨屋盖结构表面风压的影响[J].哈尔滨工业大学学报,2016,48(6):38-42

[8] 何艾狄.基于模态参数识别的高耸结构五点等效风荷载识别技术[D].武汉:武汉理工大学,2013

[9] Sarkar S, Ganguly S, Biswas G,. Effect of cylinder rotation during mixed convective flow of nanofluids past a circular cylinder[J]. Computers & Fluids, 2016,127(C):47-64

[10] 方治华,李晨.斜拉索上瞬时风速与风荷载的数值模拟[J].振动与冲击,2010,29(7):210-212

[11] 徐丰,詹昊,梁琛.基于自回归模型的桥梁脉动风场模拟[J].武汉工程大学学报,2015,37(3):25-28

[12] 周彬彬,蔡建国,冯健.基于ARMA模型和时空Kriging插值联合模拟大跨结构的脉动风速时程[J].振动与冲击,2014,33(3):29-34

[13] Zhang H, Yang JM, Xiao LF,. Large-eddy simulation of the flow past both finite and infinite circular cylinders at=3900[J]. Journal of Hydrodynamics Series B, 2015,27(2):195-203

[14] 汪权,王建国,张鸣祥.高层建筑结构随机风速场的数值模拟及风振控制[J].建筑科学与工程学报,2011,28(2):32-37

[15] 赖马树金.大跨度悬索桥分离式双箱梁涡激振动研究[D].哈尔滨:哈尔滨工业大学,2013

[16] Zeng C, Li CW. A hybrid RANS-LES model for combining flows in open-channel T-junctions[C]. Shanghai: 9thInternational Conference on Hydrodynamics, 2010:154-159

[17] 孙芳锦,梁爽.考虑流固耦合作用的大跨度桥梁风振响应研究[J].公路交通科技,2015,32(7):84-91

[18] 季小勇,陈文礼.斜拉桥斜拉索二维涡激振动的数值模拟研究[J].公路工程,2012,37(3):6-11

[19] Zheng Z, Lei J, Wu X. Numerical Simulation of the Negative Magnus Effect of a Two-Dimensional Spinning Circular Cylinder[J]. Flow Turbulence & Combustion, 2016,98(1):1-22

[20] 徐闯,孙延国.重庆永川长江大桥涡激振动性能风洞试验研究[J].世界桥梁,2015(5):25-29

[21] 李薇.斜拉索风雨激振试验与数值模拟研究[D].西安:长安大学,2012

[22] 刘克同,汤爱平,曹鹏.桥梁气动导数的格子Boltzmann大涡模拟仿真[J].工程力学,2015,32(5):111-119

[23] 乔永亮,桂洪斌,刘祥鑫.三维圆柱绕流数值模拟湍流方法的选择[J].水利水运工程学报,2016(3):119-125

[24] 谭艳,许栋,白玉川,等.关于圆柱绕流水动力特性的大涡数值模拟研究[J].港工技术,2013(3):5-8

[25] 林志兴,杨立波,李文勃,等.斜拉索顺桥向风阻系数的试验研究及数值模拟[C]//2004全国结构风工程实验技术研 讨会论文集.长沙:中国土木工程学会,2004:51-56

[26] 杨伟,顾明.高层建筑三维定常风场数值模拟[J].同济大学学报:自然科学版,2003,31(6):647-651

Numerical Simulation of Wind Field on Large Span Bridge Cable Based on Pulsating Winds

SUN Ting-ting1,2, YANG Ji-xin2*, LI Jian-hua2, ZHOU Xing-yu3

1.246003,2.430063,3.430010,

MATLAB is used to simulate the instantaneous wind speed time history of pulsating wind and import it into computational fluid dynamics software by using Kaimal fluctuating wind speed spectrum, the RNG-model is chosen to simulate the wind field of long span bridges under the action of high Reynolds number. The flow field model was established by using structural grid, and the changes of flow field and wind pressure coefficient were analyzed at the heights of 30 m and 80 m respectively. The velocity vector is analyzed at a height of 50 m, the flow rate and trace change were analyzed at the height of 80 m. The center and the end of the cable were selected to analyze the pressure distribution around the cable, and the numerical simulation was carried out. The calculation results showed that, the characteristics of fluid distribution in the flow field are in accordance with the conclusion of similar research, and a reasonable accurate numerical value can be obtained. Vortex was formed on the upper and lower sides of the cable, and the eddy current moved along the central axis and took a symmetrical shape, the pressure distribution of the cable surface was consistent with that of the wind tunnel in the literature, and in line with the engineering reality.

Kaimal pulsating wind speed spectrum; three-dimensional wind field; RNG-model; numerical simulation

O368;U441+.2; TU312+.1

A

1000-2324(2018)05-0809-06

10.3969/j.issn.1000-2324.2018.05.017

2017-06-09

2017-06-25

国家高技术研究发展计划(八六三)项目(2007AA11Z107);安徽省高校自然科学研究项目(KJ2016A448);安徽省高校自然科学研究项目(KJ2017A787)

孙亭亭(1986-),男,博士研究生,主要从事桥梁结构工程方面的研究工作. E-mail:kunpengting@163.com

通讯作者:Author for correspondence.E-mail:whutvses@163.com