Fluent 流固耦合基础教程（下）

作者 Seventy721， 2012年3月

（本文版权属于作者，欢迎转载，但是请保留作者信息）

5. 固体模型

固体变形和运动的求解可以采用很多方法。对于简单的问题，可以采用解析解。目前的这个弹性梁的问题，就可以用采用解析解的方法。对于复杂结构，有限元是比较常用的方法。建立有限元模型有两种方案。一种是采用三维实体单元完整地模拟出整个梁。这种方法的好处是流体和固体的接触面两边都是实际网格，一边是流体网格，一边是有限元网格。边界上的位移，速度，受力等计算都比较简单，缺点是固体模型计算量大。第二种方法是采用结构体单元，如梁，板壳等。固体模型简单，但是接触面上的处理稍微复杂一些。结构体单元在工程上应用得比较普遍。

为了说明如何利用结构体单元做流固耦合，这里采用三节点的梁单元。梁单元是一维单元，单元的几何形状只是一条中性线，没有体积。梁的变形是由中性线的位移和梁截面的转动描述的。在流体模型中，梁的体积是存在的。梁表面上的流体网格节点的位置需要由梁的变形确定。因此一个关键的步骤是根据梁的变形计算出梁表面的流体网格节点的位置。

这里的梁采用小变形条件下的欧拉梁。根据欧拉梁理论，梁截面为刚性面且保持与中性线垂直。图9所示为一梁单元。梁表面上有一流体网格节点A。点A所在的截面轴向坐标为ξ0。则流体网格节点的位置可以表示为

其中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语言一致。

6. UDF

UDF是用Fluent做流固耦合的关键。UDF是Fluent 用来提供可扩展功能的框架。做过windows或者linux程序开发的朋友会觉得很好理解，但不熟悉编程的朋友可能会觉得费解。在这里我简单解释一下它的意思。Fluent 是个编译好的程序，想要实现功能扩展，最方便的办法就是利用动态链接库。在程序运行的时候将指定的动态链接库调入内存，然后通过预先定义好的接口函数执行用户定义的子程序。换句话说UDF就是一些放在动态链接库里的函数。这些函数的名字和格式已经由Fluent预先定义好，但是内容是空的，需要用户来写。Fluent 会在预定的情况下调用指定的函数，去运行用户定义的代码。所以对用户来说，就是填写一个指定的函数，编译成动态链接库，在Fluent里链接上，然后等着Fluent来调用。如此简单。

这些预先指定好的函数由被称为定义宏（DEFINE macro）的宏命令来定义。这是个很拗口的说法，不过一般用户不必深究，拿它当函数来用就是。为了简单起见，下面用“UDF函数”来代替“定义宏”。 Fluent 提供了很多UDF函数，具体有哪些，都是什么功能，诸位可以看Fluent UDF手册。这里就只说跟流固耦合有关的一个。下面的UDF函数用来定义网格的运动。

DEFINE_GRID_MOTION( name, d, dt, time, dtime)

除了编译的方法，Fluent 还支持对UDF的逐行解释执行。这样做方便调试，但是动网格的一些UDF函数是不能这样做的，所以我们还得用动态链接库。

UDF是C语言编写的。Fluent 自带一个C编译器和一堆头文件。UDF的编译可以在Fluent GUI里自动进行，但是也可以手工完成。有些情况下必须要手工操作，比如有C/Fortran混合编程的时候。顺便说一句，混合程序的编译也是个头疼的问题，需要费很多周章。这跟所用的操作系统和Fortran 版本有关。如果所写的UDF比较简单，只包括普通的C文件，则自动编译就可以。如果需要手工操作，则要建立如下文件结构

work folder （工作目录，模型放这里）
  |-> libudf （UDF目录，有个Makefile）
  |-> src （UDF源文件放这里）
  |-> ultra_64 （这个可能根据所用操作系统有所不同）
  |-> 3ddp （3ddp单机版）
  |-> 3ddp_host（3ddp并行版主节点）
  |-> 3ddp_node（3ddp并行版从节点）

然后修改src/makefile文件。之后回到libudf目录执行make。

