7 自定义函数(UDF)
7.1，概述
用户自定义函数（User－Defined Functions，即UDFs）可以提高FLUENT程序的标准计算功能。它是用C语言书写的，有两种执行方式：interpreted型和compiled型。Interpreted型比较容易使用，但是可使用代码（C语言的函数等）和运行速度有限制。Compiled型运行速度快，而且也没有代码使用范围的限制，但使用略为繁琐。
我们可以用UDFs来定义：
a) 边界条件 b) 源项
c) 物性定义（除了比热外） d) 表面和体积反应速率 e) 用户自定义标量输运方程
f) 离散相模型（例如体积力，拉力，源项等）
g) 代数滑流（algebraic slip）混合物模型（滑流速度和微粒尺寸） h) 变量初始化 i) 壁面热流量
j) 使用用户自定义标量后处理
边界条件UDFs能够产生依赖于时间，位移和流场变量相关的边界条件。例如，我们可以定义依赖于流动时间的x方向的速度入口，或定义依赖于位置的温度边界。边界条件剖面UDFs用宏DEFINE_PROFILE定义。有关例子可以在5.1和6.1中找到。源项UDFs可以定义除了DO辐射模型之外的任意输运方程的源项。它用宏DEFINE_SOURCE定义。有关例子在5.2和6.2中可以找到。物性UDFs可用来定义物质的物理性质，除了比热之外，其它物性参数都可以定义。例如，我们可以定义依赖于温度的粘性系数。它用宏DEFINE_PROPERTY定义，相关例子在6.3中。反应速率UDFs用来定义表面或体积反应的反应速率，分别用宏DEFINE_SR_RATE和DEFINE_VR_RATE定义，例子见6.4。离散相模型用宏DEFINE_DPM定义相关参数，见5.4。UDFs还可以对任意用户自定义标量的输运方程进行初始化，定义壁面热流量，或计算存贮变量值（用用户自定义标量或用户自定义内存量）使之用于后处理。相关的应用见于5.3，5.5，5.6和 5.7。
UDFs有着广泛的应用，本文并不能一一叙述。如果在使用中遇到问题，可以联系FLUENT技术支部门要求帮助。在此推荐一个网站www.cfd-online.com，上面有FLUENT论坛，可进行相关询问和讨论。 7.1.1 书写UDFs的基本步骤
在使用UDFs处理FLUENT模型的过程中，我们一般按照下面五步进行： 1． 概念上函数设计 2． 使用C语言书写 3． 编译调试C程序 4． 执行UDF
5． 分析与比较结果
第一步分析我们所处理的模型，目的是得到我们要书写的UDF的数学表达式。第二步将数学表达式转化成C语言源代码。第三步编译调试C语言源代码。第四步在FLUENT中执行UDF。最后一步，将所得到的结果与我们要求的进行比较，如果不满足要求，则需要重复上面的步骤，直到与我们期望的吻合为止。 7.1.2 Interpreted型与Compiled型比较
Compiled UDFs执行的是机器语言，这和FLUENT本身运行的方式是一样
的。一个叫做Makefile的过程能够激活C编辑器，编译我们的C语言代码，从而建立一个目标代码库，目标代码库中包含有高级C语言的低级机器语言诠释。在运行的时候，一个叫做“dynamic loading”的过程将目标代码库与FLUENT连接。一旦连接之后，连接关系就会在case文件中与目标代码库一起保存，所以读入case文件时，FLUENT就会自动加载与目标代码库的连接。这些库的建立是基于特定计算机和特定FLUENT版本的，所以升级FLUENT版本后，就必须重新建立相应的库。
相反，Interpreted UDFs是在运行的时候直接装载编译C语言代码的。在这种情况下，生成的机器代码不依赖于计算机和FLUENT版本。编译后，函数信息将会保存在case文件中，所以读入case文件时，FLUENT也会自动加载相应的函数。Interpreted UDFs具有较强的可移植性，而且编译比较简单。对于简单的UDFs，如果对运行速度要求不高，一般就采用Interpreted型的。
下面列出的是两种UDFs的一些特性：
1． Interpreted UDFs
79
——独立于计算机结构；
——能够完全当作Compiled UDFs使用； ——不能与其它编译系统或用户库连接；
——只支持部分C语言，不能包含：goto语句，非ANSI－C语法，结构，联合，函数指针，函数数组等。 ! Interpreted UDFs能够使用FLUENT提供的宏，间接引用存贮于FLUENT的变量，详见2.10。
2． Compiled UDFs
——运行速度比Interpreted UDFs快； ——能够完全于C语言结合；
——能够用任何兼容ANSI－C的编辑器编译；
——对不同FLUENT版本（2D或3D）需要建立不同的共享库； ——不一定能够当作Interpreted UDFs使用。
我们在使用中首先要根据具体情况，明确使用哪种UDFs，然后才能进一步去实现，在使用中要注意上面叙述的事项。
第二节 书写 UDFs
7.2.1 概述
书写Interpreted型和Compiled型用户自定义函数的过程和书写格式是一样的。主要的区别在于与C语言的结合程度，Compiled型能够完全使用C语言的语法，而Interpreted型只能使用其中一小部分，这在前面有过论述。 7.2.2 UDF 格式
通用的UDF格式由三部分组成：
1． 定义恒定常数和包含库文件，分别由＃DEFINE和＃INCLUDE陈述（见2.3）； 2． 宏DEFINE_* 定义UDF函数（见2.4）； 3． 函数体部分（见2.9）
包含库有udf.h，sg.h，mem.h，prop.h，dpm.h等，其中udf.h是必不可少的，书写格式为#include “udf.h”；所有数值都应采用SI单位制；函数体部分字母采用小写，Interpreted型只能够包含FLUENT支持的C语言语法和函数。 7.2.3 包含库udf.h
库文件udf.h必须C函数开头包含。 7.2.4 定义函数 7.2.4.1 简介
Fluent 公司提供了一套宏，来帮助我们定义函数。这些宏都以DEFINE_开始，对它们的解释包含在udf.h文件中，所以我们必需要包含库udf.h。为了方便使用，我们把对udf.h文件中解释宏列在附录A中。
UDF使用宏DEFINE_定义，括号列表中第一个参数代表函数名。例如 DEFINE_PROFILE(inlet_x_velocity，thread，position) 定义了一个名为inlet_x_velocity的函数。 ! 所有函数名必须小写
紧接着函数名的是函数的输入参数，如上函数inlet_x_velocity有两个输入参数：thread和position，thread是一个指针，指向数据类型Thread，position是个整数（见 2.4.3）。
UDF编译和连接之后，函数名就会出现在FLUENT相应的下拉列表内。如上述函数，编译连接之后，就能在相应的边界条件面板内找到一个名为inlet_x_velocity的函数，选定之后就可以使用。 7.2.4.2 udf.h文件中对宏DEFINE_的解释
在udf.h文件中，对附录A的宏作了解释，例如：
#define DEFINE_PROFILE(name，t，I) void name(Thread *t，int i) 通用的宏解释格式为
#define macro replacement-text
在编译前，C预处理器（即 cpp）先进行宏替代。例如 DEFINE_PROFILE(inlet_x_velocity，thread，position) 替代为
void inlet_x_velocity(Thread *thread，int position)
替代后的函数返回实型值或不返回任何值。如上述函数由于是void型的，所以不返回任何值。 7.2.4.3 宏DEFINE
宏DEFINE是用来定义UDFs的，可以分为三类：通用的，离散相的和多相的。从宏DEFINE下划线的后缀，我们可以看出该宏是用来定义哪种类型函数的。
80
如DEFINE_SOURCE定义的函数用来修改输运方程源项，DEFINE_PROPERTY定义的函数用来定义物质的物理性质。通用的宏在2.5详述，离散相和多相的分别在2.6中详述。下面是附录A的简列。
通用类型：
1. DEFINE_ADJUST
2. DEFINE_DIFFUSIVITY 3. DEFINE_HEAT_FLUX 4. DEFINE_INIT
5. DEFINE_ON_DEMAND 6. DEFINE_PROFILE 7. DEFINE_PROPERTY 8. DEFINE_RW_FILE
9. DEFINE_SCAT_PHASE_FUNC 10. DEFINE_SOURCE 11. DEFINE_SR_RATE 12. DEFINE_UDS_FLUX
13. DEFINE_UDS_UNSTEADY 14. DEFINE_VR_RATE 离散相模型：
1. DEFINE_DPM_BODY_FORCE 2. DEFINE_DPM_DRAG 3. DEFINE_DPM_EROSION
4. DEFINE_DPM_INJECTION_INIT 5. DEFINE_DPM_LAW 6. DEFINE_DPM_OUTPUT 7. DEFINE_DPM_PROPERTY
8. DEFINE_DPM_SCALAR_UPDATE 9. DEFINE_DPM_SOURCE 10. DEFINE_DPM_SWITCH 多相模型：
1. DEFINE_DRIFT_DIAMETER 2. DEFINE_SLIP_VELOCITY 7.2.4.2 数据类型的定义
作为对C语言数据类型的补充，FLUENT定义了几种特殊的数据类型，最常用的是：Thread，cell_t，face_t，Node和Domain。
Thread是相应边界或网格区域的结构类型数据；cell_t表示单独一个控制体体积元，用来定义源项或物性；face_t对应于网格面，用来定义入口边界条件等；Node表示相应的网格节点；Domain是一种结构，其中包含所有的threads，cells，faces和nodes。
! Thread，cell_t，face_t，Node和Domain要区分大小写。 7.2.5 通用宏及其定义的函数
宏DEFINE用来定义UDFs，下面是通用宏的具体解释。 7.2.5.1 DEFINE_ADJUST Name DEFINE_ADJUST Arguments domain Arguments Type Domain *domain Return Type void 该函数在每一步迭代开始前，即在求解输运方程前执行。可以用来修改调节流场变量，计算积分或微分等。参数domain在执行时，传递给处理器，通知处理器该函数作用于整个流场的网格区域。
如何激活该函数请参见4.6，具体求解例子见 5.3，5.6和5.7 。 7.2.5.2 DEFINE_DIFFUSIVITY Name DEFINE_DIFFUSIVITY
Arguments c，t， i Arguments Type Return Type cell_t c，Thread *t，real int i 81 该函数定义的是组分扩散系数或者用户自定义标量输运方程的扩散系数，c代表网格，t是指向网格线的指针，i表示第几种组分或第几个用户自定义标量（传递给处理器）。函数返回的是实型数据。例子见5.3。
7.2.5.3 DEFINE_HEAT_FLUX Name DEFINE_HEAT_FLUX Arguments Arguments Type Return Type void f，t，c0，t0， face_t f，Thread *t， cid，cir cell_t c，Thread *t0， real cid[]，real cir[] 该函数定义的是网格与邻近壁面之间扩散和辐射热流量。f表示壁面，t指向壁面线，c0表示邻近壁面的网格，t0指向网格线。函数中需要给出热扩散系数（cid）和辐射系数（cir），才能求出扩散热流量（qid）和辐射热流量（qir）。在计算时，FLUENT按照下面的公式求解：
qid＝cid[0]＋cid[1]×C_T(c0, t0)－cid[2]×F_T(f, t)－cid[3]×pow(F_T(f,t), 4) qir＝cir[0]＋cir[1]×C_T(c0, t0)－cir[2]×F_T(f, t)－cir[3]×pow(F_T(f,t), 4) 该函数无返回值。如何激活函数见4.7，例子见5.6。 7.2.6 DEFINE_INIT Name DEFINE_INIT Arguments domain Arguments Type Domain *domain Return Type void 该函数用于初始化流场变量，它在FLUENT默认的初始化之后执行。作用区域是全场，无返回值。函数的激活见4.5，例子见5.4.1和5.5。
7.2.7 DEFINE_ON_DEMAND Name DEFINE_ON_DEMAND Arguments Arguments Type Return Type void 该函数不是在计算中由FLUENT自动调用，而是根据需要手工调节运行。如何执行见4.12，例子见5.8。 7.2.8 DEFINE_PROFILE Name DEFINE_PROFILE Arguments t，i Arguments Type cell_t c，Thread *t Return Type real 该函数定义边界条件。t指向定义边界条件的网格线，i用来表示边界的位置。函数在执行时，需要循环扫遍所有的边界网格线，值存贮在F_PROFILE(f，t，i)中，无返回值。
选择使用本函数见4.1，例子见5.1.1，5.1.2，5.1.3，5.3，6.1.1和6.1.2。 7.2.9 DEFINE_PROPERTY Name DEFINE_PROPERTY Arguments c，t Arguments Type cell_t c， Thread *t Return Type real 该函数用来定义物质物性参数。c表示网格，t表示网格线，返回实型值。使用见4.3，例子见6.3.1。 7.2.10 DEFINE_RW_FILE Name DEFINE_RW_FILE Arguments fp Arguments Type FILE *fp Return Type void 该函数用于读写case和data文件。fp是指向所读写文件的指针。使用见4.11，例子见2.9.8。 7.2.11 DEFINE_SCAT_PHASE_FUNC Name DEFINE_SCAT _PHASE_FUNC Arguments c，f Arguments Type real c，real *f Return Type real 该函数定义DO（Discrete Ordinate）辐射模型中的散射相函数（radiation scattering phase function）。计算两个变量：从i向到j向 散射的辐射能量分数和前向散射因子（forward scattering factor）。c表示的是i和j向夹角的余弦值，散射的能量分数由函数返回，前向散射因子存贮在指针f所指的变量中。处理器对每种物质，都会调用此函数，分别建立各物质的散射矩阵。
7.2.12 DEFINE_SOURCE Name DEFINE_SOURCE Arguments c，t，dS，i Arguments Type Return Type cell_t c ，Thread *t，real real dS[]，int i 82
该函数定义，除了DO辐射模型之外，输运方程的源项。在计算中，函数需要扫描全场网格。c表示网格，t表示
网格线，dS表示源项对所求输运方程的标量的偏导数，用于对源项的线性化；i标志所定义源项对应于哪个输运方程。使用见 4.2，例子见5.2.1，5.2.2，5.3，6.2.1。
7.2.13 DEFINE_SR_RATE Name DEFINE_SR_RATE Arguments f，t，r，mw，yi，rr Arguments Type Return Type face_t f，Thread *t，void Reaction *r，real *mw，real *yi，real *rr 该函数定义表面化学反应速率。f表示面，t表示面的线，r是结构指针，表示化学反应；mw和yi是个实型指针
数组，mw存贮物质的分子量，yi存贮物质的质量分数，rr设置函数的一个相关参数。函数无返回值，使用见4.8。
7.2.14 DEFINE_UDS_FLUX Name DEFINE_UDS_FLUX Arguments f，t，i Arguments Type Return Type face_t f，Thread *t，int real i 该函数定义用户自定义标量输运方程（user-defined scalar transport equations）的对流通量。f，t分别表示所求通量的面和面的线，i表示第几个输运方程（有处理器传递给本函数）。
7.2.15 DEFINE_UDS_UNSTEADY Name DEFINE_UDS_UNSTEADY Arguments c，t，i， apu，su Arguments Type Return Type cell_t c，Thread *t，void int i，real *apu，real *su 该函数定义用户自定义标量输运方程的非稳态项。c表示网格，t表示网格线，i表示第几个输运方程。在FLUENT中，非稳态项移到RHS中，并以下面的方式离散：
? ????dVunsteady_term???t
?????n?????n?1?
??????V ?t??
??Vn??Vn?1?????
?t?t方程右边第一项为apu，第二项为su。本函数无返回值。 7.2.16 DEFINE_VR_RATE ?Name DEFINE_VR_RATE Arguments Arguments Type Return Type void c，t，r，mw，yi，rr，cell_t c，Thread *t，rr_t Reaction *r，real *mw，real *yi，real *rr，real *rr_t 该函数定义体积化学反应速率。c表示网格，t表示网格线，r表示结构指针，表示化学反应过程，mw指针数组指向存贮物质分子量的变量，yi指向物质的质量分数；rr和rr_t分别设置层流和湍流时函数相关参数。函数无返回值，使用见4.8，例子见6.4.1。
7.2.6 离散相模型宏及其定义的函数
离散模型（DPM）的宏定义的函数与通用宏所定义的函数书写格式是一样的。对于离散相需要强调结构指针p，可以用它得到颗粒的性质和相关信息。下面是具体的宏定义。
7.2.6.1 DEFINE_DPM_BODY_FORCE Name DEFINE_DPM _BODY_FORCE Arguments p，i Argument Type Tracked_Particle *p，int i Return Type real 该函数用于定义除了重力和拉力之外的所有体积力。p为结构指针，i可取0,1,3分别表示三个方向的体积力。函数返回的是加速度。使用见4.4，例子见5.4.2。
7.2.6.2 DEFINE_DPM_DRAG
83
Name Arguments Argument Type Return Type DEFINE_DPM_DRAG Re，p Tracked_Particle *p， real real Re FD?? 224?pDp
函数返回的值是18*CD*Re/24。使用见4.4，例子见5.4.3。 7.2.6.3 DEFINE_DPM_EROSION Name Arguments 该函数定义拉力系数CD，Re为Reynolds数，与颗粒直径和相对于液相速度有关。拉力定义为： 18?CReDArgument Type Return Type DEFINE_DPM_EROSION p，t，f，normal，alpha，Tracked_Particle *p，void Vmag，mdot Thread *t， face_t f，real alpha ，real normal， real Vmag，real mdot 该函数定义颗粒撞击壁面湮灭或产生速率。t为撞击面的线，f为撞击面；数组normal存贮撞击面的单位法向量；alpha中存贮颗粒轨道与撞击面的夹角；Vmag存贮颗粒速度大小，mdot存贮颗粒与壁面撞击率。函数无返回值，颗粒湮灭或产生的计算结果存贮在面变量F_STORAGE_R(f，t，SV_DPMS_EROSION)和F_STORAGE_R(f，t，SV_DPMS_ACCRETION)中。使用见4.4。
7.2.6.4 DEFINE_DPM_INJECTION_INIT Name DEFINE_DPM _INJECTION_INIT Arguments i Argument Type Injection *I Return Type void 该函数用于定义颗粒注入轨道时的物理性质。I是指针，指向颗粒产生时的轨道。对每一次注入，该函数需要在第一步DPM迭代前调用两次，在随后颗粒进入区域前每一次迭代中再调用一次。颗粒的初始化，诸如位置，直径和速度可以通过该函数设定。函数无返回值。
7.2.6.5 DEFINE_DPM_LAW Name Arguments Argument Type Tracked_Particle *p，int ci Return Type void DEFINE_DPM_LAW p，ci 该函数定义液滴和燃烧颗粒的热和质量传输速率。p的意义如前所述，ci表示连续相和离散相是否耦合求解，取1时表示耦合，0时表示不耦合。颗粒的性质随着液滴和颗粒与其周围物质发生传热、传质而改变。函数无返回值。
2.6.6 DEFINE_DPM_OUTPUT Name Arguments Argument Type Return Type DEFINE_DPM_OUTPUT header，void f，p，pt，plane int header， FILE *fp，Tracked_particle *p， Thread *t，Plane *plane 该函数可以得到颗粒通过某一平面（见FLUENT用户手册14.10.6）时的相关变量。header在函数第一次调用时，设为1，以后都为0；fp为文件指针，指向读写相关信息的文件；p为结构指针，t指向颗粒所经过的网格线；plane为 Plane型结构指针（dpm.h），如果颗粒不是穿过一个平面而是仅仅穿过网格表面，值取为NULL。输出信息存贮于指针fp所指向的文件，函数无返回值。
例子见5.4.1。
7.2.6.7 DEFINE_DPM_PROPERTY Name Arguments Argument Type cell_t c， Thread *t，84
Return Type real DEFINE_DPM_PROPERTY c，t，p
Tracked_Particle *p 该函数用于定义离散相物质的物理性质。p为结构指针，c表示网格，t表示网格线。函数返回实型值。 7.2.6.8 DEFINE_DPM_SCALAR_UPDATE Name DEFINE_DPM _SCALAR_UPDATE Arguments c，t，initialize，p Argument Type cell_t c， Thread *t， int initialize，Tracked_particle *p Return Type void 该函数用于更新与颗粒相关的变量或求它们在整个颗粒寿命时间的积分。与颗粒相关的变量可用宏
P_USER_REAL(p，i)取出。c表示颗粒当前所处的网格，t为该网格线。initialize在初始调用本函数时，取为1，其后调用时，取为0。在计算变量对颗粒轨道的积分时，FLUENT就调用本函数。存贮颗粒相关变量的数组大小需要在FLUENT的DPM面版上指定。函数无返回值。
函数的使用见4.4，例子见5.4.1。 7.2.6.9 DEFINE_DPM_SOURCE Name Arguments Argument Type Return Type void DEFINE_DPM_SOURCE c，t，S，strength，p cell_t c，Thread *t，dpms_t *S， real strength， Tracked_Particle *p 该函数用于计算，给定网格中的颗粒在与质量、动量和能量交换项耦合DPM求解前的源项。c表示当前颗粒所在的网格，t为网格线，S为结构指针，指向源项结构dpms_t，其中包含网格的源项。strength表示单位时间内流过的颗粒数目。
p为结构指针。函数求得的源项存贮于S指定的变量中，无返回值。使用见4.4。
7.2.6.10 DEFINE_DPM_SWITCH Name Arguments Argument Type Tracked_Particle *p，int ci Return Type void DEFINE_DPM_SWITCH p，ci 该函数是FLUENT默认的颗粒定律与用户自定义的颗粒定律之间，或不同的默认定律和自定义定律之间的开关函数。p为结构指针，ci为1时，表示连续相与离散相耦合求解，0时表示不耦合求解。
! 参数类型中的Tracked_Particle，dpms_t等是FLUENT为相关模型定义的数据类型。具体含义可以参见相应章节。 7.2.7 多相模型的宏及其定义的函数 7.2.7.1 DEFINE_DRIFT_DIAM Name Arguments Argument Type cell_t c，Thread *t real real DEFINE_DRIFT_DIAM c，t 该函数用于定义代数滑流混合模型（algebraic slip mixture model）颗粒或液滴的直径。c为网格，t为网格线。函数返回颗粒或液滴的直径。使用见4.10。 7.2.7.2 DEFINE_SLIP_VELOCITY Name Arguments Argument Type Domain *domain real void DEFINE_SLIP_VELOCITY domain 该函数用于定义代数滑流混合模型（algebraic slip mixture model）的滑流速度（slip velocity）。该函数作用范围是整个网格区域，无返回值。使用见4.9。 7.2.8 特定线的指针
在很多应用UDF的场合，需要在一条特定的线上进行操作。为满足这种要求，首先，可以从Boundary Conditions 面板得到需要操作的线的ID，然后就可用宏Lookup_Thread将指针指向该条线。
在下面的例子中， C语言函数Print_Thread_Face_Centroids调用FLUENT的宏Lookup_Thread将指针指向特定的线，然后将线上所有面的质心坐标输入文件中。宏DEFINE_ON_DEMAND定义的函数get_coords取出其中的两条线，并打印线上所有面的质心坐标。
#include “udf.h”
extern Domain *domain
85
FILE *fout static void
Print_Thread_Face_Centroids(Domain *domain，int id) {
real FC[2]; face_t f;
Thread *t=Lookup_Thread(domain，id); fprintf(fout, “thread id %d\\n”,id); begin_f_loop(f,，t) {
F_CENTROID(FC，f，t);
fprintf(fout，“f%d %g %g %g\\n”，f，FC[0]，FC[1]，FC[2]); }
end_f_loop(f，t) fprintf(fout，“\\n”); }
DEFINE_ON_DEMAND(get_coords) {
fout = fopen(“faces.out”，“w”);
Printf_Thread_Face_Centroids(domain，2); Printf_Thread_Face_Centroids(domain，4); fclose(fout); }
7.2.9 函数体 7.2.9.1 介绍
用户自定义函数体部分包含在紧跟着宏DEFINE_定义的大括弧内，如下例。这和标准C语言函数体的定义是相同的。
DEFINE_PROPERTY(cell_viscosity，cell，thread) {
real mu_lam;
real temp = C_T(cell,thread); if (temp > 288.) mu_lam = 5.5e-3; else if (temp > 286.)
mu_lam = 143.2135 – 0.49725 * temp; else
mu_lam = 1.; return mu_lam; }
7.2.9.2 Interpreted UDFs的限制性
Interpreted型书写函数体时并不能完全应用C语言函数，这在前面有过论述。此外，数量的单位制必须采用国际单位制。
7.2.9.3 函数的功能
UDFs执行五种功能：
1． 返回变量值； 2． 调节参数；
3． 返回变量值并且调节参数；
4． 调节FLUENT的变量（不以参数形式传递）； 5． 向case或data文件读写信息。
宏定义的函数返回类型如果不是void型，就返回实型值。下面的例子分别说明不同功能的函数。 7.2.9.4 返回变量值
86
下面的UDF计算与温度有关的粘性系数，并返回该值。 #include “udf.h”
DEFINE_PROPERTY(cell_viscosity，cell，thread)/*定义物性的宏*/ {
real mu_lam;
real temp = C_T(cell,thread);/*变量temp存放网格的温度*/ if (temp > 288.) mu_lam = 5.5e-3; else if (temp > 286.)
mu_lam = 143.2135 – 0.49725 * temp; else
mu_lam = 1.;
return mu_lam;/*变量mu_lam存放粘性系数值，函数返回该值*/ }
7.2.9.5 调节参数
下面的UDF给出简单二元气相系统的体积反应速率。 #include “udf.h”
#define K1 2.0e-2 #define K2 5.
DEFINE_VR_RATE(user_rate，cell，thread，r，mole_weight，species_mf，rate，rr_t) {
real s1 = species_mf[0]; real mw1 = mole_weight[0];
if (FLUID_THREAD_P(thread) && THREAD_VAR(thread) .fluid. porous) *rate = K1*s1/pow((1.+K2*s1),2.)/mw1; else
*rate = 0.; }
函数名为 user_rate，函数体中用if语句判断是否处于多孔介质区，仅在多孔介质区有化学反应。 7.2.9.6 返回变量值并调节参数
下面的UDF定义的是swirl-velocity的源项。 #include “udf.h”
#define OMEGA 50 /* rotational speed of swirler */
#define WEIGHT 1.e20 /* weighting coefficients in linearized equation */ DEFINE_SOURCE(user_swirl，cell，thread，dS，eqn) {
real w_vel，x[ND_ND]，y，source; CENTROID(x，cell，thread); y = x[1];
w_vel = y*OMEGA; /* linear w_velocity at the cell */ source = WEIGHT*(w_vel – C_WSWIRL(cell，thread)); dS[eqn] = －WEIGHT; return source; }
函数名为user_swirl，函数计算了变量source并且返回其值。函数的各项参数的意义参见2.5.10。 7.2.9.7 调节FLUENT变量
下面的函数调节存贮于内存的FLUENT变量，函数定义了x方向速度的边界条件。 #include “udf.h”
DEFINE_PROFILE(inlet_x_velocity，thread，position) {
87
real x[ND_ND]; real y; face_t f;
begin_f_loop(f，thread) {
F_CENTROID(x，f，thread); y = x[1];
F_PROFILE(f，thread，position) = 20. －y*y/(.0745*.0745)*20; }
end_f_loop(f，thread) }
函数的参数position是个数字标签，标记每一步循环（loop）中被设置的变量。函数调节的FLUENT变量F_PROFILE。
7.2.9.8 读写data或case文件
下面的函数介绍了如何读写静态变量kount，如何计算静态变量请参见4.6。 #include “udf.h”
int kount = 0; /*定义静态变量kount*/ DEFINE_ADJUST(demo_calc，domain) {
kount ++;
printf(“kount = %d\\n”，kount)； }
DEFINE_RW_FILE(writer，fp) {
printf(“Writing UDF data to data file…\\n”);
fprintf(fp，”%d”，kount); /*将kount写入data文件中*/ }
DEFINE_RW_FILE(writer，fp) {
printf(“Reading UDF data from data file…\\n”);
fscanf(fp，“%d”，&kount); /*从数据文件中读取kount值*/ }
上面有三个函数。如果迭代10次，则kount值为10，然后将当前值10存贮到数据文件中，如果下次将kount值读入FLUENT继续运算，则kount将在10的基础上增加。我们可以存贮任意多的静态变量，不过读写顺序必须一致。 7.2.10 解法器函数（Solver Functions） 7.2.10.1 概述
在很多情况下，UDF需要得到FLUENT解法器中的数据。例如：
1． 所求解的变量及其导数（例如，速度，温度等）；
2． 网格和面几何性质（例如，面面积，网格体积，网格质心坐标等）； 3． 物质的物理性质（例如，密度，粘性系数，导热系数等）。
! 我们可以取出比热，但是不能修改。
我们可以利用下一节所列FLUENT提供的解法器函数，得到解法器中的数据。这里所说的函数是从广义上讲的，因为其中包括函数和宏，只有在源文件appropriate.h中定义的才是真正的函数。
! 如果使用的是Interpreted 型的UDF，则只能使用这些FLUENT提供的解法器函数。
解法器函数可以与C函数一起在函数体中混合使用。为方便起见，一些最常用的C函数列在附录B中。 下面章节列出的函数中包括它们的参数，参数类型和返回值，还有对该函数说明的源文件。例如 C_CENTROID(x，c，t)
有三个参数：x，c和t，其中c和t为输入参数，x为输出参数，输出网格的坐标值。 7.2.10.2 辅助几何关系
88
Name(Arguments) C_NNODES(c，t) C_NFACES(c，t) F_NNODES(f，t) Name(Arguments) Argument Type cell_t c, Thread *t cell_t c, Thread *t face_t f, Thread *t Returns 网格节点数 网格面数 面节点数 Source mem.h mem.h mem.h Returns x(网格坐标) x(面坐标) A(面矢量) 面矢量A大小 2D或3D网格体积；对称体网格体积/2π Source metric.h metric.h metric.h metric.h metric.h 7.2.10.3 网格坐标与面积 Argument Type real x[ND_ND]，cell_t c, Thread *t real x[ND_ND], face_t f, Thread *t A[ND_ND]，face_t f，Thread *t A[ND_ND] cell_t c, Thread *t C_CENTROID(x，c，t) F_CENTROID(x，f，t) F_AREA(A，f，t) NV_MAG(A) C_VOLUME(c ，t) 7.2.10.4 节点坐标与节点（网格）速度
列表中的节点速度与移动网格模拟有关。 Name(Arguments) NODE_X[node] NODE_Y[node] NODE_Z[node] NODE_GX[node] NODE_GY[node] NODE_GZ[node] Argument Type Node *node Node *node Node *node Node *node Node *node Node *node Returns 节点的x坐标 节点的y坐标 节点的z坐标 节点的x向速度 节点的y向速度 节点的z向速度 Source metric.h metric.h metric.h mem.h mem.h mem.h 7.2.10.5 面变量
下面列表中的函数只能在segregated solver中获取，耦合计算时不能引用的。 Name(Arguments) F_P(f，t) F_U(f，t) F_V(f，t) F_W(f，t) F_T(f，t) F_H(f，t) F_K(f，t) F_D(f，t) F_YI(f，t，i) F_UDSI(f，t，i) F_UDMI(f，t，i) Argument Type face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t face_t f, Thread *t， int i face_t f, Thread *t， int i face_t f, Thread *t， int i Returns 压力 u速度 v速度 w速度 温度 焓 湍流动能 湍流能量耗散系数 组分质量分数 用户自定义标量（i表示第几个方程） 用户自定义内存变量（i表示第几个） Source mem.h(all) 7.2.10.6 网格变量
下面列表中的是网格变量，不像面变量，网格变量在耦合与非耦合计算中都能获取（available in both the segregated and the coupled solvers）。下面三个列表中所列出的分别是：计算变量、导数和物性参数。 Name(Arguments) C_P(c，t) C_U(c，t) C_V(c，t)
Argument Type 参数类型都定义为： cell_t c, Returns 压力 u速度 v速度 89
Source mem.h C_W(c，t) C_T(c，t) C_H(c，t) C_YI(c，t) C_UDSI(c，t，i) C_UDMI(c，t，i) C_K(c，t) C_D(c，t) C_RUU(c，t) C_RVV(c，t) C_RWW(c，t) C_RUV(c，t) C_RVW(c，t) C_RUW(c，t) C_FMEAN(c，t) C_FMEAN2(c，t) C_FVAR(c，t) C_FVAR2(c，t) C_PREMIXC(c，t) C_LAM_FLAME _SPEED(c，t) C_CRITICAL_STRAIN _RATE(c，t) C_POLLUT(c，t，i) C_VOF(c，t，0) C_VOF(c，t，1) Name(Arguments) C_DUDX(c，t) C_DUDY(c，t) C_DUDZ(c，t) C_DVDX(c，t) C_DVDY(c，t) C_DVDZ(c，t) C_DWDX(c，t) C_DWDY(c，t) C_DWDZ(c，t) C_DP(c，t)[i] C_D_DENSITY(c，t)[i] Name(Arguments) C_R(c，t) C_MU_L(c，t) C_MU_T(c，t) C_MU_EFF(c，t) C_K_L(c，t) C_K_T(c，t) C_K_EFF(c，t) C_CP(c，t) C_RGAS(c，t) C_DIFF_L(c，t，i，j) Thread *t, int i, int j w速度 温度 焓 组分质量分数 用户自定义标量 用户自定义内存变量 湍流动能 湍流能量耗散系数 uu雷诺应力 vv雷诺应力 ww雷诺应力 uv雷诺应力 vw雷诺应力 uw雷诺应力 第一平均混合物分数 第二平均混合物分数 第一混合物分数偏差 第二混合物分数偏差 反应进程变量 层流火焰速度 临界应变率 污染物组分 第一相体积分数 第二相体积分数 Returns 各向速度导数 压力梯度(i表示方向) 密度梯度(i表示方向） Returns 密度 层流粘性系数 湍流粘性系数 有效粘性系数 导热系数 湍流导热系数 有效导热系数 比热 气体常数 层流组分扩散系数 90 sg_mem.h Argument Type Source mem.h Argument Type Source mem.h C_DIFF_EFF(c，t，i) C_ABS_COEFF(c，t) C_SCAT_COEFF(c，t) 有效组分扩散系数 吸收系数 散射系数 7.2.10.7 循环宏
如果要实现扫描全场的网格就需要使用循环宏。下面给出两种类型的循环宏。一种以begin开始，end结束，用来扫描线上的所有网格和面；另一种用来扫描所有的线。
cell_t c; face_t f; Thread *t; Domain *d;
begin_c_loop(c, t) { }
end_c_loop(c, t) /*循环遍历线上的所有网格*/ begin_f_loop(f, t) { }
end_f_loop(f, t) /*循环遍历线上的所有面*/ thread_loop_c(t, d) {
} /*遍历网格线*/ thread_loop_f(t, d) {
} /*遍历面上的线*/ 7.2.10.8 数据的可得性
在书写UDF的时候，要保证出现在函数中的数据是能够从解法器上获得的。为了检验数据的这种可得性，我们可以使用函数Data_Valid_p，如果数据正确该函数返回值为1，否则返回0。
Int Data_Valid_P(); /*取自文件case.h*/
当我们读取case文件时，就会装载相应的UDF。如果此时函数用到一个未初始化的变量，诸如内部网格速度，计算就会发生错误。为了避免这种类型的错误发生，解法器会执行一条if else 语句。如果数据可得，就按照正常情况执行下去；否则，就停止运算。一旦流场初始化之后，函数会被重新激活，然后读取正确的数据运行。 7.2.10.9 设置面变量值的函数
宏F_PROFILE设置的是面上的变量值，该函数在设置边界面时使用。 face_t f; Thread *t;
F_PROFILE(f,，t,，nvar) /* from mem.h */
函数中的参数nvar不需用户设定，而是由FLUENT解法器传递给UDF。nvar是特定边界的数字标签。例如，与总压和总温相关的入口边界，nvar取值为0和1；与速度三个方向标量和静态温度相关的边界，nvar取值为0，1，2和3。
我们在定义边界条件面板定义边界条件时，解法器就设定nvar的值。 ! 整型变量nvar是不能改变的，只能由FLUENT解法器传递。 7.2.11 DPM宏
下面列出的离散相模型的宏定义在源文件dpm.h中，变量p是指向结构Tracked_Particle的指针（Tracked_Particle *p）。 7.2.11.1 颗粒变量的宏
1． 当前位置颗粒变量
P_DIAM(p) 颗粒直径
P_VEL(p)[i] 颗粒速度，i＝0，1，2 P_T(p) 颗粒温度 P_RHO(p) 颗粒密度 P_MASS(p) 颗粒质量
91
P_TIME(p) 当前颗粒时间 P_DT(p) 颗粒时间步长
P_LF(p) 颗粒液相分数（只适用于湿燃烧颗粒） P_VFF(p) 颗粒挥发性馏分（只适用于燃烧颗粒） 2． 进入当前网格颗粒变量
P_DIAM0(p) 颗粒直径
P_VEL0(p)[i] 颗粒速度，i＝0，1，2 P_T0(p) 颗粒温度 P_RHO0(p) 颗粒密度 P_MASS0(p) 颗粒质量 P_TIME0(p) 颗粒时间
P_LF0(p) 颗粒液相分数（只适用于湿燃烧颗粒） 3． 刚注射入区域的颗粒变量
P_INIT_DIAM(p) 颗粒直径 P_INIT_MASS(p) 颗粒质量 P_INIT_RHO(p) 颗粒密度 P_INIT_TEMP(p) 颗粒温度
P_INIT_LF(p) 颗粒液相分数（只适用于湿燃烧颗粒）
P_EVAP_SPECIES_INDEX(p) 混合物蒸发组分数 （evaporating species index in mixture） P_DEVOL_SPECIES_INDEX(p) 混合物液化组分数 (devolatizing species index in mixture) P_OXID_SPECIES_INDEX(p) 混合物氧化物组分数 (oxidizing species index in mixture)
P_PROD_SPECIES_INDEX(p) 混合物燃烧产物分数 (combustion products species index in mixture) P_CURRENT_LAW(p) current particle law index P_NEXT_LAW(p) next particle law index
P_USER_REAL(p，i) 用户自定义变量（i表示第几个）
7.2.11.2 颗粒物性的宏
P_MATERIAL(p) 颗粒指针 DPM_SWELLING_COEFF(p) 颗粒膨胀系数
DPM_EMISSIVITY(p) 辐射模型颗粒发射率 DPM_SCATT_FACTOR(p) 辐射模型颗粒散射因子 DPM_EVAPORATION_TEMPERATURE(p) 颗粒物质蒸发温度 DPM_BOILING_TEMPERATURE(p) 颗粒物质沸腾温度 DPM_LATENT_HEAT(p) 颗粒物质潜热 DPM_HEAT_OF_PYROLYSIS(p) 颗粒物质分解热 DPM_HEAT_OF_REACTION(p) 颗粒物质反应热 DPM_VOLATILE_FRACTION(p) 颗粒物质挥发性馏分 DPM_CHAR_FRACTION(p) 颗粒物质炭化分数
DPM_SPECIFIC_HEAT(p，t) 颗粒物质在温度t时的比热
第三节 编译连接UDFs
函数的编译与连接取决于使用的UDF是Interpreted型，还是Compiled型。下面3.3和3.4节也分别介绍了Interpreted型和Compiled型UDF的编译和连接。 7.3.1 概述
我们在第一章就介绍过Interpreted型和Compiled型UDF的区别以及各自的特点，在此不再累述。为了显示区别，FLUENT中两种UDF的设置面板是不同的。Interpreted型面板有一个Compile按钮，点击时就会对源文件进行编译；Compiled型面板有一个open按钮，点击时会打开事先编译好的机器源代码库，然后连接运行。
UDF是用宏DEFINE_定义的，对宏的解释在udf.h文件中。udf.h文件存放的位置为path/Fluent.Inc/fluent5.x/src/udf.h，path为FLUENT安装路径。在编译的时候，编译器会自动到src目录下搜寻udf.h文件。注意，在更新src目录时不能移去udf.h文件，因为该文件在近期是无法更新的。如果系统搜寻不到udf.h文件，编译连接工作将无法进行。 7.3.2 Interpreted型UDFs
本节介绍了如何编译Interpreted型UDFs。编译之后，case文件会记录编译过程，FLUENT在读取case文件时会
92
自动加载UDFs。
7.3.2.1 编译Interpreted型UDFs
一般，我们按照下面的步骤编译UDF： 1． 将UDF的C源程序拷贝到当前工作目录，如果不在当前目录，则编译时需要指定该文件的具体路径。如
果使用并行运算，按照3.2.2设置；
2． 运行FLUENT程序； 3． 读取（或建立）case文件； 4． 在Interpreted UDFs面板编译UDF（例如，vprofile.c）。
Define?User-Defined?Functions?Interpreted…
(a) (b)
图3.2.1 Interpreted UDFs面板
在Source File Name中输入C源文件名（如vprofile.c）。如果源文件不在当前工作目录下，则需要输入完整路径。
在CPP Command Name中输入C预处理程序名。如果安装了/contrib组件，则会出现默认的预处理程序名：cpp。对于Windows NT用户，标准安装过程会自动安装默认的预处理程序。当然，我们还可以选择系统中其它的预处理程序。
Stack Size默认值取10000，如果我们处理的问题变量很多，10000个栈会溢出，则需要加大该值。 点击Compile键进行编译。如果选中Display Assembly Listing按钮，编译时会出现下面的信息提示：
(c) (d)
93
(e) 点击Close关闭。
! 如果我们需要编译的Interpreted UDF不止一个，可以将UDFs写入单个C源文件，如all.c，然后以同样的方式进行编译连接就可以了。
7.3.2.2 Windows NT系统网络并行运算的目录结构
如果我们使用的是Windows NT系统的并行版FLUENT，就需要以特定的方式组织文件目录。3.2.1中的第一步必须以下面的步骤代替：
1． 在Fluent.Inc目录下建立一个名为udf的可写子目录。 2． 在目录udf下建立一个目录（例如，Fluent.Inc\%udf\\myudf），将C源文件拷入其中，如myudf目录下。如果多
位用户使用，可以再建立其它目录（如Fluent.Inc\%udf\\abcudf 或xyz）。
! 由于C源文件不在当前工作目录下，编译时需要输入文件路径。如：
\\\\
3． 如果已有case文件，应保证其在当前工作目录下。 7.3.2.3 调试Interpreted UDFs
如果编译有错，在相应窗口中会出现错误信息。如果因为屏幕翻动无法看清，则应该关闭Display Assembly Listing选项。
7.3.2.4 Interpreted UDFs常见编译错误
C源文件如果不是在当前工作目录下，而在编译时没有指明具体路径，就会出现下面的错误：
如果将case文件拷贝到了其他目录，而没有将C源文件一起拷过去，会出现同样的错误。 7.3.3 Compiled型UDFs
一般，我们按照下面的步骤编译连接Compiled型UDFs： 1． 在工作目录下建立相应的目录（见3.3.1）； 2． 编译UDFs并建立共享库（见3.3.2）； 3． 运行FLUENT；
4． 读取（或建立）case文件（case文件如果事先以存在，应保证在当前工作目录下）； 5． 连接共享库到FLUENT执行（见3.3.3）。连接成功后，出现下面信息：
7.3.3.1 建立目录结构
UNIX和Windows NT系统所建目录结构是不一样的，下面分别进行说明。 UNIX系统
UNIX系统编译Compiled UDFs需要两个文件：Makefile和makefile。makefile文件是可修改的，C源文件名就在其中设定。两文件所处的相关目录为：
path/Fluent.Inc/fluent5.x/src/Makefile path/Fluent.Inc/fluent5.x/src/makefile
path为FLUENT安装目录，5.x表示FLUENT版本。
下面的步骤概括介绍了目录结构的建立过程，从图3.3.1可以看得更清楚。图3.3.1中的目录结构适用于两种版本：单精度2D运算和单精度2D并行运算。注意不需要向版本目录（2d，2d_host，2d_node）中拷入任何文件，图中所示文件在编译时会自动生成。
94
图3.3.1 UNIX目录结构说明
1． 在当前工作目录下建立一个目录用来存放库（如libudf）。
2． 拷贝Makefile.udf文件到1中建立的目录下，并重新命名为Makefile。
3． 在1中建立的目录（libudf）下，建立一个用于存放C源文件的目录，并命名为src。 4． 将C源文件（如udfexample.c）拷入src目录下。
5． 把makefile.udf文件拷入src目录下，并命名为makefile。
6． 启动FLUENT，确认FLUENT版本结构。如果是irix6.5，需要修改makefile文件（见3.3.2）。 7． 按照FLUENT版本结构建立相应的目录结构（如ultra/2d或ultra/3d等）。可能的版本有：
(a) 单精度2d或3d：2d或3d (b) 双精度2d或3d：2ddp或3ddp
(c) 单精度并行版2d或3d：2d_node和2d_host，3d_node和3d_host
(d) 双精度并行版2d或3d：2ddp_node和2ddp_host，3ddp_node和3ddp_host Windows NT系统
对于Windows NT系统需要两个文件：makefile_nt.udf和user_nt.udf。C源文件需要修改文件user_nt.udf指定。下面的步骤概括介绍了目录结构的建立过程，可参考图3.3.2。图3.3.2也只适用于两种版本：单精度2D运算和单精度2D并行运算。
图3.3.2 Windows NT目录结构说明
1． 在当前工作目录下建立一个目录（如libudf）。
2． 在新建目录下建立名为src目录，用于存放C源程序文件。 3． 将C源程序（如udfexample）拷贝入src目录。
4． 根据机器结构建立相应目录：Windows NT系统建立ntx86目录，DEC Alpha系统建立ntalpha目录。 5． 在ntx86目录或ntalpah目录下建立版本信息目录（如ntx86\\2d）。可能的FLUENT版本为：
(a) 单精度2d或3d：2d或3d (b) 双精度2d或3d：2ddp或3ddp
(c) 单精度并行版2d或3d：2d_node和2d_host，3d_node和3d_host
(d) 双精度并行版2d或3d：2ddp_node和2ddp_host，3ddp_node和3ddp_host 6． 将makefile_nt.udf和user_nt.udf拷入版本信息目录（如2d）。对于并行版本需要将这两个文件分别拷入host
95
和node目录（如2d_node和2d_host）。两文件在FLUENT中的位置： path/Fluent.Inc/fluent5.x/src/makefile_nt.udf
path/Fluent.Inc/fluent5.x/src/user_nt.udf
path为FLUENT安装目录。注意FLUENT版本更新升级后，需要把新版本的makefile_nt.udf和user_nt.udf拷入覆盖原来的文件，这和UNIX系统是不同的。 7． 把文件makefile_nt.udf重新命名为makefile。 7.3.3.2 编译建立共享库
对于UNIX和Windows NT系统编译和建立共享库的方式是不同的，下面分别叙述： UNIX系统
在建立相应的目录结构和拷入源文件之后，就可以编译生成共享库文件。 1． 编辑src目录下的makefile文件，并设置参数：
SOURCES = 要进行编译的C源文件名 FLUENT_INC = FLUENT安装路径 下面是个makefile文件的例子：
2． 如果版本结构是irix6.5，则需要进一步修改makefile文件。
(a) 在文件中找到下面部分：
CFLAGS_IRIX6R10 = －KPIC －ansi －fullwarn －0 －n32 (b) 将－ansi改为－xansi：
CFLAGS_IRIX6R10 = －KPIC －xansi －fullwarn －0 －n32 3． 在库目录（如libudf）下，输入下面的语句，执行Makefile。
Make “FLUENT_ARCH=ultra” 系统显示信息：
上述信息表明编译udfexample.c成功，同时系统会建立共享库libudf.so。
96
当然，我们可以同时编译多个C源程序文件。
Windows NT系统
在建立相应目录及拷入C源文件之后，Windows NT系统按照下面的步骤编译和建立共享库。 1． 编辑user_nt.udf文件，设置下面的参数：
SOURCES＝要进行编译的C原文件名。在文件名之前，需要前缀$(SRC) （如$(SRC)udfexample.c），对于多个源文件，文件之间要以空格间隔（如 $(SRC) udfexample1.c $(SRC) udfexample2.c）。
VERSION＝所使用解法器版本，如2d，3d，2ddp，3ddp等。 PARALLEL_NODE＝并行运算通讯库，可能的情况有： －none：单机版
－smpi：并行使用共享内存（对多处理器机器）
－vmpi：使用MPI软件实现并行使用共享内存和网络 －net： 使用RSHD软件并行使用网络通讯
! 并行运算注意修改所有的user_nt.udf文件
下面是user_nt.udf参数设置的例子：
2． 在MS－DOS操作符下，进入到所建立目录（如\\libudf\\ntx86\\2d\\）下，执行命令nmake，如果编译出现问题使
用nmake clean 清除。
7.3.3.3 连接共享库到FLUENT
按上一节编译并建立共享库之后，就可以在FLUENT中，将共享库连接使用。具体步骤如下： 1． 运行FLUENT程序
2． 读取（或建立）case文件 3． 连接共享库到FLUENT
Define?User-Defined?Functions?Compiled…
图3.3.3 Compiled UDFs面板
(a) 在Library Name中输入存储库文件的目录名（如libudf），如果不在当前工作目录下，需要输入路径。 (b) 点击Open连接共享库到FLUENT。解法器会自动找到相应的版本的共享库，一旦连接成功，将会在case
文件中保存，下次读取case文件时会自动加载共享库。
7.3.3.4 编译连接Compiled UDFs时的常见错误
当共享库不在当前工作目录下的相应目录内时，或者取出case文件而与该case文件相关联的共享库已经移动时，会显示如下错误：
97
如果在新版本下读取旧版本的case文件，如在5.0.2中读取5.0.1版本的case文件，则会显示下面的错误：
第四节 在FLUENT模型中使用UDFS
在编译连接成功之后，我们就可以在FLUENT模型中使用UDFs。下面分别说明如何在FLUENT各模型中使用UDFs。 7.4.1 边界条件
编译连接成功边界条件UDF之后，就可以在FLUENT中选择使用。例如，定义速度入口边界条件的UDF（inlet_x_velocity）可以在Velocity Inlet面板（如下所示）中找到，选择之后点击OK按钮就可以在FLUENT中使用。
边界条件UDFs由宏DEFINE_PROFILE定义，参见2.5.6。 7.4.2 源项
源项UDFs在Fluid面板中选择，如名为cell_x_source的源项UDF。
源项UDFs由宏DEFINE_SOURCE定义，参见2.5.10。 7.4.3 物质性质
98
定义物质性质的UDFs在Material面板中选择，如定义物质粘性系数的UDF cell_viscosity。
Define?Materials…
点击EDIT…
! 如果定义的是密度的UDF，密度变化越大，计算收敛性就越差。所以FLUENT限制计算过程密度变化过大的模型的使用，如可压缩流动，如果密度随压力变化太大，而压力又在计算过程中变化很大，这就不适合用FLUENT求解。
定义物质性质的UDFs由宏DEFINE_PROPERTY定义，参见2.5.7。 7.4.4 离散相模型
离散相模型的UDFs在Discrete Phase Model面板相应的下拉列表中选择。 Define?Models?Discrete Phase…
99
如果使用“变量更新”（scalar update）函数，更新与颗粒轨道有关的标量值，可以从Scalar Update下拉列表中选择。变量数目在Number of Scalars中指定。 离散相宏定义参见2.6。 7.4.5 计算初始化
计算初始化UDF在User-Define Function面板的Initialization Function中设置。 Define?User-Defined?Function Hooks…
100
计算初始化程序在FLUENT默认初始化之后调用，用于初始化流场变量。参见2.5.4。 7.4.6 更新计算变量
更新计算变量函数在User-Define Function Hooks面板的Adjust Function设定，见4.5的图示。更新计算变量函数在迭代过程中每一步都调用，用来更新全场的变量。参见2.5.1。 7.4.7 壁面热流密度
壁面热流密度函数在User-Define Function Hooks面板的Wall Heat Flux Function中设置，见2.5.3。通过壁面热流密度函数，我们可以设定壁面网格和邻近壁面网格的热流量。例如，可以改变热交换系数和温度的壁面定律（参见相应章节）。
7.4.8 反应速率
体积反应速率和表面反应速率函数分别在User-Define Function Hooks面板的Volume Reaction Rate Function和Surface Reaction Rate Function中设置。参见2.5.11。 7.4.9 Algebraic Slip Mixture 模型的Slip Velocity
选用Algebraic Slip Mixture模型之后，滑流速度(Slip Velocity)选项就可以在User-Define Function Hooks面板中进行选择。参见3.2和3.3。 7.4.10 读写case和data文件
读写case和data文件函数分别在User-Define Function Hooks面板的Read Case Function，Write Case Function，Read Data Function和Write Data Function中设定。四个函数由宏DEFINE_RW_FILE定义，参见2.5.8。 7.4.11 颗粒或液滴直径
定义颗粒或液滴直径的函数在Multiphase Model面板中设置。 Define?models?Multiphase…
7.4.12 可执行UDF
可执行UDF在Execute UDF On Demand面板中设置。 Define?User-Defined?Execute On Demand…
可执行UDF由宏DEFINE_ON_DEMAND定义。参见2.5.5。 7.4.13 内存存贮变量
在UDF中，我们可以将网格变量和参数存贮到内存变量中，以便其后取出使用。内存变量个数在User-Defines Memory面板中指定。
Define?User-Defined?Memory …
101
定义了内存变量之后，可以通过C_UDMI(c，t，i)和F_UDMI(f，t，i)引用。i取0,1,2…，表示第几个内存变量。每一个网格最多可以对应500个这样的内存变量，以在UDF中使用。例子参见5.8。
第五节 UDFs实例
本章详细分析了UDFs的各种用法的实例，具体的应用将在第六章中分析。本章5.3，5.4和5.6的例子是Compiled型的，其它都是Interpreted型的。下面各节，给出了UDF所描述的公式及其C源程序，其中所用到的C语言语法请参考相关书籍，FLUENT语法和函数请参考2.4.3和2.10，或者可以直接参照相应的*.h文件。 7.5.1 边界条件
本节给出了三个边界条件UDFs的例子，分别描述： 1． 涡轮叶片抛物型入口压力分布。 2． 完全发展的湍流流动入口条件。 3． 正弦壁面温度分布。
7.5.1.1 涡轮叶片抛物型入口压力分布
压力分布方程为：
2 y??p(y)?1.1?105－0.1?105?? 0.0745??
下面是根据该方程书写的UDF的C源程序。
/*************************************************************************/ /* profile.c */
/* UDF for specifying a steady-state pressure profile boundary condition */
/*************************************************************************/ #include \
DEFINE_PROFILE(pressure_profile, thread, position) {real x[ND_ND]; /* this will hold the position vector */ real y; face_t f;
begin_f_loop(f, thread)
{F_CENTROID(x,f,thread); y = x[1];
F_PROFILE(f, thread, position) = 1.1e5 - y*y/(.0745*.0745)*0.1e5;} end_f_loop(f, thread)}
函数名为pressure_profile，数组x通过函数F_CENTROID得到面f的质心，loop循环扫描整个定义的边界，计算每个网格面边界条件的值。参数position是边界网格面的标志，可以把求得的值与网格面一一对应。边界面上的值保存入函数F_PROFILE中，FLUENT在计算过程中处理边界时，会取出所存贮的值。x[1]对应纵坐标y，同理，x[0]和x[2]分别对应坐标x,y的值。
7.5.1.2 完全发展的湍流流动入口条件
本例使用三个UDFs分别定义：x向速度，湍流动能和耗散率，近似模拟完全发展的湍流流动。 速度采用1/7定律： 1/7y?? vx?vx,free?????假设湍流动能从壁面附近值直线变化到自由流的值，即从
ur2knw? C?
102
2线性变为 k inf ? 0 . 002 u 。 free
3/43/2耗散率为 ? ? C ? l为 ? y和 0 .085 两者之中的最小值，?? 为Von Karman常数。 k ，混和长度l摩擦速度和壁面剪切力取为：
ur??w/? ?w?f?ufree/2 摩擦因子采用Blasius方程：
2?ufree? f?0.045??v??????1/4
下面是根据上述方程编写的UDFs的C语言源程序。
/************************************************************************/ /* UDF for specifying fully-developed profile boundary conditions */
/************************************************************************/ #include \
#define YMIN 0.0 /* constants */ #define YMAX 0.4064 #define UMEAN 1.0 #define B 1./7.
#define DELOVRH 0.5 #define VISC 1.7894e-05 #define CMU 0.09 #define VKC 0.41
/* profile for x-velocity */
DEFINE_PROFILE(x_velocity, thread, position)
{real y, del, h, x[ND_ND], ufree; /* variable declarations */ face_t f;
h = YMAX - YMIN; del = DELOVRH*h; ufree = UMEAN*(B+1.); begin_f_loop(f, thread) {
F_CENTROID(x,f,thread); y = x[1]; if (y <= del)
F_PROFILE(f,thread,position) = ufree*pow(y/del,B); else
F_PROFILE(f,thread,position) = ufree*pow((h-y)/del,B); }
end_f_loop(f, thread) }
/* profile for kinetic energy */
DEFINE_PROFILE(k_profile, thread, position) {
real y, del, h, ufree, x[ND_ND]; real ff, utau, knw, kinf; face_t f;
h = YMAX - YMIN;
103
del = DELOVRH*h; ufree = UMEAN*(B+1.);
ff = 0.045/pow(ufree*del/VISC,0.25); utau=sqrt(ff*pow(ufree,2.)/2.0); knw=pow(utau,2.)/sqrt(CMU); kinf=0.002*pow(ufree,2.); begin_f_loop(f, thread) {
F_CENTROID(x,f,thread); y=x[1]; if (y <= del)
F_PROFILE(f,thread,position)=knw+y/del*(kinf-knw); else
F_PROFILE(f,thread,position)=knw+(h-y)/del*(kinf-knw); }
end_f_loop(f, thread)
}
/* profile for dissipation rate /*
DEFINE_PROFILE(dissip_profile, thread, position) {
real y, x[ND_ND], del, h, ufree; real ff, utau, knw, kinf; real mix, kay; face_t f;
h = YMAX - YMIN; del = DELOVRH*h; ufree = UMEAN*(B+1.);
ff = 0.045/pow(ufree*del/VISC,0.25); utau=sqrt(ff*pow(ufree,2.)/2.0); knw=pow(utau,2.)/sqrt(CMU); kinf=0.002*pow(ufree,2.); begin_f_loop(f, thread) {
F_CENTROID(x,f,thread); y=x[1]; if (y <= del)
kay=knw+y/del*(kinf-knw); else
kay=knw+(h-y)/del*(kinf-knw); if (VKC*y < 0.085*del) mix = VKC*y; else
mix = 0.085*del;
F_PROFILE(f,thread,position)=pow(CMU,0.75)*pow(kay,1.5)/mix; }
end_f_loop(f,thread) }
7.5.1.3 正弦壁面温度分布
壁面温度随着x呈正弦分布，靠近壁最左端温度为300K，中间为400K，最右端时又降为300K，分布式为：
104
T(x)?300?100sin??x???
?0.005? UDF的C语言源代码如下。
/*************************************************************************/ /* UDF for specifying a sinusoidal temperature profile boundary condition */
/*************************************************************************/ #include \
#define PI 3.141592654
DEFINE_PROFILE(temperature_profile, thread, position) {
real r[ND_ND]; /* this will hold the position vector */ real x; face_t f;
begin_f_loop(f, thread) {
F_CENTROID(r,f,thread); x = r[0];
F_PROFILE(f, thread, position) = 300.+100.*sin(PI*x/0.005); }
end_f_loop(f, thread) }
7.5.2 源项
本节分析源项UDFs的两个例子： 1． 与空间位置相关的多孔介质区 2． 流体区域涡产生器 7.5.2.1 多孔介质区的源项
多孔介质区域的源项方程式为：source??0.5C2?yvxvx
?eS??Avxvx ，A?0.5C2?y 假设 sourc于是，有
dSd?vx? ??Avx?Avxdvxdvx所以，UDF的返回源项值为： source??Avxvx 源项对输运方程所求变量vx的导数为：
dS??2Avx dvx源项UDF必须给出源项对所求变量的直接导数，源项的导数用于对源项的线性化（见下一节）。下面是根据上述方程编写的UDF的C源程序。
/**************************************************************/ /* UDF for specifying an x-momentum source term */
/**************************************************************/ #include \#define C2 100.0
DEFINE_SOURCE(xmom_source, cell, thread, dS, eqn)
105
{
real x[ND_ND]; real con, source;
C_CENTROID(x, cell, thread);
con = C2*0.5*C_R(cell, thread)*x[1];
source = -con*fabs(C_U(cell, thread))*C_U(cell, thread); dS[eqn] = -2.*con*fabs(C_U(cell, thread)); return source;
}
7.5.2.2 涡产生器
本函数要求处理器返回我们期望的网格变量值。考虑源项的线性化方程，令
S(?)??anb?nb?S(?old)?于是，可得差分方程：
? ? S ? ?a???? ? p???p如果线性化源项中
?S(?p??p,old) ???S??a??S(?)????nbnb?old??p,old???SS(?)?1020(?p,desired??p)，??1020
??则可得该点的参数值为， ? p ,desired相当于该点变量取我们希望的值 。 p,desired?在本例中，我们为了得到r?值，取源项：
source?1020?r??vswirl?
源项的导数：ds＝－1020，由此编写的UDF源程序如下。
/**************************************************************/ /* UDF for specifying a swirl-velocity source term */
/**************************************************************/ #include \
#define OMEGA 50. /* rotational speed of swirler */
#define WEIGHT 1.e20 /* weighting coefficients in linearized equation */ DEFINE_SOURCE(user_swirl, cell, thread, dS, eqn) {
real w_vel, x[ND_ND], y, source; C_CENTROID(x, cell, thread); y = x[1];
w_vel = y*OMEGA; /* linear w-velocity at the cell */ source = WEIGHT*(w_vel - C_WSWIRL(cell,thread)); dS[eqn] = -WEIGHT; return source; }
宏C_WSWIRL(cell，thread)用于取出2D旋转问题的切向速度分量。 7.5.3 用户自定义标量输运方程
本节例子给用户自定义标量输运方程添加源项。本节的用户自定义标量输运方程求解的是FLUENT中的P-1辐射模型，然后将求得的辐射能项加入总的能量方程求解温度。UDFs中定义了边界条件，源项和用户自定义标量输运方程的一些参数。
Poisson入射辐射能量方程，G：
106
入射辐射边界条件：
能量方程源项：
在本例中，我们使用了用户自定义标量输运方程。首先需要在FLUENT中激活用户自定义标量输运方程: Define?Models?User-Defined Scalar…
关于用户自定义函数的具体说明可以参考相应章节。
/**************************************************************/ /* Implementation of the P1 model using user-defined scalars */
/**************************************************************/ #include \#include \#include \
/* Define which user-defined scalars to use. */ enum { P1,
N_REQUIRED_UDS， };
static real abs_coeff = 1.0; /* absorption coefficient */ static real scat_coeff = 0.0; /* scattering coefficient */
static real las_coeff = 0.0; /* linear-anisotropic scattering coefficient */ static real epsilon_w = 1.0; /* wall emissivity */ DEFINE_ADJUST(p1_adjust, domain) {
107
/* Make sure there are enough user defined-scalars. */ if (n_uds < N_REQUIRED_UDS)
Internal_Error(\
}
DEFINE_SOURCE(energy_source, c, t, dS, eqn) {
dS[eqn] = -16.*abs_coeff*SIGMA_SBC*pow(C_T(c,t),3.);
return -abs_coeff*(4.*SIGMA_SBC*pow(C_T(c,t),4.) - C_UDSI(c,t,P1)); }
DEFINE_SOURCE(p1_source, c, t, dS, eqn) {
dS[eqn] = -abs_coeff;
return abs_coeff*(4.*SIGMA_SBC*pow(C_T(c,t),4.) - C_UDSI(c,t,P1)); }
DEFINE_DIFFUSIVITY(p1_diffusivity, c, t, i) {
return 1./(3.*abs_coeff + (3. - las_coeff)*scat_coeff); }
DEFINE_PROFILE(p1_bc, thread, position) {
face_t f;
real A[ND_ND],At;
real dG[ND_ND],dr0[ND_ND],es[ND_ND],ds,A_by_es; real aterm,alpha0,beta0,gamma0,Gsource,Ibw; real Ew = epsilon_w/(2.*(2. - epsilon_w)); Thread *t0=thread->t0;
/* Do nothing if areas aren't computed yet or not next to fluid. */ if (!Data_Valid_P() || !FLUID_THREAD_P(t0)) return; begin_f_loop (f,thread) {
cell_t c0 = F_C0(f,thread);
BOUNDARY_FACE_GEOMETRY(f,thread,A,ds,es,A_by_es,dr0); At = NV_MAG(A);
if (NULLP(T_STORAGE_R_NV(t0,SV_UDSI_G(P1)))) Gsource = 0.; /* if gradient not stored yet */ else
BOUNDARY_SECONDARY_GRADIENT_SOURCE(Gsource,SV_UDSI_G(P1), dG,es,A_by_es,1.);
gamma0 = C_UDSI_DIFF(c0,t0,P1); alpha0 = A_by_es/ds; beta0 = Gsource/alpha0; aterm = alpha0*gamma0/At;
Ibw = SIGMA_SBC*pow(WALL_TEMP_OUTER(f,thread),4.)/M_PI; /* Specify the radiative heat flux. */
F_PROFILE(f,thread,position) =aterm*Ew/(Ew + aterm)*(4.*M_PI*Ibw - C_UDSI(c0,t0,P1) + beta0); }
end_f_loop (f,thread) }
108
函数DEFINE_ADJUST用于检验用户自定义标量输运方程是否已经激活。由于只有一个用户自定义输运方程，引用时为C_UDSI(cell,thread,0)，所以在前面枚举类型定义中使P1为零，这样就可以使用C_UDSI(cell,thread,P1)引用用户自定义的第一个标量输运方程求解的标量。
加入用户自定义输运方程之后，UDFs的函数p1_diffusivity可以在Materials面板选择，p1_source和energy_source可以在Fluid面板选择，p1_bc可以在相应的边界条件面板中选择。
－8 24
P1_source 中宏SIGMA_SB给出Stefan-Boltzmann常数，σ＝5.672×10W/mK；p1_bc中宏FLUID_THREAD_P(t0)检测当前的网格线是否属于流体区域；
BOUNDARY_FACE_GEOMETRY(f,thread,A,ds,es,A_by_es,dr0)计算网格和各面(face area,cell,face)的质心的坐标，以及网格和各面质心之间的距离；宏NULLP(T_STORAGE_R_NV(t0,SV_UDSI_G(P1)))用于检查用户自定义标量输运方程计算的标量梯度是否已经存在。 7.5.4 离散相模型
本节包含FLUENT离散相模型的三个UDFs，由于在每一步计算中都要使用这些UDFs，所以采用Compiled型以节约时间。本节的三个实例为：
1． 沿颗粒轨道的熔化指数（Melting Index） 2． 带电颗粒的磁场力 3． 颗粒拉力系数
7.5.4.1 沿颗粒轨道的熔化指数（Melting Index）
本例使用宏DPM_SCALAR_UPDATE定义的函数，计算沿颗粒轨道的熔化指数（Melting Index）：
meltingindex??t10?dt
宏DEFINE_INIT函数初始化颗粒变量，宏DPM_OUTPUT输出特定面上的熔化指数（Melting Index），宏NULL_P，((p) = =NULL)，检测参数是否为空值。
/***********************************************************************/ /* UDF for computing the melting index along a particle trajectory */
/***********************************************************************/ #include \#include \#include \static real viscosity_0;
DEFINE_INIT(melt_setup, domain) {
/* if memory for the particle variable titles has not been allocated yet, do it now */ if (NULLP(user_particle_vars)) Init_User_Particle_Vars(); /* now set the name and label */
strcpy(user_particle_vars[0].name,\strcpy(user_particle_vars[0].label,\}
/* update the user scalar variables */
DEFINE_DPM_SCALAR_UPDATE(melting_index, cell, thread, initialize, p) {
cphase_state_t *c = &(p->cphase); if (initialize) {
/* this is the initialization call, set:
* p->user[0] contains the melting index, initialize to 0
* viscosity_0 contains the viscosity at the start of a time step */ p->user[0] = 0.;
109
viscosity_0 = c->mu;
} else {
/* use a trapezoidal rule to integrate the melting index */ p->user[0] += P_DT(p) * .5 * (1/viscosity_0 + 1/c->mu); /* save current fluid viscosity for start of next step */ viscosity_0 = c->mu; }
}
/* write melting index when sorting particles at surfaces */
DEFINE_DPM_OUTPUT(melting_output, header, fp, p, thread, plane) {
char name[100]; if (header) {
if (NNULLP(thread))
cxprintf(fp,\else
cxprintf(fp,\cxprintf(fp,\\
\} else {
sprintf(name,\cxprintf(fp,
\\p->state.pos[0], p->state.pos[1], p->state.pos[2], p->state.V[0], p->state.V[1], p->state.V[2],
p->state.diam, p->state.temp, p->flow_rate, p->state.time, p->user[0], name); } }
7.5.4.2 带电颗粒的磁场力
本例计算带电颗粒的磁场力。带电颗粒从上游进入层流流动的流体内，然后在t＝tstart时进入下游的磁场区域，做近似圆周运动（并不完全是圆周运动，因为颗粒所受的磁场力大小是变化的）。
宏P_TIME(p)给出颗粒在流场中运动的时间，p为指针，指向颗粒的轨道。C源程序给出如下。 /***********************************************************************/ /* UDF for computing the magnetic force on a charged particle */
/***********************************************************************/ #include \#include \
#define Q 1.0 /* particle electric charge */
#define BZ 3.0 /* z component of magnetic field */ #define TSTART 18.0 /* field applied at t = tstart */
/* Calculate magnetic force on charged particle. Magnetic */
110
/* force is particle charge times cross product of particle */
/* velocity with magnetic field: Fx= q*bz*Vy, Fy= -q*bz*Vx */ DEFINE_DPM_BODY_FORCE(particle_body_force, p, i) {
float bforce;
if(P_TIME(p)>=TSTART) {
if(i==0) bforce=Q*BZ*P_VEL(p)[1];
else if(i==1) bforce=-Q*BZ*P_VEL(p)[0]; } else
bforce=0.0;
/* an acceleration should be returned */ return (bforce/P_MASS(p)); }
7.5.4.3 颗粒拉力系数
本例计算颗粒的拉力。流场条件5.4.2相同，颗粒所受的拉力与雷诺数有关。 /***********************************************************************/ /* UDF for computing particle drag coefficient (18 Cd Re/24) */ /* curve as suggested by R. Clift, J. R. Grace and M.E. Weber */ /* \
/***********************************************************************/ #include \#include \
DEFINE_DPM_DRAG(particle_drag_force, Re, p) {
float w, drag_force; if (Re < 0.01) {
drag_force=18.0; return (drag_force); }
else if (Re < 20.0) {
w = log10(Re);
drag_force = 18.0 + 2.367*pow(Re,0.82-0.05*w) ; return (drag_force); } else
/* Note: suggested valid range 20 < Re < 260 */ {
drag_force = 18.0 + 3.483*pow(Re,0.6305) ; return (drag_force); } }
7.5.5 计算初始化
初始化函数DEFINE_INIT在FLUENT默认的初始化之后调用。有些流场变量无法在FLUENT默认的初始化中实现，则可以通过初始化函数实现。
/***********************************************************************/ /* UDF for initializing flow field variables */
111
/***********************************************************************/ #include \
DEFINE_INIT(my_init_function, domain) {
cell_t c;
Thread *thread; real xc[ND_ND];
thread_loop_c (thread,domain) {
begin_c_loop_all (c,thread) {
C_CENTROID(xc,c,thread);
if (sqrt(ND_SUM(pow(xc[0] - 0.5,2.), pow(xc[1] - 0.5,2.),
pow(xc[2] - 0.5,2.))) < 0.25) C_T(c,thread) = 400.; else
C_T(c,thread) = 300.; }
end_c_loop_all (c,thread) } }
函数ND_SUM（a，b，c）用于求解函数参数列表中前面两参数（2D），或三参数（3D）之和。两嵌套的loop循环，能够扫描全场，给全场赋初始值。 7.5.6 壁面热流量
宏DEFINE_HEAT_FLUX能够定义壁面与邻近网格之间的热流量。热流量是按照下面两式计算，本例只求解qid 。 qid＝cid[0]＋cid[1]×C_T(c0, t0)－cid[2]×F_T(f, t)－cid[3]×pow(F_T(f,t), 4) qir＝cir[0]＋cir[1]×C_T(c0, t0)－cir[2]×F_T(f, t)－cir[3]×pow(F_T(f,t), 4) C语言源程序为：
/***********************************************************************/ /* UDF for specifying the diffusive heat flux between a wall and neighboring cells */
/***********************************************************************/ #include \
static real h = 0.; /* heat transfer coefficient (W/m^2 C) */ DEFINE_ADJUST(htc_adjust, domain) {
/* Define the heat transfer coefficient. */ h = 120; }
DEFINE_HEAT_FLUX(heat_flux, f, t, c0, t0, cid, cir) {
cid[0] = 0.; cid[1] = h; cid[2] = h; cid[3] = 0.; }
7.5.7 使用用户自定义标量进行后处理
宏DEFINE_ADJUST在每一步迭代执行前调用，所以可以在其中求解标量的微积分，更新与计算结果有关的流场
112
变量等。下例中计算求解温度四次方的导数，并将值存于用户自定义标量之中，以用于后处理。
/***********************************************************************/ /* UDF for computing the magnitude of the gradient of T^4 */
/***********************************************************************/ #include \
/* Define which user-defined scalars to use. */ enum { T4,
MAG_GRAD_T4, N_REQUIRED_UDS }；
DEFINE_ADJUST(adjust_fcn, domain) {
Thread *t; cell_t c; face_t f;
/* Make sure there are enough user-defined scalars. */ if (n_uds < N_REQUIRED_UDS)
Internal_Error(\/* Fill first UDS with temperature raised to fourth power. */ thread_loop_c (t,domain) {
if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4))) {
begin_c_loop (c,t) {
real T = C_T(c,t);
C_UDSI(c,t,T4) = pow(T,4.); }
end_c_loop (c,t) } }
thread_loop_f (t,domain) {
if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4))) {
begin_f_loop (f,t) {
real T = 0.;
if (NULL != THREAD_STORAGE(t,SV_T)) T = F_T(f,t);
else if (NULL != THREAD_STORAGE(t->t0,SV_T)) T = C_T(F_C0(f,t),t->t0); F_UDSI(f,t,T4) = pow(T,4.); }
end_f_loop (f,t) }
113
}
/* Fill second UDS with magnitude of gradient. */ thread_loop_c (t,domain) {
if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)) && NULL != T_STORAGE_R_NV(t,SV_UDSI_G(T4))) {
begin_c_loop (c,t) {
C_UDSI(c,t,MAG_GRAD_T4) = NV_MAG(C_UDSI_G(c,t,T4)); }
end_c_loop (c,t) } }
thread_loop_f (t,domain) {
if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)) && NULL != T_STORAGE_R_NV(t->t0,SV_UDSI_G(T4))) {
begin_f_loop (f,t) {
F_UDSI(f,t,MAG_GRAD_T4) = C_UDSI(F_C0(f,t),t->t0,MAG_GRAD_T4); }
end_f_loop (f,t) } } }
条件语句if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)))检测用户自定义标量中是否已存贮温度的四次方；条件语句NULL != T_STORAGE_R_ NV (t,SV_UDSI_G(T4))检测温度四次方的梯度是否已存入用户自定义标量之中。本例使用了用户自定义标量，所以首先要激活两个用户自定义标量输运方程。 7.5.8 可执行UDFs
本例名为demo_calc的UDF，计算打印当前数据库中的最低温度，最高温度和平均温度，并计算下面的函数，存于用户定义的内存变量中。
f?T??T?Tmin
Tmax?Tmin/**********************************************************************/ /* demand.c */
/* UDF to calculate temperature field function and store in */
/* user-defined memory. Also print min, max, avg temperatures. */
/**********************************************************************/ #include \
extern Domain* domain;
DEFINE_ON_DEMAND(demo_calc) {
float tavg = 0.; float tmax = 0.; float tmin = 0.;
float temp,volume,vol_tot; Thread *t;
114
cell_t c;
/* Loop over all cells in the domain */ thread_loop_c(t,domain) {
/* Compute max, min, volume-averaged temperature */
begin_c_loop(c,t) {
volume = C_VOLUME(c,t); /* get cell volume */ temp = C_T(c,t); /* get cell temperature */ if (temp < tmin || tmin == 0.) tmin = temp; if (temp > tmax || tmax == 0.) tmax = temp; vol_tot += volume; tavg += temp*volume;} end_c_loop(c,t) tavg /= vol_tot;
printf(\/* Compute temperature function and store in user-defined memory */ /* (location index 0) */
begin_c_loop(c,t) {temp = C_T(c,t);
C_UDMI(c,t,0) = (temp-tmin)/(tmax-tmin);} end_c_loop(c,t)}}
函数中使用了嵌套的loop循环，循环内部的函数体实现温度最值，平均值和函数f(T)值的求解。由于函数中用到了用户自定义内存变量，所以需要事先激活用户自定义内存变量，内存变量个数不应小于所使用的内存变量个数。
第六节 UDFs的应用
本章列举了UDFs的几个应用实例。在例子中，我们将边界条件UDFs，源项UDFs，物理性质UDFs和反应速率UDFs应用于实际的相对较简单的物理模型。我们分析这些实例的目的，是为了读者更好的了解UDFs的功能，并能熟练的掌握UDFs。
7.6.1 边界条件UDFs的应用
本节包含两个Interpreted型的例子： 1． 涡轮叶片入口速度分布 2． 管道入口瞬态速度分布 7.6.1.1 涡轮叶片入口速度分布
如图6.1.1，涡轮周围流体网格采用三角无结构网格，定义涡轮腔体左端为速度入口，右端为压力出口边界条件。我们分别计算恒定和抛物型入口速度边界分布的情况，比较两种情况下流场速度的分布。
图6.1.2和图6.1.3为恒定速度入口边界条件计算的结果。入口速度大小为20 m/sec，可以看出恒定速度场在涡轮叶片周围产生扭曲变形。
115
图6.1.1 涡轮叶片周围的网格分布
图6.1.2 恒定速度入口条件下的速度大小分布
图6.1.3 恒定速度入口条件下的速度矢量分布
抛物型速度入口边界条件速度分布满足下面的关系式：
?y? vx?20?20??
?0.0745?其中，y在入口中央取值为0.0，在入口上部和底部分别为＋0.0745，－0.0745，根据入口速度分布式速度中央大小为20 m/sec，上下部值都为20 m/sec。
定义抛物型速度入口边界条件的C源程序如下：
/*************************************************************************/ /* vprofile.c */
/* UDF for specifying a steady-state velocity profile boundary condition */
/*************************************************************************/ #include \
DEFINE_PROFILE(inlet_x_velocity, thread, position) {real x[ND_ND]; /* this will hold the position vector */ real y;
116
2face_t f;
begin_f_loop(f, thread)
{F_CENTROID(x,f,thread); y = x[1];
F_PROFILE(f, thread, position) = 20. - y*y/(.0745*.0745)*20.;} end_f_loop(f, thread) }
本例使用Interpreted型UDFs，按照前面相关章节叙述对C源程序进行编译连接，操作在Interpreted UDFs面板中完成。然后在Velocity Inlet面板的X_Velocity下拉列表中，选择udf inlet_x_velocity，在计算的时候FLUENT就能按照抛物性速度入口边界条件进行计算，计算结果见图6.1.4和图6.1.5。
图6.1.4 抛物型速度入口条件下的速度大小分布
117
图6.1.5 抛物型速度入口条件下的速度矢量分布
7.6.2.1 管道入口瞬态速度分布
在本例中，入口速度满足关系式：vx?20?5sin(10t)，管道长1m，半径0.2m，管道内流体为空气，密度1 kg/m3，粘性系数2×105 kg/m-s。入口速度分布与时间有关，随时间正弦变化。C源程序如下： /***********************************************************************/ /* unsteady.c */
/* UDF for specifying a transient velocity profile boundary condition */
/***********************************************************************/ #include \
DEFINE_PROFILE(unsteady_velocity, thread, position) {face_t f;
begin_f_loop(f, thread)
{real t = RP_Get_Real(\
F_PROFILE(f, thread, position) = 20. + 5.0*sin(10.*t);} end_f_loop(f, thread)}
函数名为unsteady_velocity，变量flow-time存贮流场当前历经时间，函数RP_GET_REAL得到当前时间，由于本例为非稳态问题，需要选择非稳态解法器，如下激活非稳态解法器：
Define?Models?Solver…
本例使用Interpreted型UDFs，在面板Interpreted UDFs中编译连接UDF，在面板Velocity Inlet中X－Velocity下拉列表选择函数udf unsteady_velocity。进行FLUENT默认初始化之后，就可以在Iterate面板进行迭代计算。
Solve?Iterate…
－
118
本例取时间步长为0.0314s，总计算次数为60，则计算总时间为0.0314×60s。每时间步长最多迭代20次，每次迭代前都要更新边界条件并且输出结果。
计算60步之后，我们就可以检查压力出口的速度大小。如果要在计算过程中，查看相关信息，需要事先设定。打开Surface Monitors面板：
Solve?Monitors?Surface…
把Surface Monitors设为1，monitor－1就可以使用了。我们可以输入新文件名代替monitor－1，然后选择是Plot，Print，还是Write。Every有两个选项，分别为Iteration，Flow Time或Flow Time，用来定义监测流场相关变量的时间间隔。点击Define…之后，出现Define Surface Monitor 面板：
在相应下拉列表中选择Velocity…和Velocity Magnitude，选择所要输出参数的面pressure –outlet－5。监测输出参数类型为平均值（Average），x轴取流动时间。
设置好之后，就会在每一步时间步长都输出速度大小，以供计算时检测。我们也可以通过File XY Plot面板，手工绘制文件monitor-1.out的图形。
Plot?File…
119
在Files里选择文件，如果需要添加文件点击按钮add…，选择文件之后
点击Plot按钮，可以得到文件输出的图形。图6.1.6输出不同时刻本例压力出口的速度大小。
图6.1.6 不同时刻压力出口的速度大小
可以看出，正如我们所料，压力出口速度也是以20 m/sec为平衡位置，振幅为5 m/sec的周期性分布。 7.6.2 源项UDFs的应用
本节实例讲述源项UDFs的应用，执行方式为Interpreted型。 7.6.2.1 添加管内流动源项
本例添加2D管道流动动量方程的源项。管道长4m，宽2m，液态金属（物性见表6.2.1）从左端以1mm/s的速度进入，温度为290K。进入管道0.5m之后，管道上壁开始对液态金属冷却，上壁保持恒温280K。为了模拟液态金属的凝固，当液态金属温度低于288K时，就添加动量方程的源项。在本例中，我们仅以液态金属的速度大小来表征是否凝固，如果速度等于零，则认为液态金属已经凝固，而不是求解能量方程。温度低于288K的液态金属，由于添加了源项，速度会逐渐降为零。
我们添加的是x向动量方程的源项，源项方程为：
Sx??Cvx ，
?Sx??C ?vx 表6.2.1 液态金属物性参数
120
Property Density Viscosity Specific Heat Thermal Conductivity Value 8000 kg/m3 5.5×10-3 kg/m-s 680 J/kg-K 30 W/m-K 添加源项UDF的C源程序如下： #include \#define CON 20.0
DEFINE_SOURCE(cell_x_source, cell, thread, dS, eqn) {real source;
if (C_T(cell,thread) <= 288.) {/* source term */
source = -CON*C_U(cell,thread);
/* derivative of source term w.r.t. x-velocity. */
dS[eqn] = -CON;} else
source = dS[eqn] = 0.; return source;}
本例使用Interpreted Udfs，编译连接之后，在Fluid面板中的X Momentum下拉列表可以选择函数cell_x_xource。
计算收敛之后，我们可以绘出管道内的温度、速度、流量函数，如图6.2.1，6.2.2和6.2.3。
121
图6.2.1 管道内温度分布
图6.2.2 管道内速度大小分布
图6.2.3 管道内流函数分布
122
如图可以明显看出管道内凝固的区域。本例只是凝固情况的一种简单近似，更精确的计算需要添加能量方程的源项等。
7.6.3 物性UDFs的应用
本节实例讲述物性UDFs的应用，执行方式为Interpreted型。 7.6.3.1 粘性系数与温度相关流体的凝固
－
对于暖流体（T > 288K）粘性系数为5.5×103kg/m-s，冷流体（T < 286K）粘性系数为1.0kg/m-s，处于中间温度（286K ≤T≤ 288K）的流体，粘性系数满足关系式： 。本例以粘性系数的值反应流体的凝
??143.2135?0.49725T固，因为流体从液态变为固态固粘性系数会陡增。
UDF的C语言源程序为：
/***********************************************************************/ /* viscosity.c */
/* UDF for specifying a temperature-dependent viscosity property */
/***********************************************************************/ #include \
DEFINE_PROPERTY(cell_viscosity, cell, thread) {real mu_lam;
real temp = C_T(cell, thread); if (temp > 288.) mu_lam = 5.5e-3; else if (temp > 286.)
mu_lam = 143.2135 - 0.49725 * temp; else
mu_lam = 1.; return mu_lam;}
变量temp存放温度，mu_lam返回粘性系数值。物性定义UDFs编译连接之后，可以在Material面板中选择设置。
点击Edit…
123
选择cell_viscosity。
结果见图6.3.1，6.3.2和6.3.3，可以看出粘性系数在极小的区域内从0.0055增至1.0 kg/m-s。从图6.3.2还可以看出，速度并不像6.2.1中的结果一样陡降为零，而是逐渐降低，存在一个很大的“模糊区”。
图6.3.1 粘性系数分布
图6.3.2 速度大小分布
124
图6.3.3 流函数分布
7.6.4 反应速率UDFs的应用
本节实例讲述反应速率UDFs的应用，执行方式为Compiled型。 7.6.4.1 体积反应速率
对于简单的双气相组分系统，组分a生成组分b的反应速率公式为：
R?K1Xa?1?K2Xa?2
其中，Xa为组分a的质量分数，K1和K2为常数。
2D管道长114英寸，宽16英寸，在中部有90度的拐角（如图6.4.1）。管道底部和右端覆盖6英寸厚的多孔介质，
－
反应发生在多空介质区。管道内组分有相同物理性质，密度为1.0 kg/m，粘性系数为1.72×105 kg/m-s。
图6.4.1 2D管道略图
气体a以0.1 m/s的速度从左端进入管道，经过多空介质区和管道空间，多孔介质区各向的惯性阻力都为5 m-1。从图6.4.2的流场可以看出大多数流体是从管道空间流出的，经过多孔介质区的只是很小一部分。
125
图6.4.2 具有多空介质区的2D管道内的流线分布
从图6.4.3速度矢量可以看出，多孔介质内的流动比管道空间的流动要慢得多。
图6.4.3 具有多空介质区的2D管道内的速度矢量分布 本例的C源程序为：
/**************************************************************/ /* rate.c */
/* UDF for specifying a reaction rate in a porous media */
/**************************************************************/ #include \#define K1 2.0e-2 #define K2 5.
DEFINE_VR_RATE(user_rate, cell, thread, r, mole_weight, species_mf, rate, rr_t) {real s1 = species_mf[0];
real mw1 = mole_weight[0];
if (FLUID_THREAD_P(thread) && THREAD_VAR(thread).fluid.porous) *rate = K1*s1/pow((1.+K2*s1),2.0)/mw1; else
*rate = 0.;}
126
函数FLUID_THREAD_P(thread)判断当前网格线是否在多空介质区，变量 THREAD_VAR(thread).fluid.porous用于检测当前网格是否属于多空介质区。
按照第三章，编译建立共享库之后，就可以在Compiled UDFs面板打开共享库连接UDFs。连接成功之后，可以在User－Defined Function Hooks面板选择使用反应速率UDFs。
Define?User-Defined?Function Hooks…
初始化和计算之后，我们可以得到收敛的结果，图6.4.4给出了组分a的质量分数。从图中可以看出，多空介质区组分a占0％，说明组分a都已经生成了组分b；管道空间组分a占100％，说明组分a没有进行反应。
图6.4.4 具有多空介质区的2D管道内的组分a质量分数分布
附录A udf.h中的宏解释
A.1 General-Use DEFINE Macros
#define DEFINE_DIFFUSIVITY(name, c, t, i)
real name(cell_t c, Thread *t, int i)
#define DEFINE_PROFILE(name, t, i) \\
void name(Thread *t, int i)
#define DEFINE_PROPERTY(name, c, t) \\
real name(cell_t c, Thread *t)
#define DEFINE_SOURCE(name, c, t, dS, i) \\
real name(cell_t c, Thread *t, real dS[], int i)
#define DEFINE_INIT(name, domain) \\
void name(Domain *domain)
#define DEFINE_ADJUST(name, domain) \\
void name(Domain *domain)
#define DEFINE_UDS_FLUX(name, f, t, i) \\
127
real name(face_t f, Thread *t, int i)
#define DEFINE_UDS_UNSTEADY(name, c, t, i, apu, su) \\
void name(cell_t c, Thread *t, int i, real *apu, real *su)
#define DEFINE_HEAT_FLUX(name, f, t, c0, t0, cid, cir) \\
void name(face_t f, Thread *t, cell_t c0, \\ Thread *t0, real cid[], real cir[])
#define DEFINE_VR_RATE(name, c, t, r, mw, yi, rr, rr_t) \\
void name(cell_t c, Thread *t, \\ Reaction *r, real *mw, real *yi, \\ real *rr, real *rr_t)
#define DEFINE_SR_RATE(name, f, t, r, mw, yi, rr) \\
void name(face_t c, Thread *t, \\
Reaction *r, real *mw, real *yi, real *rr)
#define DEFINE_SCAT_PHASE_FUNC(name, c, f) \\
real name(real c, real *f)
#define DEFINE_RW_FILE(name, fp) \\
void name(FILE *fp)
#define DEFINE_ON_DEMAND(name) \\
void name(void)
A.2 Discrete Phase DEFINE Macros
#define DEFINE_DPM_BODY_FORCE(name, p, i) \\
real name(Tracked_Particle *p, int i)
#define DEFINE_DPM_DRAG(name, Re) \\
real name(real Re)
#define DEFINE_DPM_SOURCE(name, c, t, S, strength, p) \\
void name(cell_t c, Thread *t, dpms_t *S, \\ real strength, Tracked_Particle *p)
#define DEFINE_DPM_PROPERTY(name, c, t, p) \\
real name(cell_t c, Thread *t, Tracked_Particle *p)
#define DEFINE_DPM_OUTPUT(name, header, fp, p, t, plane) \\
void name(int header, FILE *fp, \\
Tracked_Particle *p, Thread *t, Plane *plane)
#define DEFINE_DPM_EROSION(name, p, t, f, normal, alpha, Vmag, mdot) \\
void name(Tracked_Particle *p, Thread *t, \\ face_t f, real normal[], real alpha, \\ real Vmag, real mdot)
#define DEFINE_DPM_SCALAR_UPDATE(name, c, t, initialize, p) \\
void name(cell_t c, Thread *t, int initialize, \\ Tracked_Particle *p)
#define DEFINE_DPM_LAW(name, p, ci)
void name(Tracked_Particle *p, int ci)
#define DEFINE_DPM_SWITCH(name, p, ci) \\
void name(Tracked_Particle *p, int ci)
#define DEFINE_DPM_INJECTION_INIT(name, I) \\
void name(Injection *I)
A.3 Multiphase DEFINE Macros
#define DEFINE_DRIFT_DIAM(name, c, t) \\
real name(cell_t c, Thread *t)
#define DEFINE_SLIP_VELOCITY(name, domain) \\
void name(Domain *domain)
128