udf 调试

40534427

悬赏一个1技术分。
这是一个自定义湍流的 udf，很遗憾作者没有给全。放上来给有兴趣和这方面专长的坛友一起研究一下，填写出正确的/* do something... */

#include "udf.h"
/* Turbulence model constants */
#define  C_MU  0.09
#define  SIG_K 1.0
#define  SIG_D 1.3
#define  C1_D  1.44
#define  C2_D  1.92
/* User-defined scalars */
enum
{
  K,
  D,
  N_REQUIRED_UDS
};
/* Reynolds number definitions */
real Re_y(cell_t c, Thread *t)
{
  real c_k = MAX(C_UDSI(c,t,K), EPSILON);
  return C_R(c,t)*sqrt(c_k)*C_WALL_DIST(c,t)/C_MU_L(c,t);
}
real Re_t(cell_t c, Thread *t)
{
  real c_k, c_d;
  
  c_k = MAX(C_UDSI(c,t,K), EPSILON);
  c_d = MAX(C_UDSI(c,t,D), EPSILON);
  return C_R(c,t)*SQR(c_k)/C_MU_L(c,t)/c_d;
}
/* Damping Functions */
real f_mu(cell_t c, Thread *t)
{
  return tanh(0.008*Re_y(c,t))*(1.+4/pow(MAX(Re_t(c,t),EPSILON),0.75));
}

real f_1(cell_t c, Thread *t)
{
  return 1.;
}
real f_2(cell_t c, Thread *t)
{
  return (1.-2/9*exp(-MIN(Re_t(c,t)*Re_t(c,t)/36,10)))
    *(1.-exp(-MIN(Re_y(c,t)/12,10)));  
}
DEFINE_SOURCE(k_source, c, t, dS, eqn)
{
  real c_k, c_d;
  real mu_t, G_k;
  
  c_k = MAX(C_UDSI(c,t,K), EPSILON);
  c_d = MAX(C_UDSI(c,t,D), EPSILON);
  mu_t = C_MU_T(c,t) + C_MU_L(c,t)*1.0E-10;
  G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t));
  
  dS[eqn] = -2.*C_R(c,t)*C_R(c,t)*C_MU*f_mu(c,t)*c_k/mu_t;
  return G_k - C_R(c,t)*C_R(c,t)*C_MU*f_mu(c,t)*SQR(c_k)/mu_t;
}
DEFINE_SOURCE(d_source, c, t, dS, eqn)
{
  real c_k, c_d;
  real G_k;
  
  c_k = MAX(C_UDSI(c,t,K), EPSILON);
  c_d = MAX(C_UDSI(c,t,D), EPSILON);
  G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t));
  dS[eqn] = C1_D*f_1(c,t)*G_k/c_k
     - 2.*C2_D*f_2(c,t)*C_R(c,t)*c_d/c_k;
  return C1_D*f_1(c,t)*G_k*c_d/c_k
     - C2_D*f_2(c,t)*C_R(c,t)*SQR(c_d)/c_k;
}
DEFINE_DIFFUSIVITY(ke_diffusivity, c, t, eqn)
{
  real diff;
  real mu = C_MU_L(c,t);
  /* real mu_t =  C_MU_T(c,t); */
  real c_k = MAX(C_UDSI(c,t,K), 0.);
  real c_d = MAX(C_UDSI(c,t,D), EPSILON);
  real mu_t = C_R(c,t)*C_MU*f_mu(c,t)*SQR(c_k)/c_d;
  switch (eqn)
    {
    case K:
      diff = mu_t/SIG_K + mu;
      break;
    case D:
      diff = mu_t/SIG_D + mu;
      break;
    default:
      diff = mu_t + mu;
    }
  return diff;
}
DEFINE_TURBULENT_VISCOSITY(ke_mut, c, t)
  /* Set the turbulent viscosity */
{
  real c_k = MAX(C_UDSI(c,t,K), 0.);
  real c_d = MAX(C_UDSI(c,t,D), EPSILON);
  real mu_t = C_R(c,t)*C_MU*f_mu(c,t)*SQR(c_k)/c_d;
  return mu_t;
}
DEFINE_ADJUST(turb_adjust, domain)
{
  /* do something... */
}

DEFINE_INIT(turb_init, domain)
{
  /* do something... */
}
DEFINE_PROFILE(wall_d_bc, t, position)
{
  face_t f;
  cell_t c0;
  Thread *t0 = t->t0; /* t0 is cell thread */
  real xw[ND_ND], xc[ND_ND], dx[ND_ND], rootk, dy, drootkdy;
  
  /* Do nothing if wall distance not computed or not in fluid zone */
  if (!Data_Valid_P() || !FLUID_THREAD_P(t0)) return;
  begin_f_loop(f,t)
    {
      c0 = F_C0(f,t);
      rootk = sqrt(MAX(C_UDSI(c0,t0,K), 0.));
      F_CENTROID(xw,f,t);
      C_CENTROID(xc,c0,t0);
      NV_VV(dx, =, xc, -, xw);
      dy = ND_MAG(dx[0], dx[1], dx[2]);
      drootkdy = rootk/dy;
      F_PROFILE(f,t,position) = 2.*C_MU_L(c0,t0)/C_R(c0,t0)*drootkdy*drootkdy;
    }
  end_f_loop(f,t)   
}

fox000002

从代码功能上看，UDS 相关的定义是完整的

DEFINE_ADJUST 和 DEFINE_INIT 其实并不是必需的

可以用于检查一下 UDS 数目或做些后处理的功能

比如

1. 
   
2. /*  Make  sure  there  are  enough  user  defined-scalars.  */
   
3. if  (n_uds  <  N_REQUIRED_UDS)
   
4.      Internal_Error("not  enough  user-defined  scalars  allocated");
   

复制代码