恢复交通首日的虎门大桥现场监控截图作者 | 陈东阳博士  仿真秀专栏作者首发 | 仿真秀公众号(ID:fangzhenxiu2018)导读：据广东省交通集团发布消息，5 月 15 日 9 时，虎门大桥正式恢复交通。虎门大桥异常抖动的原因基本查明：沿桥梁边护栏连续设置的水马改变了钢箱梁气动外形，在特定的风况条件下，诱发悬索桥发生竖向涡激共振；持续较长时间、较大振幅的涡振导致桥梁结构阻尼下降，即使水马拆除，特定风况下空载桥梁又发生了涡振。自5月5日虎门大桥异常抖动以来，引发了社会和专业人士的强烈关注，虎门大桥异常抖动有关的撰文数量达到9722篇，传播范围极广。其中西北工业大学张伟伟博士团队撰稿《冯卡门也犯错？大桥振动和卡门涡街到底什么关系？》和仿真秀APP专栏作者徐伟博士原创的《虎门大桥异常抖动背后：流固耦合仿真方法解析》受到了专业人士的高度好评。为了帮助用户更好的理解流固耦合与工程应用，推荐您关注仿真秀APP专栏作者，墨尔本皇家理工大学航空工程学院访问学者——陈东阳博士与我们分享弹性支撑的二维柱体结构涡激振动仿真和基于CFD/FEM安装螺旋列板的立管涡激振动研究。感兴趣的朋友，还可报名陈东阳博士5月21日19时30分在仿真秀平台的公开直播，详情见后文。一、弹性支撑的二维柱体结构涡激振动仿真1、柱体结构涡激振动定义对于海洋工程、风工程上普遍采用的圆柱形断面结构物，流体绕过柱体时会产生交替发放的泻涡，这种交替发放的泻涡又会在柱体上生成顺流向及横流向周期性变化的脉动压力。如果此时柱体是弹性支撑的，或者柔性管体允许发生弹性变形，那么脉动流体力将引发柱体(管体)的周期性振动，这种规律性的柱状体振动反过来又会改变其尾流的泻涡发放形态。这种流体-结构物相互作用的问题被称作"涡激振动"(Vortex-Induced Vibration :VIV)，它是典型的一种流固耦合振动现象。2、柱体结构涡激振动动力学模型及仿真方法简介进行二维弹性支撑柱体的数值模拟是研究海洋工程和风工程中柱体结构涡激振动现象和机理的重要手段。根据牛顿第二定律，2-DOF弹性支撑的柱体运动的控制方程可以写为：(1)

式中， m为圆柱体的质量， c为结构阻尼系数， k为结构刚度系数。

式(1)又可以写为：(2)式(2)中，柱体固有频率为：

阻尼比为：对于涡激振动的数值模拟大多以柱体的横向振动研究为主，同时考虑横向和来流向的耦合振动研究相对较少。其主要原因在于，对柱体结构涡激振动数值模拟，必须保证柱体周围网格质量非常好，才能有效预测涡激振动响应。一般使用网格质量高的结构化网格，但是当柱体结构发生较大振动位移时，周围流场网格会发生畸变，甚至产生负网格，导致计算失败。如果再考虑流向耦合振动，计算难度将非常大，计算成功率也将明显降低。如果采用非结构化网格，且采用网格重构技术，可以吸收柱体较大的振动位移，但是非结构化网格质量相比结构化网格质量较差，且采用非结构化网格必定需要大大增加网格量，从而大大增加了计算时间。此外，如果对于涡激振动抑制装置设计研究，也就是在柱体结构表面形状变复杂的情况下，网格划分难度和计算难度也将大大增加。因此，寻求一种既可以保证网格质量，又能不大幅度增加网格数量，且可以避免网格畸变或者负网格问题的方法十分重要。基于CFD商业软件FLUENT和结构动力学原理，通过用户自定义函数(UDF)及嵌套网格技术，建立了2-DOF弹性支撑柱体结构VIV数值模型。非定常不可压缩流体RANS方程为：(4)(5)式(5)中，式中，为不可压缩流体的密度；表示 方向上的瞬时速度分量， 为 方向上速度脉动量，为速度的时间平均值；分别表示笛卡尔坐标系、时间、压力、运动粘度；为湍流黏度，下标“t”表示湍流；为湍动能；是“Kronecker delta”符号，就是当时，，当时，。湍流模型选用SST k–ω湍流模型。通过计算流场，可以得到二维柱体表面的压力分布，进而可以得到作用在二维柱体上的升力和阻力系数：

(6)

(7)结合方程(2)，2-DOF弹性支撑的柱体运动的控制方程可以写为：

