k-w方程加入源项后出现的问题

duanhuishen

#include "udf.h"
real Re_k = 6.0;
real Beta_inf=0.09;
real a_inf=0.52;
real a0_star=0.024;
real a0=1/9.0;
real Re_w = 2.95;
real Re_Beta=8;
real Beta_inf_star=0.09;
real Re_t(cell_t c, Thread *t)
{return C_R(c,t)*C_K(c,t)/(C_MU_L(c,t)*C_O(c,t));}
real alpha(cell_t c, Thread *t)
{return (0.024+Re_t(c,t)/Re_k)/(1+Re_t(c,t)/Re_k);}
real alpha1(cell_t c, Thread *t)
{ return (a0+Re_t(c,t)/Re_w)/(1+Re_t(c,t)/Re_w);}
real Beta_star(cell_t c, Thread *t)
{ return (4/15+pow(Re_t(c,t)/Re_Beta,4))/(1+pow(Re_t(c,t)/Re_Beta,4));}
real F_Mt(cell_t c, Thread *t)
{   real FM,Mt,a;
    a=sqrt(Gama*R*C_T(c,t));
    Mt=sqrt(2*C_K(c,t))/a;
    if(Mt<Mt0) FM=0;
else FM=SQR(Mt)-SQR(Mt0);
return FM;
}
real F_Beta(cell_t c, Thread *t)
{
  real FB_star,X_k;
  X_k=(1/pow(C_O(c,t),3))*(C_K_G(c,t)[0]*C_O_G(c,t)[0]+C_K_G(c,t)[1]*C_O_G(c,t)[1]);
  if(X_k>0) FB_star=(1+680*SQR(X_k))/(1+400*SQR(X_k));
  else FB_star=1;
  return FB_star;
}
// Set k equation Source
DEFINE_SOURCE(k_source,c,t,dS,eqn)
{
real G_k,D_k,D_k1,D_k2,source;
dS[eqn]=-Beta_star(c,t)*C_R(c,t)*Beta_inf_star*F_Beta(c,t)*C_O(c,t);
source=(1-Beta_star(c,t))*C_R(c,t)*Beta_inf_star*F_Beta(c,t)*C_K(c,t)*C_O(c,t);
return source;
}
// Set w equation Source
DEFINE_SOURCE(w_source,c,t,dS,eqn)
{
real G_w,D_w,G_k,source,mu_t0;
G_k=C_MU_T(c,t)*pow(C_STRAIN_RATE_MAG(c,t),2);
G_w=a_inf*G_k*C_O(c,t)/C_K(c,t);
dS[eqn]=-1.0*a_inf*G_k/C_K(c,t);
source=(alpha1(c,t)/alpha(c,t)-1.0)*G_w;
return source;
}
// Set Turbulent viscosity
DEFINE_TURBULENT_VISCOSITY(user_mu_t,c,t)
{
real mu_t;
real rho=C_R(c,t);
real k=C_K(c,t);
real o=C_O(c,t);
real a=alpha(c,t);
mu_t=a*rho*k/o;
return mu_t;
}
上述UDF通过，这个UDF实际上就是不选者transitonal flow的情况下，按照帮助中低雷诺数修正方法，验证一下转捩模拟(为进一步的工作奠定基础)，但是在计算的过程中出现问题如下图，希望同行能给与指点！

duanhuishen

计算发散了(不可压)，请大家帮忙看一下上述UDF右什么问题？