UDF的宏用法及相关算例

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

(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源文件不在当前工作目录下，编译时需要输入文件路径。如：

\\\\\\Fluent.Inc\%udf\\myudf\\example.c fileserver表示安装FLUENT的计算机名。

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

函数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

本文来源：https://www.bwwdw.com/article/xubg.html