大跨索承桥梁流线型钢箱梁抖振响应流固耦合数值模拟

徐梓栋 王浩 刘震卿

摘要 大跨索承桥梁抖振性能评估是桥梁抗风安全设计的重要环节。但当前抖振分析理论在揭示抖振基本物理成因时的“描述性”强于“解释性”，且风洞试验难以重现抖振过程中的流固耦合细节。对此，以苏通大桥标准主梁节段为研究对象，在Fluent软件中开展了模型抖振响应数值模拟，并将抖振响应数值解与理论解进行对比。结果表明，对高湍流度风场而言，模型风振以湍流引起的强迫振动为主，考虑流固耦合效应的抖振响应理论解与数值解的RMS较为接近，模型上、下方涡核发展交替进行，特征湍流对模型抖振响应的影响较小。

关键词 大跨桥梁; 抖振响应; 计算流体力学; 流固耦合; 大涡模拟

引 言

大跨索承桥梁抖振是由空气中湍流引起的限幅强迫振动，对大跨桥梁抖振性能进行评估是桥梁抗风安全设计的重要环节［1?2］。目前，工程中常用的抖振分析理论均引入了大量的气动参数对抖振荷载进行建模，且各类气动参数需由风洞试验或数值模拟等方法确定，因而既有抖振分析理论被认为具有半试验性质与统计性质［3?4］。此外，既有抖振分析理论以最大程度刻画桥梁抖振现象并服务于工程实际为目的，在工程领域得到广泛应用，但其难以揭示抖振现象的基本物理成因，且实际中重要的流固耦合（Fluid?Structure Interaction， FSI）细节容易被忽视［3，5?6］。因而，有必要发展更加精细化的大跨桥梁抖振响应分析方法，尽量避免使用试验气动参数与经验模型。

为充分把握空气与桥梁结构的FSI细节，实现大跨桥梁气动响应的精细化模拟，研究者们通过结合计算流体动力学（Computational Fluid Dynamics， CFD）与结构动力学方法，对桥梁?空气耦合系统动力方程进行求解［6］。例如，Nariman［7］采用双向FSI对大跨索承桥梁的涡激振及其“锁定”现象进行了分析；詹昊等［8］基于用户自定义函数（User Defined Function， UDF）对Fluent软件进行二次开发，实现了大跨桥梁主梁颤振的FSI数值仿真，并据此开展了桥梁颤振主动控制研究；董国朝等［9］依托ANSYS CFX软件及User Fortran工具开发了流固耦合分析程序，并据此开展了斜拉桥拉索风致驰振FSI分析；Kavrakov等［10］采用二维模型开展了大跨桥梁主梁节段抖振响应FSI分析，并将数值解与理论解进行了对比。总体上，现有大跨索承桥梁数值模拟研究多针对如颤振、涡振等自激振动，而受大气边界层湍流难以准确模拟、数值计算开销巨大、气固耦合机理难以准确把握等因素制约［10?11］，对由湍流引发的大跨桥梁抖振响应数值模拟研究并不多见。

大跨索承桥梁风致抖振是由大气边界层湍流引起的强迫振动，因而合理高效地模拟大气边界层湍流风环境是开展风致抖振分析的重要前提。大涡模拟（Large Eddy Simulation， LES）是CFD中重要的湍流处理方式，其对大涡部分进行直接求解，而对亚尺度涡旋进行模化。利用LES开展大气边界层湍流风环境模拟的关键在于生成满足特定条件的入口湍流［11?12］。当前，LES入口湍流生成技术包括层?湍流过渡技术、循环域技术、驱动域技术与合成湍流技术四类［13］。其中，合成湍流技术又包括谱方法、涡方法与滤波法三类，因计算开销小且可实现入口边界平均风剖面与湍流特性的精准控制，该技术被广泛应用于湍流大气边界层模拟［14］。综合考虑各类LES湍流入口生成技术在计算效率、湍流结构可控性、大气边界层模拟适用性等方面的优缺点，拟采用合成湍流技术中的窄带合成（Narrowband Synthetic Random Flow Generation， NSRFG）法［14］生成满足给定风特性的湍流风场，并据此实现主梁节段抖振响应的数值模拟。