Fluent 的数据结构是cell - face - node，如图12所示。每个网格的数据点是中心点 (cell center)。

图12

6.1 UDF - DEFINE_GRID_MOTION

这一节讲一下DEFINE_GRID_MOTION 的大体结构。

// 首先引用几个头文件。这几个头文件在Fluent的src 目录下。

#include "udf.h"

#include "dynamesh_tools.h"

#include "sg.h"

// 定义计算梁响应的有限元函数。此函数是Fortran函数。

extern void beam_response_(

  int *nfgrid,

  double (*fgrid)[ARR_LIMIT_FACE_BEAM][3],

  double (*force)[ARR_LIMIT_FACE_BEAM][3],

  double *gamm,

  double *beta,

  double *dt,

  double *t,

  int *run,

  int *ndgrid,

  double (*dgrid)[ARR_LIMIT_FACE_BEAM*NNODE_ON_FACE][3],

  double (*disp)[ARR_LIMIT_FACE_BEAM*NNODE_ON_FACE][3],

  double (*velo)[ARR_LIMIT_FACE_BEAM*NNODE_ON_FACE][3],

  int *info

);

// ... (定义其它全局变量和函数)

/* ------------------------------------------------------------------------- */

/* --------------------- macros are defined below here --------------------- */

/* ------------------------------------------------------------------------- */

/* Macro for defining dynamic mesh udf for the beam */

DEFINE_GRID_MOTION(beam,domain,dt,time,dtime)

{

  // beam - UDF 的名字

  // domain - 当前的domain数据结构，存储了有关网格的信息，但不是网格本身

  // dt - 指向网格数据结构（thead）的指针

  // time - 当前时刻

  // dtime - 时间步长

  // Define variables: pointers and utilities

  Thread *t = DT_THREAD(dt);

  cell_t c;

  face_t f;

  Node *node;

  int idNodeL; // local index of a node inside a face

  int countF; // number of faces

  int countN; // number of nodes

  int index;

  int i,j,k;

  int run;

  int id;

  // 定义计算结构响应需要的变量

  double xi; // the axial coordinate of a grid node

  double area; // area

  double pressure; // pressure

  double distance; // distance between two points

  double a_by_es; // A'A/(A'e)

  double gamm=0.5;

  double beta=0.25;

  // 定义主/从节点间数据传递用的数组... （省略）

  // Define constants

  const double PI=3.14159265358979323846;

  const double VISCOSITY = 0.001;

  const double DENSITY = 1000.0;

  const double TOLERANCE = 1e-15;

  const double LOWERLIMIT = 1e-8; 

  // 向量初始化 ... (省略)

  // 这一段用来取得梁表面的流体网格节点位置

  #if !RP_HOST // SERIAL OR NODE

  countF=0;

  countN=0;

  begin_f_loop(f,t) // 扫描当前domain的所有face

  {

    if PRINCIPAL_FACE_P(f,t)

    {

      countF++;

      if(countF>ARR_LIMIT_FACE_BEAM)

      {

        // 错误处理 ... (省略)

      }

      f_node_loop(f,t,idNodeL) // 扫描当前face的所有node

      {

         node=F_NODE(f,t,idNodeL);

         // NODE_X, NODE_Y, NODE_Z存储了节点坐标，但是注意node不是整数类型，实际上是联合变量。

         arrNode[countN][0]=NODE_X(node);

         arrNode[countN][1]=NODE_Y(node);

         arrNode[countN][2]=NODE_Z(node);

         countN++;

       }

    }

  }

  end_f_loop(f,t);

  #endif

  // 将数据传输到主节点 ... (省略)

  // 将相邻区域设置为动网格

  #if !RP_HOST // SERIAL OR NODE

  SET_DEFORMING_THREAD_FLAG(THREAD_T0(t));

  #endif

  // 计算表面力（见6.2节）

  // 将数据传输到主节点 ... (省略)

  // 计算梁的响应（见6.3节）

  // 将数据传输到计算节点 ... (省略)

  // 更新网格 （见6.4节）

}

