Standard k-eps with UDF/UDS
 

Hi,
I am trying to implement turbulence models in FLUENT using UDS/UDF. As a first step I tried with just implementing the Standard k-epsilon. However I cannot get it to work. In the first step I add two UDS variables and let them have the same diffusivites and sources as k and epsilon but I use the ordinary k and epsilon in the equations for UDS0 and UDS1. Then I patch UDS-0 and UDS-1 with k and epsilon respectively. What I would expect is that my UDS-0 and UDS-1 would keep like k and epsilon. Instead they quikly get very different and I even get negative values for UDS-1 (“epsilon”).
I would be very thankful if someone could have a look and see if I am missing something in my implementation.

#include "udf.h"
#include "math.h"

/* DEFINES CONSTANTS FOR THE TURBULENCE MODELS */
const real sigmaK = 1;
const real sigmaEps = 1;
const real C1eps = 1.44;
const real C2eps = 1.92;


DEFINE_TURBULENT_VISCOSITY(user_mu_t, c, t)
{
real mu_t;
real rho = C_R(c,t);
real k = C_K(c,t);
real eps = C_D(c,t);
//real k = C_UDSI(c,t,0);
//real eps = C_UDSI(c,t,1);

mu_t = M_keCmu*rho*SQR(k)/eps;
//mu_t = 0.12*rho*SQR(k)/eps;

return mu_t;
}

DEFINE_DIFFUSIVITY(k_diff, c, t, eqn) // Diffusiviteten för k
{
real mu = C_MU_L(c,t);
real mu_t = C_MU_T(c,t);
real diff;

diff = mu+mu_t/sigmaK;

return diff;
}

DEFINE_DIFFUSIVITY(eps_diff, c, t, eqn) // Diffusiviteten för epsilon
{
real mu = C_MU_L(c,t);
real mu_t = C_MU_T(c,t);
real diff;

diff = mu+mu_t/sigmaEps;

return diff;
}

/* 2D-VERSION OF THE k SOURCE-TERM*/
DEFINE_SOURCE(k_source_twod, c, t, dS, eqn)
{
real rho = C_R(c,t);
real eps = C_D(c,t);
//real eps = C_UDSI(c,t,1);
real mu_t = C_MU_T(c,t);

real Gk = 2*mu_t*(C_DUDX(c,t)*C_DUDX(c,t) + 0.5*(C_DVDX(c,t)+C_DUDY(c,t))*(C_DVDX(c,t)+C_DUDY( c,t)) + C_DVDY(c,t)*C_DVDY(c,t) );

real source = Gk - rho*eps;
//dS[eqn] = 0.0;
return source;
}

/* 2D-VERSION OF THE EPSILON SOURCE-TERM*/
DEFINE_SOURCE(eps_source_twod, c, t, dS, eqn)
{
real k = C_K(c,t);
real eps = C_D(c,t);
//real k = C_UDSI(c,t,0);
//real eps = C_UDSI(c,t,1);
real mu_t = C_MU_T(c,t);
real rho = C_R(c,t);

real Gk = 2*mu_t*(C_DUDX(c,t)*C_DUDX(c,t) + 0.5*(C_DVDX(c,t)+C_DUDY(c,t))*(C_DVDX(c,t)+C_DUDY( c,t)) + C_DVDY(c,t)*C_DVDY(c,t) );

real source = C1eps*eps/k*Gk-C2eps*rho*eps*eps/k;
//dS[eqn] = 0.0;

return source;
}