本文依托Fluent软件平台，基于Fluent中的用户自定义函数（User Defined Function， UDF）功能植入NSRFG法，从而生成满足给定风特性的湍流风场。此外，通过UDF编制苏通大桥主梁标准节段的抖振响应FSI模拟程序，开展该桥三维主梁节段抖振响应数值模拟。为确保模拟结果合理可靠，以大桥钢箱梁气动参数的风洞试验结果为参考，对数值模型进行验证，并将抖振响应数值解与半经验模型获得的理论解进行对比分析。

1 苏通大桥标准主梁节段数值模型

1.1 几何建模与网格划分

取模型几何缩尺比1∶50，并参考文献［15］构造三维计算域，其中计算域z方向厚度16 m为苏通大桥标准主梁节段长度，计算域x?y平面尺寸及边界条件设置如图1所示。

图1中，钢箱梁重心位置坐标为（x， y=0 mm， z=160 mm）。为便于抖振响应数值模拟，计算域网格划分采用多域网格划分方案以降低对模拟时间步长的要求。各网格区域由interior交界面分隔，各区域网格在交界面上共节点。在充分利用计算资源的前提下，计算域采用几何适应性较强的非结构网格划分，边界层采用三棱柱网格划分，第一层边界层网格厚度为0.001 m，壁面y+平均值45，边界层中的黏性亚层与缓冲层则通过壁面函数进行模化，计算域网格总数超460万（4611692个），动网格区域及边界层网格划分情况如图2所示。

图2中网格划分在Hypermesh软件中完成，计算域网格最大等效扭曲度（Equivalent Skew）小于0.7，满足Fluent计算要求。

1.2 数值模型气动参数分析与验证

节段模型的风致抖振主要由抖振荷载引起，抖振理论表明抖振荷载与湍流风场及断面三分力系数密切相关［10］。其中，湍流風场可通过LES实现，数值模型的三分力系数则需要依据风洞试验结果进行验证。

数值模型三分力系数分析时，关闭动网格设置，采用雷诺平均N?S方程（Reynolds?Averaged Navier?Stokes， RANS）进行瞬态（Unsteady）计算，湍流模型取Realizable k ?ε模型，设置变量与连续性方程残差小于10–5为收敛标准，时间步长取0.002 s，时间步总数13000步，离散与插值格式设置如表1所示。以1.1节中数值模型为分析对象，分别对模型-3°到+3°攻角下（步长1°）的三分力系数进行计算，并与风洞试验结果进行对比，如图3所示。

由图3可知，升力系数和扭矩系数与风攻角正相关，而阻力系数在风攻角为负时，随着风攻角的增大而降低，反之增加。受边界层网格划分和风特性参数设置等因素影响，模拟所得阻力系数的绝对误差相对较大，但各工况下三分力系数模拟值与风洞试验值的相对误差均不超过20%，且变化趋势大体相同，表明该模型可进一步用于后续抖振响应数值模拟。

2 抖振响应数值模拟流程与参数设置

2.1 模拟流程

对于二自由度主梁断面，考虑竖向与扭转自由度的动力学方程可写为［10］：

式中 h与α分别表示竖向与扭转位移；x˙表示变量x对时间求导；mh为单位长度主梁质量；mα为单位长度主梁转动惯量；ci=2miξiωi （i=h，α）为阻尼系数，其中ωh=1.1567为实测主梁一阶正对称竖弯频率，ωα=25.244为实测主梁一阶正对称扭转频率；ξi（i=h，α）表示阻尼比，各方向均取0.4%；ki=miω2i （i=h，α）为刚度系数；L与M分别表示作用于主梁的气动升力与扭矩，其线性非定常模型表达如下［10］：

式中 χiu与χiw （i=1，2，3，4）为气动导纳函数，本文偏于保守的取1；K=ωB/U表示折减频率；H*i与A*i （i=1，2，3，4）表示颤振导数；u表示顺风向脉动风速；w表示竖向脉动风速；U表示平均风速；ρ表示空气密度；B表示梁宽。苏通大桥主梁断面抖振响应数值模拟在Fluent中直接完成，流体方程与结构动力学方程耦合求解流程如图4所示。