6.2 力的计算

梁表面所受到的激励分为法向力和切向力。法向力由表面压力引起；切向力就是剪力。梁的表面覆盖着流体网格的面(face)。数据的计算在面中点(face centroid)上进行。法向力的大小等于压力乘以网格面积。剪力等于剪应力乘以网格面积。剪应力等于粘度乘以速度的导数。速度的导数可以简单近似为为cell centre velocity 除以cell centre 到 wall的距离。当然也可以采用更高阶的近似方法。需要注意的是在并行系统上，网格被分为多个区，每个区之间的交界面上的face会被重复计算。为了防止这种情况发生，需要用PRINCIPAL_FACE_P判断是否为该face实际存在于当前的区里。

// Obtain the centroid location and the force on the centroids

index=0;

begin_f_loop(f,t)

{

  if PRINCIPAL_FACE_P(f,t)

  {

  // Save the face id

  arrFaceID[index]=f;

  // Obtain the centroid coordinates and save in arrCentr

  F_CENTROID(vecXf,f,t);

  arrCentr[index][0]=vecXf[0]; // x

  arrCentr[index][1]=vecXf[1]; // y

  arrCentr[index][2]=vecXf[2]; // z

  // Obtain the area vector and the area

  F_AREA(vecArea,f,t);

  area=sqrt(pow(vecArea[0],2)+pow(vecArea[1],2)+pow(vecArea[2],2));

  // Obtain the pressure and calculate the pressure force

  pressure=F_P(f,t);

  vecLift[0]=vecArea[0]*pressure;

  vecLift[1]=vecArea[1]*pressure;

  vecLift[2]=vecArea[2]*pressure;

  arrPForce[index][0]=vecLift[0];

  arrPForce[index][1]=vecLift[1];

  arrPForce[index][2]=vecLift[2];

  // Obtain the shear stress and calculate the viscous force

  // get the face parameters

  BOUNDARY_FACE_GEOMETRY(f,t,vecArea,distance,vecEs,a_by_es,vecDr0)

  if(distance < TOLERANCE) // distance is always positive

  {

    distance=LOWERLIMIT;

  }

  // -- get the centroid velocity of the cell

  vecVelo[0]=C_U(c,t);

  vecVelo[1]=C_V(c,t);

  vecVelo[2]=C_W(c,t);

  // -- calculate the viscous force (= shear stress * area)

  vecDrag[0]=(VISCOSITY/distance)*vecVelo[0]*area;

  vecDrag[1]=(VISCOSITY/distance)*vecVelo[1]*area;

  vecDrag[2]=(VISCOSITY/distance)*vecVelo[2]*area;

  arrVForce[index][0]=vecDrag[0];

  arrVForce[index][1]=vecDrag[1];

  arrVForce[index][2]=vecDrag[2];

  // Increase index

  index=index+1;

  }

}

end_f_loop(f,t);

// Calculate the total force

for(i=1;i<=countF;i++)

{

  arrForce[i][1] = arrPForce[i][1] + arrVForce[i][1];

  arrForce[i][2] = arrPForce[i][2] + arrVForce[i][2];

  arrForce[i][3] = arrPForce[i][3] + arrVForce[i][3];

}

6.3结构响应

结构响应由Fortran函数求得，采用Newmark时间积分。注意Fortran 函数的末尾要加上一个额外的下划线。另外变量名前要加&，因为Fortran函数参数都是地址传递而非值传递。

// Perform the Newmark scheme.

if (flagInitial!=1)

{

  // Calculate the beam response

  info=0;

  for (i=1;i<=countN;i++)

  {

    arrDisp[i][0]=0.0;

    arrDisp[i][1]=0.0;

    arrDisp[i][2]=0.0;

    arrVelo[i][0]=0.0;

    arrVelo[i][1]=0.0;

    arrVelo[i][2]=0.0;

  }

  run=1;

  beam_response_(

  &countF, // nfgrid,

  &arrCentr, // fgrid,

  &arrForce, // force,

  &gamm, // gamma=0.5,

  &beta, // beta=0.25,

  &dtime, // dt,

  &time, // t,

  &run, // run, 0-preparation only, 1-run

  &countN, // ndgrid,

  &arrNode, // dgrid,

  &arrDisp, // disp,

  &arrVelo, // velo,

  &info // info

  );

}

