Fluent

Fluent流固耦合基础教程（上）作者Seventy721，2012年2月（本文版权属于作者，欢迎转载，但是请保留作者信息）最近用fluentUDF做了个弹性梁的流固耦合问题。这里把经验跟大家分享一下。这个帖子主要是介绍基本操作和程序编写，算是扫盲贴；对于CFD，有限元、梁理论、流固耦合原理不做过多探讨，请大家见谅。另外作者水平有限，时间精力也有限，写得不好也请见谅了。作者假设读者熟悉Fluent，对CFD的基本概念有所了解，具备一定的C/Fortran语言编程技能，并且了解有限元的基本理论和程序编写。不过即使背景有所欠缺也无妨，点是了解基本流程和操作方法。1.流固耦合问题的种类流固耦合问题可以分为很多种。按照耦合程度分类，可以分为强耦合问题和弱耦合问题。但是这种分类方法有两个流行的版本。第一版本根据流体和固体之间相互影响的程度来划分，在工业应用中谈得比较多。如果固体运动很小，对流场的影响不大，则认为是弱耦合。比如金属管道內的水流引起的管壁运动，机翼的振颤等，可以属于这类。在计算机技术不够发达的时代，这类问题一般采用简化的流场模型，如简单几何形状的理想流体，解析出由于固体运动而引起的流场作用力，然后再将这些作用力施加到固体控制方程上，从而得到附加质量，附加粘度，附加刚度，以及其他非线性项。随着计算机技术的发展，对这类问题的直接数值模拟逐渐成为可能。如果固体变形或运动比较大，其对流场的影响不可忽略，则属于强耦合。强耦合问题必须对流场和固体进行细致的计算，充分考虑固体变形对流场的影响。阀门开闭，血液通过血管瓣膜，旗帜在风中飘舞，等都属于这类问题。解决这类问题需要依靠数值模拟，传统的简化模型很难应用在这类问题上。第二个版本根据流体和固体求解的模式来划分，在学术领域用的比较多。强耦合问题被认为是流场和固体变形必须同时求解的问题。带有渗流的多孔介质变形问题可以算这类问题。弱耦合问题则是流场和固体变形可以分开求解，但是二者的信息交换通过循环迭代完成。多数工程问题都可以算做这类问题。这里我要讲的是第一类划分标准中的弱耦合问题，但所用的数值方法也可以应用到强耦合问题中，只是结构体的有限元方法需要采用大变形理论。按照第二种划分标准来说，这里解决的是流场和固体变形/运动分别求解，两场之间的相互作用通过迭代实现。2.流固耦合问题的数值模拟既然是数值模拟，则需将系统方程在时域上离散。在每个时间步上分别解算流体域和固体域。固体的变形或位移导致流固边界的运动。这个运动以动网格的形势传递给流场，然后求解流场。流场的解则包括了这个运动所导致的反力。将这个反力反馈给固体，求解其下一时间步的响应。如图1所示。图1以上的概念恐怕大家都清楚，只是谈到实现的时候，需要做很多细致的工作。当然，如果大家有ANSYS14WorkBench，则流固耦合问题很容易解决，因为大部分的工作都自动化处理了。但是如果只用Fluent，则需要做一些工作。即便使用先进的数值模拟工具，如果能弄清楚整个流程以及其中的一些技巧，也会对提高效率以及保证精度有所帮助。下面就讲讲具体怎么做。3.流固耦合例题既然是具体讲，就不如找个例题来说。图2是一根柔软的弹性梁，浸泡在在直径为Dt的圆截面管道里。梁的截面也是圆形，直径为D，长度为L。梁的两端为铰接。整个结构是三维结构。梁可以在三维空间里弯曲。流体为水，梁中点所在截面上的平均流速为U。管道入口处的流速为Uin。上游长度为Lup，下游长度为Ldw。上下游的长度分别为10D。图24.流体模型流体模型的建立是关键的一环，因为流体模型的好坏直接影响到所得解的合理性。建立好的流体模型需要利用流体力学知识对所处理的问题进行综合分析。根据流体的雷诺数和其它特征建立合适的网格，选择合理的湍流模型和算法。我们希望流体模型能够尽量准确地捕捉到梁表面的力，因此准确的near-wall处理方法很重要。采用wallfunction还是精确地模拟边界层，需要给出适当的理由。由于我们处理的是动力问题，希望模型对梁表面力的变化给出准确的值。采用动态k-epsilon或者k-omega湍流模型计算所得的激励幅值明显偏小，这是由于雷诺平均造成的。这里我们采用对流体动力特性捕捉比较好的大涡湍流模型(LargeEddySimulation,LES)。LES对near-wall处理方法要求比较高，需要边界上的第一层网格厚度满足y+~1。在Gambit里可以方便地建立边界层网格，达到以上要求。如何建立良好的网格是个复杂的话题，这里且不涉及。由于固体-液体边界要移动，流体网格必须改变，需要运用动网格技术。此时建模的时候需要将动网格区域和不动网格区域分开，如图3所示。图3固体和液体的交界面附近的区域是动网格区域（fluid_dyn），远离固体的区域可以定义为不动网格区域。动网格区域的大小应当视具体问题而定，但起码要能够囊括固体的最大位置变化范围。梁表面的网格节点位置随着梁的运动而变化，因此需要将梁的表面单独定义为一个区域wall_mov_beam。梁的两个端点也需要进行控制，因此也可以分别定义独立区域：wall_mov_end1和wall_mov_end2。在划分流体网格的时候，接近固体表面的地方需要建立边界层。Fluent帮助文档上说边界层应该单独作为一个动网格区域，这是为了便于控制边界层的变形，防止出现边界层的异化。实际上如果将边界层的厚度增加到物理边界层厚度的两倍以上，并且只采用smoothing网格控制，则可以不必将边界层单独定义为一个区域。但是如果变形较大，流体网格需要重新划分，则最好将边界层单独分离为一个区域，并且在这个区域内使用smoothing控制。对于管壁附近的边界层，我们并不关心管壁上的压力和剪力，而且考虑到出了边界层以后这个区域的局部解对整体解影响不是很大，因此可以适当降低网格要求，甚至采用wallfunction近似逼近，这样可以节省不少计算量。下面两张图分别是图4几何模型（部分）和图5网格划分（部分）。图4图5流体为不可压缩流体，求解器采用一次隐式瞬态算法。LES的动网格计算上Fluent不支持二次格式，这就需要时间步设定要小一些。时间步长和求解器参数的选取也是个复杂的问题，这里也不展开讨论了。在这个算例中，我们取较大的步长，以加快计算速度。具体参数如下：Timestepsize:0.0005sSolversettings:unsteady/pressure-based/1storderimplicitViscousmodel:LES/Smagorinsky-Lilly/nodynamicstressSolutionsettings:velocity-pressurecoupling:SIMPLERelaxationfactor:pressure0.3/momentum0.7/density1.0/bodyforce1.0Discretization:pressurestandard/momentumnboundedcentraldifference在不加入流固耦合的情况下，计算结果收敛很好，稳定以后每个时间步上的循环（iteration）为两步。每个时间步循环终止时的残差为：continuityx-velocityy-velocityz-velocity4.2946e-048.3873e-068.1793e-061.3922e-04流体模型解算成功之后就需要考虑动网格和固体变形的问题了。梁表面流体网格节点的运动需要利用UDF（UserDefinedFunction）来控制。这将在下面的章节里说明。这里首先设置动网格的参数。动网格区域fluid_dyn应该被设置为Deforming，采用的网格控制方法为Smoothing，但是对于复杂结构，可能需要Remeshing。具体选项的含义请参见Fluent帮助文档中的用户手册第11章（ModelingFlowsUsingSlidingandDeformingMeshes）。这里只说明大体思路。这个区域的网格将会随着这个区域的边界（梁表面wall_mov_beam）变化。我们只要控制wall_mov_beam，则Fluent可以自动计算出这个区域内的流体网格的新位置。参数设置如图6所示。图6面区域wall_mov_end1和wall_mov_end2也设置为动网格区域，跟fluid_dyn类似，具体的网格变形由Fluent自动处理。这是由于我们的梁是两端铰接的，因此端面上的网格变化不大。对于其他的情况，比如悬臂梁，这两个面上的点应该跟据梁的运动用UDF来控制，不能交给Fluent自动处理。具体设置参数如图7所示。图7梁表面网格节点的位移必须通过UDF来控制。因此在选择的时候应当选User-Defined。在MotionAttributes里面选择相应的UDF。这一点在下面的章节中说明。这里先给出MeshingOptions的值，参看图8。这里提醒读者，wall_mov_beam是驱动网格变化的“源头”。我们通过UDF控制wall_mov_beam的运动，而周边的网格将随着wall_mov_beam的运动而改变。图85.固体模型固体变形和运动的求解可以采用很多方法。对于简单的问题，可以采用解析解。目前的这个弹性梁的问题，就可以用采用解析解的方法。对于复杂结构，有限元是比较常用的方法。建立有限元模型有两种方案。一种是采用三维实体单元完整地模拟出整个梁。这种方法的好处是流体和固体的接触面两边都是实际网格，一边是流体网格，一边是有限元网格。边界上的位移，速度，受力等计算都比较简单，缺点是固体模型计算量大。第二种方法是采用结构体单元，如梁，板壳等。固体模型简单，但是接触面上的处理稍微复杂一些。结构体单元在工程上应用得比较普遍。为了说明如何利用结构体单元做流固耦合，这里采用三节点的梁单元。梁单元是一维单元，单元的几何形状只是一条中性线，没有体积。梁的变形是由中性线的位移和梁截面的转动描述的。在流体模型中，梁的体积是存在的。梁表面上的流体网格节点的位置需要由梁的变形确定。因此一个关键的步骤是根据梁的变形计算出梁表面的流体网格节点的位置。这里的梁采用小变形条件下的欧拉梁。根据欧拉梁理论，梁截面为刚性面且保持与中性线垂直。图9所示为一梁单元。梁表面上有一流体网格节点A。点A所在的截面轴向坐标为ξ。则流体网格节点的位置可以表示为其中uA是A点的位移向量，uO是O点的位移向量。r是O点到A点的位置向量。Φ是旋转矩阵。对于一般情况，Φ矩阵的分量是欧拉转角的函数，在小变形小转动条件下，可以简化为O点处梁截面的转角的函数:O点处的位移和梁截面的转角可以通过有限元的形函数求得。图9图10给出了梁变形后其表面上流体网格节点的位置。两个端点处的流体网格没有加以控制，而是交给fluent自动处理，这只是对当前这种支撑情况有效。对于一般情况，如悬臂梁，则需要人工计算。图10与网格位置计算类似。流体计算所得到的力需要传递给固体模型。梁表面的流体网格上的压力和剪力也需要利用有限元的形函数离散到有限元节点上。有限元方面的基本知识，我在这里就不罗嗦了。任何一本教材上都写得很清楚。流固耦合问题是瞬态问题，因此需进行时间积分。常用的时间积分算法有Newmark，Wilson-theta，Runge-Kutta等。Newmark方法在多自由度的线性问题中应用比较普遍。这里我们也采用Newmark方法。但是在程序编写的时候需要注意Fluent求解流固耦合问题的流程。Fluent必须作为整个过程的主导程序，如图11所示。图11这种流程模式给程序编写带来一些麻烦，因为Newmark算法的循环被拆解开进行，必须按照单个时间步考虑。这就要求每个时间步之间的数据传递不能用通常的变量存储。解决的办法有两个：一个办法是用硬盘存储，但是这样耽误在文件I/O上的时间很多；另外一个办法是利用常驻内存的global变量存储，但是具体操作要看所用的编程语言环境。这里我用Fortran90编写的有限元程序，global变量保存在独立的module里。如果只用UDF的话，跟C语言一致。本帖最后由