图4中，首先通过NSRFG法生成LES入口湍流，再编制UDF提取每一时间步上作用于主梁断面的气动升力与扭矩，并使用4阶Runge?Kutta法求解式（1）中主梁断面动力方程，从而确定主梁断面的运动状态。最后，根据Fluent中自带的Remeshing与Smoothing方法实现网格更新。上述耦合过程由CFD求解器与UDF_BUFFETING间的数据交换体现，属于双向弱耦合。

2.2 模拟参数设置

为使数值模型满足相似准则，设模型几何缩尺比n=50，则模型中其他参数的相似关系与取值如表2所示［1，16］。LES入口湍流风特性定义如表3所示［14］。

开展抖振数值模拟时，取表3中α=0（常数风剖面），湍流强度则参考苏通大桥桥址区长期风特性实测结果取值（10 min基本时距湍流强度平均值0.2051）［17］，风场湍流度较高。模拟在配置为96核Intel（R） Xeon（R） Platinum 8163 CPU @ 2.50 GHz 的工作站上完成，采用32核并行计算。首先，关闭动网格，采用LES进行瞬态分析使流场发展稳定。再打开动网格，进行主梁断面抖振响应FSI模拟，LES时间步长为0.002 s，时间步总数28000。设置变量与连续性方程残差小于10-5为收敛标准，模拟总耗时超1800 h（2个半月）。LES离散与插值格式设置如表4所示。

3 模拟结果分析与验证

为开展模拟结果的量化分析与验证，在主梁节段附近设置监测点BP（-0.5 m， 0 m， 0.16 m），测点处顺风向与竖向脉动风速时程如图5所示，相对应的湍流谱密度如图6所示。

由图5和6可知，随着湍流在流域内不断发展，模型位置处功率谱密度监测值仅在低频区与目标值吻合良好，而在高频区迅速衰减。该现象一方面与数值算法产生的湍动能耗散有关；另一方面由于LES采用滤波函数过滤了风速时程中的高频成分，从而使LES得到的湍流功率谱在高频部分偏低。尽管如此，主梁节段模型竖弯与扭转频率仍位于LES解析尺度内，故生成的湍流场满足模型抖振模拟要求。将BP测点脉动风速时程代入式（2），并采用4阶Runge?Kutta法求解式（1），获得考虑与不考虑气动自激力时主梁节段抖振响应的理论解，其中模型气动参数根据文献［18］取值。将其与数值解进行对比，以验证模拟结果的合理性，如图7所示。相对应的位移响应功率谱密度分析结果如图8所示。

由图7和8可知，当模拟风场湍流度较高时，主梁节段模型风振响应兼具强迫振动（随机性）与自激振动（周期性）特征，但以强迫振动为主。就时程结果而言，由于抖振力模型存在假设简化，且数值模型气动参数与风洞试验结果存在差异，故数值解与理论解并不完全吻合。当不考虑自激力时（即不考虑流固耦合效应），竖向位移理论解明显偏大，而扭转位移理论解则无明显变化。此外，抖振响应数值解与理论解所包含的频率成分相近，且主要频率成分与主梁节段竖弯及扭转频率相对应。对随机振动而言，统计量较时程结果更具代表性，故分别计算了抖振响应数值解与理论解的均方根值（Root Mean Square， RMS），结果如表5所示。

由表5可知，不考虑自激力时的竖向位移RMS相对误差高达60.7%，而考虑自激力时，竖向位移RMS相对误差仅1.7%。扭转位移RMS较小，是否考虑自激力对扭转位移RMS的影响也较小。实际工程中不考虑自激力获得的抖振响应总体上偏于安全。

为进一步探究抖振发生时的绕流场细节，以计算域z=160 mm平面（Plane 1）为监测面，任意选取计算稳定时主梁断面附近瞬时速度场进行展示，并基于Q准则对主梁断面附近的涡核进行识别，速度场与瞬時涡核分布如图9所示。其中，位移和加速度的数值模拟结果（缩尺后结果）也一并给出，并将对应时刻的位移结果标注于图中。

