中科大FLUENT讲稿 第七章 自定义函数

宏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 */ /* 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 */

/***********************************************************************/ #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在每一步迭代执行前调用，所以可以在其中求解标量的微积分，更新与计算结果有关的流场变量等。下例中计算求解温度四次方的导数，并将值存于用户自定义标量之中，以用于后处理。

/***********************************************************************/ /* 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) } }

/* Fill second UDS with magnitude of gradient. */ thread_loop_c (t,domain) {

if (NULL != THREAD_STORAGE(t,SV_UDS_I(T4)) &&