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;
}