干货 | OpenFOAM涡激振动仿真详解

11月前浏览1633

涡激振动

当雷诺数超过40时，流体在圆柱体上流动会发生涡脱落现象，从而在圆柱体上施加周期性的作用力。如果圆柱体受到弹性支撑，则该力的横向分量就会导致圆柱体的横向振动，这一现象称为涡激振动现象。

弹性安装圆柱形结构的涡激振动现象在许多工程应用中经常碰到，如换热器管、海洋立管、海上平台结构、海底管道、海上风力涡轮机柱基以及输电线路等等。当涡脱落的频率接近系统的固有频率时，系统会发生大幅度横向振动，这就是锁定或同步现象。锁定现象发生时，流体在圆柱结构上施加大幅值周期性流体力，会导致这些柔性结构发生疲劳失效，造成一定的风险隐患。

近几十年，科学家们对涡激振动现象展开了大量的研究。在进行设计时，难免要与已有的文献研究结果进行比较。本文将从笔者自身经验出发，分享如何使用OpenFOAM模拟涡激振动，以及如何取得与文献相符合的结果的一些个人经验，供初学者参考。

OpenFOAM具体设置

OpenFOAM版本为esi的v2006版本，使用其动网格功能来实现涡激振动计算。本文所用到的算例在gitee上可供下载：

https://gitee.com/chen-jingle219/vortex_induced_vibration/blob/master/README.md

算例基本设置
雷诺数Re	150
圆柱直径D	0.0016m
运动粘度nu	1e^-6 m²/s（水）
来流速度U	0.09375 m/s
圆柱质量比m*	3
网格量	45k

流场网格

使用OpenFOAM自带的动网格技术（dynamicMesh），用到的网格是在ICEM画的结构网格。如下图所示，在圆柱附近设置了一个很大范围的圈，包裹着圆柱。正常的圆柱绕流算例不会那么复杂，但为了应对涡激振动在特定约化速度Ur下幅值较大的情况，特地留了充足的网格变形范围。如果处理不当，很容易造成局部网格变得极端小，从而计算发散。这里，圆柱直径D=0.0016m（沿用一些文献的尺寸），网格范围为30D X 50D。

0：边界条件以及初始条件

OpenFOAM边界条件以及初始条件在case/0文件夹中面统一设置，可以在我的gitee中看到。值得一提的是，要使用动网格，需要加pointDisplacement 文件。然后在U文件中，圆柱表面的边界类型由 fixedValue 改为 movingWallVelocity：

其它边界条件以及初始条件与圆柱绕流算例基本相同。

constant: 动网格设置

constant里面除了接下来着重讲的dynamicMeshDict以外，还有常见的transportPropertiesDict和turbulencePropertiesDict两个Dict文件。将运动粘度nu设置为水的粘度1e^-6，湍流模型设置为层流。对于dynamicMeshDict，主要讲解需要改动的部分：

（1）运动边界以及网格变动范围

如上图，WALL表示圆柱面对应的边界，如果自己的算例是别的名字，那么这里就得改。innerDistance代表以运动边界中心为圆心，网格变形的最小距离；outerDistance代表网格变形最大的距离。这两者之间就是一个网格变形环。这两个距离需要根据自己网格的实际尺寸来设置。如果想要全局变形，那么可以将innerDistance设置为0，outerDistance设置成超过网格最大尺寸的一个数。

（2）质量、惯量以及密度

质量由质量比、圆柱长度、圆柱直径和流体密度计算得到，公式在图里。质量比可以理解为固体的质量和其排开水的质量的比值，代表其在水里有多“重”。momentOfInertia为固体在这个惯性系里面三个方向的惯量，因为本算例不涉及转动，所以这个无关紧要，但不能设为0，否则会报错。密度设置为水的密度1000kg/m³。

（3）约束constraints

Constraints可以有约束的意思，且约束的程度更强，有限定的意思。在这里，constraints就是给圆柱的运动方向进行设定。如图所示，yline表示其只能在y方向上运动，zAxis则是不可绕z轴转动。这样，圆柱只能沿y轴进行平移。其它的constraints类型也在图中给出：fixedPoint可以把圆柱固定在某一点处；zPlane可以让圆柱在z平面移动，也就是x、y两个方向都可以运动。

（4） 约束restraints

Restraints约束的程度较弱，在这里，是给圆柱运动施加一定的阻力。Restraints有很多种，适合涡激振动的就是图中所示的linearSpring，也就是一个弹簧-阻尼系统。这个弹簧有定点（anchor）和动点（refAttachmentPt）两个end points。计算时，anchor始终固定，而refAttachmentPt则是动点的一个参考点，用来计算弹簧长度的变化量的。Stiffness、damping和restLength顾名思义分别是刚度、阻尼和弹簧原长，根据需要进行改动。至于应该怎么给？这些都是根据约化速度Ur来定的，在任何一篇涡激振动的文章里头都有对应的公式。

（5） system/fvSolution

fvScheme基本与圆柱绕流的一致。fvSolution则需要加上pcorr的矩阵求解器：

其它基本与pimpleFoam算例里的一致。

运行-后处理-使用脚本提取圆柱位移信息

运行

按自己机器的核数并行计算即可。需要注意，提取圆柱运动的位移，需要在运行时加上> log.pimpleFoam，这样，输出的信息就能存在log文件里面。

后处理

使用cat命令来提取圆柱的位移信息：

最后，圆柱的位移曲线就会存在plotfile.dat里面，就可以用matlab处理或者用tecplot进行查看。

计算结果

Ur=3~10共8组数据，与文献结果符合较好，说明之前的设置是比较好的。每个人算的结果都可能会有一点不一样，这与网格质量、数值格式、矩阵求解器设置等等相关。

分享一些经验

尽量用现实的变量进行计算而非无量纲数。写求解器计算时都会根据无量纲数设置一些参数，如将U、D设为1，那么nu就为Re的倒数；但在OpenFOAM里不推荐这样做。一方面，不可压求解器pimpleFoam等压力方程的p实际上是p/rho，在量纲上就会比较乱；另一方面，设计不同算例时，换算也容易出错。建议物理常数如nu、m、rho等都按照实际对应的物理量来设置，再根据公式来换算其他物理量。最后再对结果进行无量纲化，减小出错率。

直接改刚度，一个算例算N遍。在计算涡激振动时，一般都是固定了Re，改变圆柱刚度来改变Ur的。这意味着，流场的情况基本不变，只需要改变圆柱的刚度即可。在计算不同Ur值时，笔者是用了1个算例来算N次，期间不断修改dynamicMeshDict里面的stiffness。本次的算例继承了上一个算例的流场，节省了大量拉涡-形成卡门涡街-相对稳定的运算时间，提高实验效率。

来源：多相流在线