#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)    }