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

/* 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; 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，分布式为：

x?? T(x)?300?10s0in???

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 假设 sourc?eS??Avxvx ，A?0.5C2?y 于是，有

dSd?vx? ??Avx?Avxdvxdvx所以，UDF的返回源项值为： sourc?e?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) {

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 涡产生器

本函数要求处理器返回我们期望的网格变量值。考虑源项的线性化方程，令

?SS(?)??anb?nb?S(?old)?(?p??p,old)

??于是，可得差分方程：

? ? S ? ?S???p,old??ap???p??anb?nb??S(?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：

入射辐射边界条件：

能量方程源项：

在本例中，我们使用了用户自定义标量输运方程。首先需要在FLUENT中激活用户自定义标量输运方程: