UDF案例

nullFLUENT UDF案例FLUENT UDF案例飞昂软件技术（上海）有限公司 ANSYS　FLUENT CHINA 2008-12内容概要内容概要UDS——低Re湍流模型添加 动网格 Degassing Boundary（多相流） 用户自定义的Drag Law （多相流） UDS——低Re湍流模型添加UDS——低Re湍流模型添加问题描述问题描述通过UDF添加低雷诺数模型，AKN model (Abe, Kondoh, Nagano(IJHMT, 1995)) 需要添加k和epsion两个方程壁面边界条件：需要用到的数据宏需要用到的数据宏密度 C_R(c,t) 层流黏度 C_MU_L(c,t) 湍动能 C_UDSI(c,t,TKE) 湍流耗散率 C_UDSI(c,t, TDR) 平均速度应变率 C_STRAINRATE_MAG(c,t) 湍流黏度 C_MU_T(c,t) 程序实现程序实现#include “udf.h” const real C_MU=0.09; const real SIG_TDR=1.4; const real SIG_TKE=1.4; const real C1_D=1.5; const real C2_D=1.9; Const real f_1=1.; /* User-defined Scalars */ enum { TKE, TDR, N_REQUIRED_UDS }; real y_star (cell_t c, Thread *t) { return C_R(c,t)*pow(C_UDSI(c,t,TDR)*C_MU_L(c,t)/C_R(c,t), 0.25)*C_WALL_DIST(c,t)/C_MU_L(c, t); } real Re_t (cell_t c, Thread *t) { return C_R(c,t)*C_UDSI(c,t,TKE)*C_UDSI(c,t,TKE) /C_MU_L(c,t)/C_UDSI(c,t,TDR); } 程序实现程序实现real f_mu (cell_t c, Thread *t) { return SQR(1. - exp(-y_star(c,t)/14.))*(1.+ 5./pow(Re_t(c,t), 0.75)* exp( -SQR(Re_t(c,t)/200.) )); } real f_2(cell_t c, Thread *t) { return SQR( 1.0 - exp(-y_star(c,t)/3.1 )) *(1.-0.3*exp(-SQR(Re_t(c,t)/6.5))); } 程序实现程序实现K方程源项 线性化DEFINE_SOURCE(tke_source, c, t, dS, eqn) { real rho_Pk, rho_epsilon; rho_Pk=C_MU_T(c,t)*C_STRAIN_RATE_MAG(c,t) *C_STRAIN_RATE_MAG(c,t); rho_epsilon = SQR(C_R(c,t))*C_MU*f_mu(c,t) * C_UDSI(c,t,TKE)*C_UDSI(c,t,TKE) /C_MU_T(c,t); dS[eqn]=-2.*SQR(C_R(c,t)) *C_MU*f_mu(c,t)*C_UDSI(c,t,TKE)/C_MU_T(c,t); return rho_Pk - rho_epsilon; } 程序实现程序实现Epsilon方程源项及其线性化DEFINE_SOURCE(tdr_source, c, t, dS, eqn) { real Pk, source; Pk = C_MU_T(c,t)/C_R(c,t)* SQR(C_STRAIN_RATE_MAG(c,t)); dS[eqn]=C1_D*f_1*Pk/C_UDSI(c,t,TKE) -2.*C_R(c,t)*C_UDSI(c,t,TDR)* f_2(c,t)*C2_D/C_UDSI(c,t,TKE); source=C_R(c,t)*C_UDSI(c,t,TDR)/ C_UDSI(c,t,TKE)*( f_1*C1_D*Pk- C2_D*f_2(c,t)*C_UDSI(c,t,TDR)); return source; } 程序实现程序实现扩散项 DEFINE_DIFFUSIVITY(ke_diffusivity, c, t, i) real diff; real mu_t; int nscalar=i; mu_t=C_R(c,t)*C_MU*f_mu(c,t)*C_UDSI(c,t,TKE)* C_UDSI(c,t,TKE)/C_UDSI(c,t,TDR); if ( nscalar==TKE ) diff=mu_t/SIG_TKE + C_MU_L(c,t); else diff=mu_t/SIG_TDR + C_MU_L(c,t); return diff; } 程序实现程序实现返回湍流黏度 DEFINE_TURBULENT_VISCOSITY(turbVis,c,t) { real mu_t; mu_t=C_R(c,t)*C_MU*f_mu(c,t)* SQR(C_UDSI(c,t,TKE))/C_UDSI(c,t,TDR); return mu_t; } 程序实现程序实现壁面边界条件 DEFINE_PROFILE(tdr_wall, thread, position) { real xw[ND_ND], xc[ND_ND], delta[ND_ND], dy; face_t f; cell_t c0; Thread *t0=THREAD_T0(thread); begin_f_loop(f,thread) { c0=F_C0(f,thread); C_CENTROID(xc, c0, t0); F_CENTROID(xw, f, thread); NV_VV(delta,=, xc, -, xw); dy=NV_MAG(delta); F_PROFILE(f, thread, position)=2.* ( C_MU_L(c0,t0)/C_R(c0,t0) )*C_UDSI(c0,t0,TKE)/SQR(dy); }end_f_loop(f,thread); } 程序实现程序实现动网格UDF案例动网格UDF案例问题运动问题运动绿色的阀体作定速转动 黄色膜片作随时间做正/余旋的变形 动网格区域 阀体——DEFINE_CG_MOTION 膜片—— DEFINE_GRID_MOTION 上部壁面——DEFINE_GEOMDEFINE_CG_MOTIONDEFINE_CG_MOTION阀体转动 #define omega 1.0 DEFINE_CG_MOTION(butterfly_flex_UDF, dt, cg_vel, cg_omega, time, dtime) { cg_vel[0] = 0.0; cg_vel[1] = 0.0; cg_vel[2] = 0.0; cg_omega[0] = 0.0; cg_omega[1] = 0.0; cg_omega[2] = omega; } DEFINE_GRID_MOTIONDEFINE_GRID_MOTIONDEFINE_GRID_MOTION(moving_arc, domain, dt, time, dtime) { Thread *tf = DT_THREAD (dt); face_t f; Node *node_p; real alpha, theta, x, ymag, yfull, y; int n; SET_DEFORMING_THREAD_FLAG (THREAD_T0 (tf)); alpha = omega * CURRENT_TIME; theta = 2.0 * alpha + 3.0 * M_PI / 2.0; begin_f_loop (f, tf) { f_node_loop (f, tf, n) { node_p = F_NODE (f, tf, n); if (NODE_POS_NEED_UPDATE (node_p)) { NODE_POS_UPDATED (node_p); x = NODE_X (node_p); ymag = sqrt (R*R - x*x) + 0.03; yfull = ymag - 0.1; y = - 0.1 + yfull * sin(theta); NODE_Y (node_p) = y; } } } end_f_loop (f, tf); } DEFINE_GEOMDEFINE_GEOM#define R 0.109 DEFINE_GEOM(plane, domain, dt, position) { position[1] = R; }Degassing BoundaryDegassing Boundary问题描述问题描述多相流问题中，出于计算的目的，同一个边界，有时希望对其中某个相是出口，对其它相则是墙或是进口等边界。 本例中： 下部边界是气体的进口，液滴的出口思路思路对于气相，底部面是标准的进口。 对液滴相，靠近底部的最近一层网格需要添加一个汇，以达到质量守恒。 主相无需作出调整。程序实现程序实现DEFINE_SOURCE(degassing_source, cell, thread, dS, eqn) { real source; Thread *tm = THREAD_SUPER_THREAD(thread); source = -C_R(cell,thread)*C_VOF(cell,thread)/CURRENT_TIMESTEP ; C_UDMI(cell,tm,0) = source; dS[eqn] = -C_R(cell,thread)/CURRENT_TIMESTEP; return source; }守恒判断守恒判断相当结果结果3 m/s7 m/sairdrop反作用力反作用力没有考虑反作用力。 DEFINE_SOURCE(x_prim_recoil, cell, tp, dS, eqn) { real source; Thread *tm = THREAD_SUPER_THREAD(tp); Thread *ts; ts = THREAD_SUB_THREAD(tm,1); source = -C_R(cell,ts)*C_VOF(cell,ts)/CURRENT_TIMESTEP*C_U(cell,tp); dS[eqn] =-C_R(cell,ts)*C_VOF(cell,ts)/CURRENT_TIMESTEP; return source; } DEFINE_SOURCE(x_sec_recoil, cell, ts, dS, eqn) { real source; Thread *tm = THREAD_SUPER_THREAD(ts); source = -C_R(cell,ts)*C_VOF(cell,ts)/CURRENT_TIMESTEP*C_U(cell,ts); dS[eqn] = -C_R(cell,ts)*C_VOF(cell,ts)/CURRENT_TIMESTEP; return source; }改进后的结果改进后的结果2.6 m/s2.7 m/sairdrop用户自定义的Drag Law用户自定义的Drag Law问题描述问题描述欧拉多相流模型中的Drag Law多为半经验模型，需要根据具体问题进行调整。 Syamlal-O‘brien模型的最小流化速度大于21cm/s。 本例中，最小流化速度为8cm/s。 Drag Law为如下形式： 其中程序实现程序实现#include "udf.h" #include "sg_mphase.h" #define diam2 3.e-4 DEFINE_EXCHANGE_PROPERTY(custom_drag_syam, cell, mix_thread, s_col, f_col) { Thread *thread_g, *thread_s; real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y, rho_g, rho_s, mu_g, reyp, afac, bfac, void_g, vfac, fdrgs, taup, k_g_s; thread_g = THREAD_SUB_THREAD(mix_thread, s_col thread_s = THREAD_SUB_THREAD(mix_thread, f_col); x_vel_g = C_U(cell, thread_g); y_vel_g = C_V(cell, thread_g); x_vel_s = C_U(cell, thread_s); y_vel_s = C_V(cell, thread_s); slip_x = x_vel_g - x_vel_s; slip_y = y_vel_g - y_vel_s; rho_g = C_R(cell, thread_g); rho_s = C_R(cell, thread_s); mu_g = C_MU_L(cell, thread_g); /*compute slip*/ abs_v = sqrt(slip_x*slip_x + slip_y*slip_y); /*compute reynolds number*/ reyp = rho_g*abs_v*diam2/mu_g; /* compute particle relaxation time */ taup = rho_s*diam2*diam2/18./mu_g; void_g = C_VOF(cell, thread_g);/* gas vol frac*/ /*compute drag and return drag coeff, k_g_s*/ afac = pow(void_g,4.14); if(void_g<=0.85) bfac = 0.281632*pow(void_g, 1.28); else bfac = pow(void_g, 9.076960); vfac = 0.5*(afac-0.06*reyp+sqrt(0.0036*reyp*reyp+0.12*reyp*(2.*bfac- afac)+afac*afac)); fdrgs = void_g*(pow((0.63*sqrt(reyp)/vfac+4.8*sqrt(vfac)/vfac),2))/24.0; k_g_s = (1.-void_g)*rho_s*fdrgs/taup; return k_g_s; }