abid.c

fluent湍流模型代码,基于fluent软件。主要是利用标量方程修改里面的湍流模型。

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