(8)

(a) 两自由度弹性支撑刚性柱体(b) 二维两自由度弹性支撑刚性柱体VIV模型

图1 2-DOF弹性支撑圆柱体VIV模型示意图

图2 2-DOF弹性支撑圆柱体流场计算网格

两自由度弹性支撑刚性柱体在流体作用下的结构示意图如图1(a)所示，二维2-DOF振动柱体VIV模型示意图如图1(b)所示。一般柱体流场的尾迹区域需要大于等于22.5D(D为柱体直径)，整体局域高度一般需要大于等于20D，柱体振动才不受流体区域边界的影响。

因此，综合考虑计算条件的情况下，流场域的尺寸大小如图1(b)中标注所示，尾迹区域30D，柱体前端和上下距离柱体都是10D。包围这柱体的组分网格外边界直径大小为3D。流场入口边界条件为速度入口，出口为压力出口，上下壁面为滑移壁面，柱体表面即动边界为无滑移壁面。

流场随着柱体边界的改变而改变，通过动网格技术来实现流场中柱体边界的运动。嵌套网格技术是最新的动网格技术，主要适用于刚性边界运动问题。如图2所示，流场域网格划分采用的是嵌套网格。

如图2(a)所示，背景网格和嵌套网格都使用结构化网格，靠近柱体表面部分为边界层网格()，较好的保证了网格质量。采用嵌套网格技术，可以无需担心网格畸变以及负网格导致求解失败等问题。同时，不会较多的增加计算量。嵌套网格即多重网格相互重叠组合成的一组网格。有可能存在两套或者两套以上的网格相互重叠。

嵌套网格求解的大致思路为：首先划分包裹柱体的组分网格(组分网格数量为5262)，和外流场的背景网格(背景网格数量为15731)，求解器识别嵌套网格边界，对被组分网格遮蔽的背景网格部分进行“挖洞”，然后对嵌套区域边界单元进行插值，将背景区域的边界单元变量信息插值到嵌套区域的边界单元(如图2(b)所示)，最后进行流场计算。

整个流场的计算网格如图2(c)所示。对于流场的数值计算，时间项采用全隐式积分方法，对流项则采用二阶迎风离散格式。控制方程中速度分量与压力的耦合则采用COUPLED算法进行处理。初始条件为：

图3 2-DOF弹性支撑圆柱体VIV计算流程图

流场域求解基于CFD商业软件FLUENT，根据边界条件获得流场和二维柱体表面的压力、速度等信息。提取作用在柱体表面的力，然后代入柱体的结构运动方程，通过求解二维柱体的运动方程，得到当前时间步长下的柱体运动的位移和速度。同时利用得到的柱体位移和瞬时速度，更新流场网格，然后进行下一个时间步的计算。这个双向流固耦合仿真过程是通过FLUENT软件的用户自定义函数(UDF)来实现的。

UDF中可以使用标准C语言的库函数，也可使用FLUENT中预定义的宏。通过预定义宏可以获得FLUENT计算过程中的流场数据。FLUENT中用户自定义函数是通过DEFINE宏来实现的。基于CFD的2-DOF弹性支撑柱体VIV数值求解的计算流程图如图3所示。图中的虚线框内为通过C语言编制的UDF程序来实现。

3、部分仿真结果展示

图4 基于CFD模型的2-DOF柱体VIV计算结果

柱体涡激振动的最大幅值及频率比随不同约化速度的变化如图4所示，与实验数据对比，误差较小，验证了本文计算方法的正确性。

从图4(a)可以看出，数值仿真出3种响应分支，在在3~4之间的时候，原始分支向上端分支转变；当=5~6之间的时候，出现下端分支。在上端分支中振幅达到最大值0.98，而在下端分支中振幅最大值0.642；从=11的时候开始，圆柱体的响应位移又回落到一个很小的数值。从图4(b)可以看出，在频率“锁定”区间=4~10内，柱体的实际振动频率与固定柱体的泄涡频率 分离，不再符合st与Re数关系图。同时，柱体的实际涡泻频率与柱体固有频率比值稳定在1.15附近，而在解锁区域，柱体的实际振动频率 与固定柱体的涡脱频率相同，这与前人的实验结果大致相同。

图5 时，不同时刻的涡量云图(周期T=3.34s)

图5给出了=5时弹性支撑柱体的75s~78.5s的涡量云图，包含了一个周期的运动，从图中可以看出，=5时的涡脱模式为P+S模式(即一个涡脱周期内有一个单个涡+一对涡形成)。Govardhan 和 Williamson 的实验研究表明一般在柱体振幅较大时候涡脱模式为P+S或者2P(即在一个涡脱周期内有2对尾涡形成)，在振幅较小的时候涡脱模式为2S(即在一个涡脱周期内有2个单独的尾涡形成)。图5中的黑虚线为柱体的原始位置，红点为柱体当前时刻的中心位置，从图中可以看出，柱体振动游走的轨迹是一个“8”字形。