6.4网格更新

得到结构响应以后，需要更新梁表面流体网格节点的位置。可以采用NV_V命令赋值。NODE_COORD(node)对应的是节点node的坐标。由于节点循环的时候会重复访问两个单元共有的节点，因此更新的时候需要用NODE_POS_NEED_UPDATE检查当前节点是否已经被更新过。

if (flagInitial!=1)

{

  // Update the grid nodes

  index=0;

  begin_f_loop(f,t)

  {

    if PRINCIPAL_FACE_P(f,t)

    {

      f_node_loop(f,t,idNodeL)

      {

        node=F_NODE(f,t,idNodeL); NV_V(NODE_COORD(node)

        // Update node if the current node has not been previously

        // visited when looping through previous faces

        if (NODE_POS_NEED_UPDATE (node))

        {

          // Indicate that node position has been update

          // so that it's not updated more than once

          NODE_POS_UPDATED(node);

          // Retrieve the nodel displacement

          vecNew[0]=arrNode[index][0]; // x

          vecNew[1]=arrNode[index][1]; // y

          vecNew[2]=arrNode[index][2]; // z

          // Update the fluid grid node

          NV_V(NODE_COORD(node),=,vecNew);

        }

        index++;

      }

    }

  }

  end_f_loop(f,t);

}

7. 计算结果

为了保证模型的正确性，需要查看数值模拟的结果。这里不讲如何诠释结果，只讲如何检查流固耦合过程是否正确。最直接的方法是观察网格的变形。简单的做法是在Write/Autosave中定义自动保存的频率。一般按照时间步来保存，比如选择每1000个时间步保存一次。然后查看每次保存的时刻上网格的变形程度。保存的时候一定要将case和data文件都保存下来。另外保存频率的选择要保证一个振动周期内保存至少5次以上。这样才能看出周期运动。图13是结构网格变形的例子。

图13

查看流体网格是比较麻烦的工作，但也是必须的工作。如果只查看固体变形的过程，并不能说明流场网格也正确变化。笔者第开始运行程序的时候就遇到这样的情况。；流体网格没有更新，但是流体力被正确地传递给，而固体振动也是正常的。这样的结果得到是接耦合的流致振动解，比如湍流激励引起的振动。但是这样的解是不正确的，因为没有包括附加质量，流体阻尼，和附加刚度。所以一个很重要的方法是观察结构的响应，然后根据响应分析这些附加质量等效果是否存在。对于这个梁，没有附加质量的自由振动周期为0.4秒左右，带有附加质量以后上升到0.57秒左右。与理论估计得到的结果一致。

图14

8. 后记

用每天挤牙膏的方式，终于写完了。水平有限，时间有限，只能写到这个层次了。这篇文章写得很粗略，只是作为大概的介绍，很多细致的内容还要读者自己去探索。流固耦合还是个很有意思的话题的。现在ANSYS做得好，恐怕大家以后都不用这么费劲写UDF了。耦合问题现在主要需要面对的是计算量问题。解决实际工程问题，没有并行计算恐怕很难完成。

有限体积法和有限元已经都是很成熟的东西。关键是湍流问题不好解决。大涡理论是个很好的方向，只是near wall条件不好处理。如果能做出好的wall function 就可以大大降低计算量。流固耦合的难题还是在强耦合问题上，据说ADINA做得不错。

希望这篇文章对打算进入流固耦合领域的朋友有点启发，起码算是个快速入门。这篇文章里的代码是具有代表性的段落，但不是完整的程序，仅供读者参考，请不要不加思考地照搬。另外请不要找我要完整代码，我是不会给你的。做这个东西，我是有钱挣的。坦白地说一共八千。如果你非要完整代码也可以，我给你打到极限折，只收十分之一。顺便说一句，货币单位是美元。

朋友们好好学，这年头知识就是钱。