流固耦合动力学方法在工程中的应用

weixin

导读：现代航空、航天、船舶与海洋工程领域，航行器日益追求高速度，高机动性。工程装备追求质量轻，性能好。这些需求就使得航行器以及工程装备呈现出轻结构、大柔性的特点。当这些柔性结构在流体的作用下，就会诱发流固耦合振动的问题。在流体力的作用下，结构会发生弹性变形，而结构变形又改变了流场分布，这种相互耦合作用会使得弹性体的振动逐渐达到平衡状态，或者发散。这种发散现象就会导致结构的破坏。概括地讲，流固耦合动力学就是研究流体和结构相互耦合作用而产生的各种动力学问题。

晃动的赛格大厦

一、写在前面

以气体为流体介质的流固耦合问题我们称之为气动弹性问题，以水为流体介质的流固耦合问题我们称之为水弹性问题。许多的气动弹性、水弹性问题涉及流体力，弹性力和惯性力，这类问题称为动气动弹性或者动水弹性问题。另外一些气动弹性、水弹性问题只涉及流体力和弹性力，称之为静气动弹性或静水弹性问题。

随着计算机技术和数值方法的不断发展，流固耦合高保真仿真技术得到了大大提高，但计算代价依然很高，计算非常耗时。因此，建立适用于工程实用的流固耦合快速建模和仿真方法具有重要意义。

图1 柱体结构涡激振动二维、三维仿真模型

以柔性柱体结构、海洋立管涡激振动、水下航行器舵系统为例（图1、2、3），笔者从二维模型机理研究、到基于二维模型的流动控制、振动控制研究、到三维模型的动力学快速建模、快速仿真，再到高保真的全三维CFD/FEM双向耦合方法都做了一定的研究，总结了较多的仿真和科研经验。

图2 基于二维柱体涡激振动模型的振动模型、流动控制模型

图3 水下航行器舵系统的二维、三维流固耦合动力学模型

二、二维涡激振动模型如何建立？

我们知道二维弹性支撑柱体的数值模拟是研究海洋工程和风工程中柱体结构涡激振动现象和机理的重要手段。那么，二维涡激振动模型如何建立呢？

根据牛顿第二定律，2-DOF弹性支撑的柱体运动的控制方程可以写为：

式中，m 为圆柱体的质量，c 为结构阻尼系数，k 为结构刚度系数。此式又可以写为：

式中，柱体固有频率ω0=√k/m，阻尼比ζ=c/2√km。

对于涡激振动的数值模拟大多以柱体的横向振动研究为主，同时考虑横向和来流向的耦合振动研究相对较少。其主要原因在于，对柱体结构涡激振动数值模拟，必须保证柱体周围网格质量非常好，才能有效预测涡激振动响应。一般使用网格质量高的结构化网格，但是当柱体结构发生较大振动位移时，周围流场网格会发生畸变，甚至产生负网格，导致计算失败。

如果再考虑流向耦合振动，计算难度将非常大，计算成功率也将明显降低。如果采用非结构化网格，且采用网格重构技术，可以吸收柱体较大的振动位移，但是非结构化网格质量相比结构化网格质量较差，且采用非结构化网格必定需要大大增加网格量，从而大大增加了计算时间。

此外，如果对于涡激振动抑制装置设计研究，也就是在柱体结构表面形状变复杂的情况下，网格划分难度和计算难度也将大大增加。因此，寻求一种既可以保证网格质量，又能不大幅度增加网格数量，且可以避免网格畸变或者负网格问题的方法十分重要。

基于CFD商业软件FLUENT和结构动力学原理，通过用户自定义函数 (UDF) 及嵌套网格技术，建立了2-DOF弹性支撑柱体结构VIV数值模型。非定常不可压缩流体RANS方程为：

上式中，

式中，ρf 为不可压缩流体的密度；ui 表示i 方向上的瞬时速度分量，ui' 为i 方向上速度脉动量， 为速度的时间平均值；xi、t、p、μ 分别表示笛卡尔坐标系、时间、压力、运动粘度；μt 为湍流黏度，下标“t”表示湍流；kt 为湍动能；δij 是“Kronecker delta”符号，就是当i=j 时，δij =1，当i≠j 时，δij =0。湍流模型选用SST k–ω 湍流模型。通过计算流场，可以得到二维柱体表面的压力分布，进而可以得到作用在二维柱体上的升力和阻力系数：

结合方程

2-DOF弹性支撑的柱体运动的控制方程可以写为：

(a) 两自由度弹性支撑刚性柱体

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

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

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

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

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

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

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

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

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

初始条件为：

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

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

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

(a) 不同约化速度下的振幅分布

(b) 频率比随约化速度变化图

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

柱体涡激振动的最大幅值及频率比随不同约化速度的变化如图7所示，与实验数据对比，误差较小，验证了本文计算方法的正确性。从图7(a)可以看出，数值仿真出3种响应分支，当Ur 在3~4之间的时候，原始分支向上端分支转变；当Ur =5~6之间的时候，出现下端分支。在上端分支中振幅达到最大值0.98，而在下端分支中振幅最大值0.642；从Ur =11的时候开始，圆柱体的响应位移又回落到一个很小的数值。从图7(b)可以看出，在频率“锁定”区间Ur =4~10内，柱体的实际振动频率fv 与固定柱体的泄涡频率fst 分离，不再符合St 与Re数关系图。同时，柱体的实际涡泻频率fv 与柱体固有频率fn 比值稳定在1.15附近，而在解锁区域，柱体的实际振动频率fv 与固定柱体的涡脱频率fst 相同，这与前人的实验结果大致相同。

图8 不同时刻的涡量云图（周期T=3.34s）

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

那么，这个二维涡激振动的仿真方法我们同样可以用来设计涡激振动抑制装置。基于这种弹性振动的二维涡激振动模型，就可以合理地设计出涡激振动抑制装置的外形。比如图9为笔者参与设计的螺旋列板抑制装置结构。

图9 涡激振动抑制装置

再比如复杂系统的流固耦合动力学快速建模，往往采用一些势流理论可以试验快速计算。比如Theodorsen 时域水/气动力模型就比较适用于薄翼型结构的水动力/气动力结算，结合结构动力学相关理论，模态叠加就可以求解带有翼型结构的复杂系统的流固耦合振动响应。Theodorsen 时域模型为：

式中，

比如图10所示，为平板机翼时域模型仿真结果，就是采用Theodorsen 时域模型与平板的模态，进行模态叠加快速计算出振动响应。

(a) U=40m/s

(b) U=51m/s

(c) U=55m/s

图10 平板机翼时域模型仿真结果

三、流固耦合动力学的学习建议

同学们在学习流固耦合动力学方法、解决实际工程问题的时候不应一味追求高保真的全三维CFD/FEM双向耦合方法，这种方法计算量太大，短期内很难得到大量数据。建议：

  · 从二维模型开始做起，结合工程实际，合理简化模型，利用简单模型掌握基本方法、获得大量计算结果，总结经验，分析动力学机理。

  · 在掌握二维模型之后，同学们可以结合工程实际，建立势流理论/结构动力学方程全耦合的三维流固耦合动力学快速仿真模型，势流理论往往仿真速度快，精度往往符合工程需求，可以用于工程快速计算，解决工程实际问题。

以上所讲的两类方法比较适用于机理研究，模型优化，快速计算。全三维CFD/FEM双向耦合方法往往用于最后的模型验证，或者较少量的模型计算。总而言之，针对实际的工程问题，选用合适的流固耦合仿真方法是关键。

来源：仿真秀App微信公众号（ID：fangzhenxiu2018），作者：陈东阳。