二、基于CFD/FEM安装螺旋列板的立管涡激振动研究

1、 基于CFD/FEM双向耦合的立管涡激振动模型验证

本节以Bijan Sanaati的实验模型为研究对象，立管的两端边界为铰-铰约束，同时立管顶端顶张力大小为T，如图6所示。文献所用实验模型为PVC管，详细结构参数、周围流体环境参数、内部流体介质参数根据文献换算出的参数如表1所示。

图6 柔性立管模型

建立PVC立管的外流场，其中PVC管的直径为D。入口处到立管的距离为10D，两边到立管的距离也为10D，出口处到立管的距离为30D。进口处采用速度进口边界条件，出口处采用压力出口边界条件，其他面都采用对称边界条件。流体域模型如图6(b)所示。流体与壁面间相互作用，很多因变量具有较大的梯度，而且粘度对传输过程有很大的影响。采用SST湍流模型要求近壁处有y+<1。计算网格如图7所示，在沿着管长方向分布了100个节点，网格数量为2776653。

图7 柔性立管流场域计算网格

阻力系数为，升力系数为。分别为升力、阻力，L是立管的长度。斯特鲁哈数，为漩涡脱落的频率，通过对升力系数时域历程作FFT变换得到。约化速度，其中是立管的湿模态频率。由于立管周围的尾涡具有三维特性、随机特性，采用CFD/CSD全三维数值模拟可以克服二维模拟或者切片法带来的误差。

本节基于CFD/CSD双向流固耦合方法，分别计算了顶张力为60N和260N时铰-铰约束的PVC立管动态响应。CFD/CSD双向耦合流程图如图8所示。

图8 柔性立管双向耦合流程图

从图9中可以看出，立管的质量比为1， 之间为本节PVC立管的频率锁定区间，该区间内幅值相对周围区间有显著的增大。本文计算了八组结果，与文献实验数据对比，误差较小，说明基于CFD/CSD双向流固耦合仿真精度较高，也验证了CFD/CSD双向耦合计算模型的准确性。

图10为T=60N, U=0.15m/s时的部分仿真结果。图10(a)为立管横向和纵向振幅值均方根沿着立管长度方向的分布，从图中可以看出，横向振幅大于来流向振幅，立管最大幅值发生在立管中部，幅值无量纲量为0.42，与处实验数据很接近。

图10(b)为立管横向振动幅值无量纲量的历时云图，图中红色表示正值，蓝色表示负值，从图中可以看出，该流速情况下，立管主要以第1阶模态振动为主。图10(c)为立管中间点处的运动轨迹，从图中可以看出，运动轨迹比较杂乱，横向和来流向的振动位移响应都比较大。

图10(d)、(e)分别为沿立管长度方向横向和来流向的振动响应频率谱分析结果，横向振动频率几乎都是1.42Hz，与文献的实验数据1.52Hz接近，说明仿真结果较好。来流向振动频率几乎都为2.7Hz。来流向振动响应频率大概是横向振动响应频率的两倍。此处采用CFD/FEM双向流固耦合的计算量非常大，仅仅该八组数据采用16核32线程工作站并行计算共耗时近一个月时间。

图9 立管Z方向幅值均方根随着约化速度变化图

图10 仿真结果 (T=60N, U=0.15m/s)

2、安装有螺旋列板的立管涡激振动响应研究

本节研究不同结构参数的螺旋列板对立管涡激振动特性的影响规律。目前采用全三维CFD/CSD计算立管的涡激振动响应，计算量太大。本节采用从细长立管上截取一小段作为研究对象，用短立管来研究安装和不安装螺旋列板两种情况下立管的涡激振动特性。前人对于螺旋列板的研究大量的工作大多基于实验和二维仿真。本节提供了一种可以考虑螺旋列板三维效应的方法，且计算量适中。

图11 带螺旋列板的立管结构示意图

表2 立管数值仿真工况

计算短立管采用如下参数, Re=200，短立管两端采用固定-固定约束。因此，立管的第1阶固有频率

在Re=200左右时，

其中是漩涡泻放频率，因此

