Fluent 流固耦合基础教程（上）
Fluent 流固耦合基础教程（上）作者 Seventy721， 2012年2月
（本文版权属于作者，欢迎转载，但是请保留作者信息）
最近用fluent UDF做了个弹性梁的流固耦合问题。这里把经验跟大家分享一下。这个帖子主要是介绍基本操作和程序编写，算是扫盲贴；对于CFD，有限元、梁理论、流固耦合原理不做过多探讨，请大家见谅。另外作者水平有限，时间精力也有限，写得不好也请见谅了。
作者假设读者熟悉Fluent，对CFD的基本概念有所了解，具备一定的C/Fortran语言编程技能，并且了解有限元的基本理论和程序编写。不过即使背景有所欠缺也无妨，点是了解基本流程和操作方法。
1. 流固耦合问题的种类
流固耦合问题可以分为很多种。按照耦合程度分类，可以分为强耦合问题和弱耦合问题。但是这种分类方法有两个流行的版本。
第一版本根据流体和固体之间相互影响的程度来划分，在工业应用中谈得比较多。如果固体运动很小，对流场的影响不大，则认为是弱耦合。比如金属管道內的水流引起的管壁运动，机翼的振颤等，可以属于这类。在计算机技术不够发达的时代，这类问题一般采用简化的流场模型，如简单几何形状的理想流体，解析出由于固体运动而引起的流场作用力，然后再将这些作用力施加到固体控制方程上，从而得到附加质量，附加粘度，附加刚度，以及其他非线性项。随着计算机技术的发展，对这类问题的直接数值模拟逐渐成为可能。如果固体变形或运动比较大，其对流场的影响不可忽略，则属于强耦合。强耦合问题必须对流场和固体进行细致的计算，充分考虑固体变形对流场的影响。阀门开闭，血液通过血管瓣膜，旗帜在风中飘舞，等都属于这类问题。解决这类问题需要依靠数值模拟，传统的简化模型很难应用在这类问题上。
第二个版本根据流体和固体求解的模式来划分，在学术领域用的比较多。强耦合问题被认为是流场和固体变形必须同时求解的问题。带有渗流的多孔介质变形问题可以算这类问题。弱耦合问题则是流场和固体变形可以分开求解，但是二者的信息交换通过循环迭代完成。多数工程问题都可以算做这类问题。
这里我要讲的是第一类划分标准中的弱耦合问题，但所用的数值方法也可以应用到强耦合问题中，只是结构体的有限元方法需要采用大变形理论。按照第二种划分标准来说，这里解决的是流场和固体变形/运动分别求解，两场之间的相互作用通过迭代实现。
2. 流固耦合问题的数值模拟
既然是数值模拟，则需将系统方程在时域上离散。在每个时间步上分别解算流体域和固体域。固体的变形或位移导致流固边界的运动。这个运动以动网格的形势传递给流场，然后求解流场。流场的解则包括了这个运动所导致的反力。将这个反力反馈给固体，求解其下一时间步的响应。如图1所示。
http://forum.vibunion.com/data/attachment/album/201202/25/104743cxdvxpdwxbvaflxf.png
图1
以上的概念恐怕大家都清楚，只是谈到实现的时候，需要做很多细致的工作。当然，如果大家有ANSYS14 WorkBench，则流固耦合问题很容易解决，因为大部分的工作都自动化处理了。但是如果只用Fluent，则需要做一些工作。即便使用先进的数值模拟工具，如果能弄清楚整个流程以及其中的一些技巧，也会对提高效率以及保证精度有所帮助。下面就讲讲具体怎么做。
3. 流固耦合例题
既然是具体讲，就不如找个例题来说。图2是一根柔软的弹性梁，浸泡在在直径为Dt的圆截面管道里。梁的截面也是圆形，直径为D，长度为L。梁的两端为铰接。整个结构是三维结构。梁可以在三维空间里弯曲。流体为水，梁中点所在截面上的平均流速为U。管道入口处的流速为Uin。上游长度为Lup，下游长度为Ldw。上下游的长度分别为10D。
http://forum.vibunion.com/data/attachment/album/201202/25/104744pjjpngy5ybyjghgt.png
图2
4. 流体模型
流体模型的建立是关键的一环，因为流体模型的好坏直接影响到所得解的合理性。建立好的流体模型需要利用流体力学知识对所处理的问题进行综合分析。根据流体的雷诺数和其它特征建立合适的网格，选择合理的湍流模型和算法。我们希望流体模型能够尽量准确地捕捉到梁表面的力，因此准确的near-wall处理方法很重要。采用wall function还是精确地模拟边界层，需要给出适当的理由。由于我们处理的是动力问题，希望模型对梁表面力的变化给出准确的值。采用动态k-epsilon 或者k-omega湍流模型计算所得的激励幅值明显偏小，这是由于雷诺平均造成的。这里我们采用对流体动力特性捕捉比较好的大涡湍流模型(Large Eddy Simulation, LES)。LES对near-wall处理方法要求比较高，需要边界上的第一层网格厚度满足y+ ~ 1。在Gambit里可以方便地建立边界层网格，达到以上要求。如何建立良好的网格是个复杂的话题，这里且不涉及。
由于固体-液体边界要移动，流体网格必须改变，需要运用动网格技术。此时建模的时候需要将动网格区域和不动网格区域分开，如图3所示。
http://forum.vibunion.com/data/attachment/album/201202/25/104745xkcv0v4a9xvqx4vd.png
图3
固体和液体的交界面附近的区域是动网格区域（fluid_dyn），远离固体的区域可以定义为不动网格区域。动网格区域的大小应当视具体问题而定，但起码要能够囊括固体的最大位置变化范围。梁表面的网格节点位置随着梁的运动而变化，因此需要将梁的表面单独定义为一个区域 wall_mov_beam。梁的两个端点也需要进行控制，因此也可以分别定义独立区域：wall_mov_end1 和 wall_mov_end2。在划分流体网格的时候，接近固体表面的地方需要建立边界层。Fluent帮助文档上说边界层应该单独作为一个动网格区域，这是为了便于控制边界层的变形，防止出现边界层的异化。实际上如果将边界层的厚度增加到物理边界层厚度的两倍以上，并且只采用smoothing网格控制，则可以不必将边界层单独定义为一个区域。但是如果变形较大，流体网格需要重新划分，则最好将边界层单独分离为一个区域，并且在这个区域内使用smoothing控制。对于管壁附近的边界层，我们并不关心管壁上的压力和剪力，而且考虑到出了边界层以后这个区域的局部解对整体解影响不是很大，因此可以适当降低网格要求，甚至采用wall function近似逼近，这样可以节省不少计算量。下面两张图分别是图4几何模型（部分）和图5网格划分（部分）。
http://forum.vibunion.com/data/attachment/album/201202/25/104746h2dhjjs69jduzj4j.png
图4
http://forum.vibunion.com/data/attachment/album/201202/25/104748sseuoxeuo70oespu.png
图5
流体为不可压缩流体，求解器采用一次隐式瞬态算法。LES的动网格计算上Fluent 不支持二次格式，这就需要时间步设定要小一些。时间步长和求解器参数的选取也是个复杂的问题，这里也不展开讨论了。在这个算例中，我们取较大的步长，以加快计算速度。具体参数如下：
Time step size: 0.0005 s
Solver settings: unsteady / pressure-based / 1st order implicit
Viscous model: LES / Smagorinsky-Lilly / no dynamic stress
Solution settings:
velocity-pressure coupling: SIMPLE
Relaxation factor: pressure 0.3 / momentum 0.7 / density 1.0 / body force 1.0
Discretization: pressure standard / momentumn bounded central difference
在不加入流固耦合的情况下，计算结果收敛很好，稳定以后每个时间步上的循环（iteration）为两步。每个时间步循环终止时的残差为：
continuity x-velocity y-velocity z-velocity
4.2946e-04 8.3873e-06 8.1793e-06 1.3922e-04
流体模型解算成功之后就需要考虑动网格和固体变形的问题了。梁表面流体网格节点的运动需要利用UDF（User Defined Function）来控制。这将在下面的章节里说明。这里首先设置动网格的参数。
动网格区域fluid_dyn应该被设置为Deforming，采用的网格控制方法为Smoothing，但是对于复杂结构，可能需要Remeshing。具体选项的含义请参见Fluent 帮助文档中的用户手册第11章（Modeling Flows Using Sliding and Deforming Meshes）。这里只说明大体思路。这个区域的网格将会随着这个区域的边界（梁表面wall_mov_beam）变化。我们只要控制wall_mov_beam，则Fluent可以自动计算出这个区域内的流体网格的新位置。参数设置如图6所示。
http://forum.vibunion.com/data/attachment/album/201202/25/104749k37rn76f337y5pp7.png
图6
面区域wall_mov_end 1和wall_mov_end2也设置为动网格区域，跟fluid_dyn类似，具体的网格变形由Fluent自动处理。这是由于我们的梁是两端铰接的，因此端面上的网格变化不大。对于其他的情况，比如悬臂梁，这两个面上的点应该跟据梁的运动用UDF来控制，不能交给Fluent自动处理。 具体设置参数如图7所示。
http://forum.vibunion.com/data/attachment/album/201202/25/1047509xl9x117zl14318q.png
图7
梁表面网格节点的位移必须通过UDF来控制。因此在选择的时候应当选User-Defined。在Motion Attributes里面选择相应的UDF。这一点在下面的章节中说明。这里先给出Meshing Options的值，参看图8。这里提醒读者，wall_mov_beam是驱动网格变化的“源头”。我们通过UDF控制wall_mov_beam的运动，而周边的网格将随着wall_mov_beam的运动而改变。
http://forum.vibunion.com/data/attachment/album/201202/25/1047515pmpfg30s06s9p6c.png
图8
（～待续～）
内容预告：
5. 固体模型
6. 流固耦合UDF
6.1 力的传递
6.2网格控制
7. 结果的解读
8. 后记 真是营养贴，拜读，
赞！ 回复 1 # Seventy721 的帖子
非常感谢LZ，此贴必火，火速留名~ 谢谢楼主分享，学到些东西 怎么没有更新啊？期待~~ 抱歉，刚写好第5节。实在是没有时间写。最近一直在加班，有紧急任务。 拜读。。 受益匪浅！期待下部！ 赞！谢谢楼主分享。。。 下部写好了，在这里
http://forum.vibunion.com/thread-109534-1-1.html 强帖，收藏学习之 收藏了！楼主辛苦了！ 感谢 学习一下 好贴留名啊。。。。。。 拜读了，楼主很强