由图9可知，模型上方涡旋形成于迎风侧桥面端点处，而下方涡旋则形成于迎风侧风嘴端点处。模型上、下方涡核经充分发展后，均在模型中心线附近达到最大尺寸。特别地，当模型下方涡核破碎消散后，又将在底板与背风侧斜腹板交点处重新形成，并向尾涡区发展。此外，模型上、下方涡核发展交替进行。当上方涡核不断发展时，下方涡核则逐渐破碎消散，反之亦然。而模型附近涡核演化（特征湍流引起）与模型抖振响应间未见明显关系，表明特征湍流对模型抖振响应的影响较小。

4 结 论

本文基于NSRFG法实现了LES入口湍流生成，通过编制UDF实现了大跨索承桥梁流线型钢箱梁抖振响应数值模拟，并将抖振响应数值解与理论解进行了对比分析，得到主要结论如下：

（1）三分力系数的CFD计算值与风洞试验值相对误差不超过20%，且变化趋势大体相同，考虑流固耦合效应的抖振响应理论解与数值解的RMS较为接近。

（2）是否考虑自激力对扭转位移RMS的影响都较小，但不考虑自激力时的竖向抖振位移的RMS则明显偏大，实际工程中不考虑流固耦合效应获得的抖振响应总体上偏于安全。

（3）当模拟风场湍流度较高时，主梁节段模型风振响应兼具强迫振动（随机性）与自激振动（周期性）特征，但以强迫振动为主。模型上、下方涡旋分别形成于迎风侧桥面端点和风嘴端点处，两侧涡核经过充分发展后，均在模型中心线附近达到最大尺寸，且下方涡核破碎消散后，又将在底板与背风侧斜腹板交点处重新形成，并向尾涡区发展。

（4）模型上、下方涡核发展交替进行。当上方涡核不断发展时，下方涡核则逐渐破碎消散，反之亦然。模型附近涡核演化（特征湍流引起）与模型抖振响应间未见明显关系，特征湍流对模型抖振响应的影响较小。

与理论分析及风洞试验相比，抖振响应数值模拟充分捕捉了模型的绕流场特性，并清晰展示了抖振发生时模型绕流场的演化规律，研究成果有望进一步用于主梁断面抖振性能优化及桥面風环境评估等。然而，受计算资源限制，数值模拟时的LES解析尺度尚无法完全涵盖实际风场中的湍流尺度，模拟所得数值风场与实际风场仍存在一定差异。因而，在计算资源允许时，仍可开展更加精细化的数值模拟，以进一步验证本文结果。

参考文献

1陈政清. 桥梁风工程［M］. 北京： 人民交通出版社， 2005.

2Ge Y， Xia J， Zhao L， et al. Full aeroelastic model testing for examining wind-induced vibration of a 5，000 m spanned suspension bridge［J］. Frontiers in Built Environment， 2018， 4： 20.

3Simiu E， Scanlan R H. Wind Effects on Structures［M］. New York： John Wiley and Sons， 1996.

4葛耀君. 大跨度桥梁抗风的技术挑战与精细化研究［J］. 工程力学， 2011， 28（增刊 Ⅱ）： 11-23.

Ge Y J. Technical challenges and refinement research on wind resistance of long-span bridges［J］. Engineering Mechanics， 2011， 28（Sup Ⅱ）： 11-23.

5Ahsan K. Advanced Structural Wind Engineering［M］. New York： Springer， 2013.

6Tamura T. Reliability on CFD estimation for wind-structure interaction problems［J］. Journal of Wind Engineering and Industrial Aerodynamics， 1999， 81（1-3）： 117-143.

7Nariman N A. Influence of fluid-structure interaction on vortex induced vibration and lock-in phenomena in long span bridges［J］. Frontiers of Structural and Civil Engineering， 2016， 10（4）： 363-384.

8詹昊， 廖海黎. 桥梁主梁端部翼板颤振主动控制流固耦合计算［J］. 振动工程学报， 2018， 31（2）： 276-282.

Zhan H， Liao H L. Flutter active control studies of bridge with flaps attached to deck edge by FSI calculation［J］. Journal of Vibration Engineering， 2018， 31（2）： 276-282.

9董国朝， 陈政清， 罗建辉， 等. 安装亮化灯具导致的斜拉桥拉索风致驰振流固耦合分析［J］. 中国公路学报， 2012， 25（1）： 67-75.