此时，漩涡泻放频率与立管固有频率接近，发生频率“锁定”现象，立管振幅变大。在不改变立管其他尺寸的情况下，取短立管的长度为10D。通过降低短立管的弹性模量和增加短立管的密度直到第1阶湿模态频率等于0.1238Hz，且质量比等于10。采用这样的参数，就可以保证立管的频率锁定区间较小，计算更容易实现。带有螺旋列板的立管结构示意图如图11所示，螺旋列板的齿高为h ，螺距为p。所研究的螺旋列板的列数为1、2、3；螺距为 为5D、10D 、15D；齿高为0.05D、 0.1D、0.2D；螺旋列板覆盖率为100% 和1/3 (可安装在短立管的上部、中间、底部)。螺旋列板采用和立管相同的材料，表2给出了十一组立管数值仿真工况的计算条件。

图12 立管中点的运动轨迹

图13 频谱图

当流体绕过立管时，漩涡以固定的频率泻放。计算结果表明，螺旋列板对立管涡激振动有很好的抑制作用。图14(a)为裸管的不同约化速度情况下的立管中点运动轨迹。从图中可以看出，运动轨迹成“8”字形，锁定区间为，在此区间内振幅较大。图14(b)为不同约化速度情况下带有螺旋列板(螺旋列板覆盖率100%，3列， p=10D, h=0.1D)的立管中点运动轨迹。从图中可以看出，加装了螺旋列板后，横向振动的振幅相比裸管明显减小很多。对来流向的振动没有明显的影响。为了更好的理解带有和不带有螺旋列板的立管的动力学行为，对带有和不带有螺旋列板的立管动力学响应作频谱分析。图13为情况下立管带有和不带有螺旋列板的频谱分析图，从图中可以看出，带有螺旋列板的立管主响应频率都接近于0，次频率都小于裸管的主响应频率。因此可以得出结论，螺旋列板抑制涡激振动机理是改变了漩涡泻放的频率。

图14 立管中点的运动轨迹

图15 涡核图和涡量云图

图14给出了螺旋列板齿高、螺距、列数、不同位置的1/3覆盖率对立管横向和来流向的涡激振动的影响。所有计算结果都与裸管的结果作对比。当齿高h = 0.1D，不同列数对立管的涡激振动响应影响如图14(a)所示。计算结果表明增加螺旋列板的列数可以大大减小横向振动的振幅，但是增加了立管来流向的弯曲变形。2列和3列螺旋列板对横向振动振幅抑制效果接近。

从图14(b)可以看出h=0.1D时不同螺距对立管涡激振动的影响，螺距越小，横向振动的振幅抑制效果越好，但同时增加了立管来流向的弯曲变形。图14(c)为3列，h=0.05D情况下不同齿高对立管涡激振动的影响。从图中可以看出，齿高h=0.1D和 h=0.2D对立管横向振动的振幅抑制作用差不多，但是h= 0.2D时大大增加了立管来流向的弯曲变形。

综合考虑，h=0.1D时，既可以较好的抑制立管的横向振动振幅，又不会引起立管来流向太大的弯曲变形。如图14(d)所示，将 1/3覆盖率的螺旋列板分别置于立管的上部、中间和底部，计算结果表明，将1/3覆盖率的螺旋列板置于立管的中部对横向振动振幅的抑制效果相比置于立管额上部和底部稍微好一点，置于立管上部，立管的来流向振幅也相对减小，综合考虑，将其置于立管上部相比置于立管中部和底部效果要好一点。

图15为立管带有和不带有螺旋列板情况下的涡核和涡量云图。从图中很明显可以看到，螺旋列板破坏了漩涡结构，这对漩涡泻放频率产生了较大影响，导致了漩涡泻放频率大大减小。

三、公开直播《流固耦合动力学仿真与工程应用》

受仿真秀平台邀请，笔者将与2020年5月21日19时30分在仿真秀直播间公开直播《流固耦合动力学仿真与工程应用》，感兴趣的朋友可以识别下方二维码报名，

注意：报名且参与直播的用户，本课程可以在仿真秀官网和APP开通长期反复免费观看。且提供课程资料和模型下载等服务。如有任何疑问请添加小助手微信咨询【fangzhenxiu999】

相关阅读

虎门大桥异常抖动背后：流固耦合仿真方法解析FLUENT流固耦合柱体结构涡激振动仿真案例解析工程设计中如何考虑流固耦合问题?

即将直播：流-热-固多场耦合计算关键技术及典型应用(2月20日)

作者:陈东阳博士，仿真秀专栏作者

声明：原创文章，本文首发仿真秀公众号(ID:fangzhenxiu2018),本文已经被授权，部分图片源自网络，如有不当请联系我们，欢迎分享，禁止私自转载，转载请联系我们

点击阅读原文，立即体验陈博士主讲的《Fluent流固耦合方法与技能20讲》