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