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[1] = arrPForce[1] + arrVForce[1];
arrForce[2] = arrPForce[2] + arrVForce[2];
arrForce[3] = arrPForce[3] + arrVForce[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[0]=0.0;
arrDisp[1]=0.0;
arrDisp[2]=0.0;
arrVelo[0]=0.0;
arrVelo[1]=0.0;
arrVelofor(i=1;i<=countF;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_N{
arrForce[1] = arrPForce[1] + arrVForce[1];
arrForce[2] = arrPForce[2] + arrVForce[2];
arrForce[3] = arrPForce[3] + arrVForce[3];
}
|
[复制链接]
发表于 2012-3-6 08:50
提升卡
变色卡
显身卡
楼主
发表于 2012-3-6 09:20
强帖，拜读下，LZ辛苦了
楼主厉害！