24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

进击的imp

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1.5
帖子: 7
在线: 15.5小时
虫号: 2834096
注册: 2013-11-28
专业: 内流流体力学

[求助] 新手自己写了个UDF，欲替代FLUENT里的S-A湍流模型，总是出现问题，求高手指点

#include "udf.h"

/*Turbulence model constants*/
#define C_B1 0.1355
#define C_B2 0.622
#define SIG_V  0.66667
#define C_V1 7.1
#define C_W2 0.3
#define C_W3 2.0
#define K    0.4187
#define C_W1 3.206

/*Modeling the turbulent viscosity*/

real X(cell_t c, Thread *t)
{ return C_UDSI(c,t,0)/C_MU_L(c,t)*C_R(c,t);}

real f_v1(cell_t c, Thread *t)
{ return X(c,t)*X(c,t)*X(c,t)/(X(c,t)*X(c,t)*X(c,t)+C_V1*C_V1*C_V1);}

DEFINE_TURBULENT_VISCOSITY(sa_mut,c,t)
{
return C_R(c,t)*f_v1(c,t)*C_UDSI(c,t,0);
}

/*Modeling the turbulent production*/

real f_v2(cell_t c, Thread *t)
{ return 1.0-X(c,t)/(1.0+X(c,t)*f_v1(c,t));}

real ome(cell_t c, Thread *t)
{ return sqrt(SQR(C_DUDY(c,t)-C_DVDX(c,t))+SQR(C_DUDZ(c,t)-C_DWDX(c,t))+SQR(C_DVDZ(c,t)-C_DWDY(c,t)));
}

real S(cell_t c, Thread *t)
{return ome(c,t)+2.0*MIN(0.0,(C_STRAIN_RATE_MAG(c,t)-ome(c,t)));}

real S_W(cell_t c, Thread *t)
{ return S(c,t)+C_UDSI(c,t,0)/K/K/C_WALL_DIST(c,t)/C_WALL_DIST(c,t)*f_v2(c,t);}

real G_V(cell_t c,Thread *t)
{return C_B1*C_R(c,t)*S_W(c,t)*C_UDSI(c,t,0);}

/*Modeling the turbulent destruction*/
real r(cell_t c, Thread *t)
{return C_UDSI(c,t,0)/S_W(c,t)/K(c,t)/K(c,t)/C_WALL_DIST(c,t)/C_WALL_DIST(c,t);}

real g(cell_t c, Thread *t)
{return r(c,t)+C_W2*(pow(r(c,t),6)-r(c,t));}

real f_w(cell_t c, Thread *t)
{return g(c,t)*pow((1.0+pow(C_W3,6)/(pow(r(c,t),6)+pow(C_W3,6))),0.166667);}

real Y_V(cell_t c, Thread *t)
{return C_W1*C_R(c,t)*f_w(c,t)*SQR(C_UDSI(c,t,0)/C_WALL_DIST(c,t));}

/*Source term*/
DEFINE_SOURCE(sa_source,c,t,dS,eqn)
{
dS[eqn]=0;
return  G_V(c,t)-Y_V(c,t)+C_B2*C_R(c,t)/SIG_V*(SQR(C_UDSI_G(c,t,0)[0])+SQR(C_UDSI_G(c,t,0)[1])+SQR(C_UDSI_G(c,t,0)[2]));
}

/*Diffusivity*/

DEFINE_DIFFUSIVITY(sa_diffusivity,c,t,eqn)
{
   return (C_MU_L(c,t)+C_R(c,t)*C_UDSI(c,t，0))/SIG_V;
}
基本是按照斯坦福的一个ppt里写KE模型的udf改写的，参考FLUENT帮助文件里SA模型的方程。但是编译时有很多问题。另外，dS[eqn]没法确切表达，师兄说可以设成0，不知道是否合理。不胜感激。

回复此楼