Dong G C， Chen Z Q， Luo J H， et al. Fluid-structure interaction analysis of wind-induced galloping of cables with lamps of cable-stayed bridge［J］. China Journal of Highway and Transport， 2012， 25（1）： 67-75.

10Kavrakov I， Morgenthal G. A synergistic study of a CFD and semi-analytical models for aeroelastic analysis of bridges in turbulent wind conditions［J］. Journal of Fluids and Structures， 2018， 82： 59-85.

11周桐， 杨庆山， 闫渤文， 等. 大气边界层大涡模拟入口湍流生成方法综述［J］. 工程力学， 2020， 37（5）： 15-25.

Zhou T， Yang Q S， Yan B W， et al. Review of inflow turbulence generation methods with large eddy simulation for atmospheric boundary layer［J］. Engineering Mechanics， 2020， 37（5）： 15-25.

12刘震卿， 张冲， 吴晓波， 等. 考虑粗糙地表的近地层风场仿真模拟研究及误差分析［J］. 空气动力学学报， 2020， 38（4）： 715-722.

Liu Z Q， Zhang C， Wu X B， et al. Simulation and error analysis of wind field considering surface roughness［J］. ACTA Aerodynamic SINICA， 2020， 38（4）： 715-722.

13Yan B W， Li Q S. Inflow turbulence generation methods with large eddy simulation for wind effects on tall buildings［J］. Computers and Fluids， 2015， 116： 158-175.

14Yu Y， Yang Y， Xie Z. A new inflow turbulence generator for large eddy simulation evaluation of wind effects on a standard high-rise building［J］. Building and Environment， 2018， 138： 300-313.

15Grinderslev C， Lubek M， Zhang Z. Nonlinear fluid-structure interaction of bridge deck： CFD analysis and semi-analytical modeling［J］. Wind and Structures， 2018， 27（6）： 381-397.

16許福友， 马如进， 陈艾荣， 等. 苏通大桥全桥气弹模型设计与模态调试［J］. 工程力学， 2009， 26（12）： 150-154.

Xu F Y， Ma R J， Chen A R， et al. Full aeroelastic model design and modal test for Sutong Bridge［J］. Engineering Mechanics， 2009，26 （12）： 150-154.

17Wang H， Li A， Niu J， et al. Long-term monitoring of wind characteristics at Sutong Bridge site［J］. Journal of Wind Engineering and Industrial Aerodynamics， 2013， 115： 39-47.

18Wang H， Hu R， Xie J， et al. Comparative study on buffeting performance of Sutong Bridge based on design and measured spectrum［J］. Journal of Bridge Engineering， 2013， 18（7）： 587-600.

Buffeting numerical simulation of streamlined steel box girder of long-span cable-supported bridge using fluid-structural interaction

XU Zi-dong 1，2 ?WANG Hao 2 ?LIU Zhen-qing 3

1. Key Laboratory of Concrete and Prestressed Concrete Structures of Ministry of Education， Southeast University， Nanjing 211189，China；

2. School of Civil Engineering， Inner Mongolia University of Science and Technology， Baotou 014010，China；

3. School of Civil and Hydraulic Engineering， Huazhong University of Science and Technology， Wuhan 430074，China

Abstract Evaluation on the buffeting performance of the long-span cable-supported bridge is the significant step for the bridge wind-resistant safety design. However， current buffeting theories are more descriptive， but not explanatory in the sense of revealing basic buffeting physical causes. In addition， the details of the fluid-structure interaction （FSI） during the buffeting are hard to be reproduced. In this work， the Sutong Bridge is selected as the research object， the buffeting responses of its numerical model are simulated in the software Fluent. And the numerical buffeting responses are compared with the theoretical ones. Results show that the forced vibration is the major component of the wind-induced vibration in the wind field with high turbulence intensity. The theoretical and numerical RMS of the buffeting responses are close when considering FSI effect. Vortex core development conducts alternatively between the upper and bottom of the deck. Whats more， signature turbulence has little effect on the buffeting responses of the model.

Keywords long-span bridge; buffeting response; computational fluid dynamics; fluid-structural interaction; large eddy simulation