[logoswort]

Hello Community,

for my master thesis, I want to program my own transition model with UDF Fluent and run it in a parallel compiled process.
Attached is the .c file of my model. It is build regarding the k-kL-w model from Walters in 2008 in JoF. However, so far the model is not reproducing the results as in the publications regarding the flat plate test cases ,T3x.

source code errors:
- momentum source terms x,y,and z for the Boussinesq modeling are not written correctly
-the logic of the model is not correct

hooking errors:
- sources, diffusivity, turbulent viscosity, equations are not correctly set and hooked

Could you give me please some advice regarding my model and the implementation into Fluent. Later on, I want to parallize the model in order to calculate bigger 3D meshes.

[logoswort]

I hooked this model as follows to Fluent:

- Set Models-Viscous- SST k-omega
- Set Models-Viscous - User-Defined Functions: Turbulent Viscosity (mod_mu_t)
- Set Materials-Fluid-Air USD Diffusivity user-defined user_diffusivity
- Set Cell Zone Conditions - fluid(fluid) - Source Terms - Add X-,Y-,Z- Momentum Equation (udf x,y,z_source), UDS Scalar 0,1,2 for udf k,KL,omega_source
- Set Boundary Conditions :
1.inlet(velocity-inlet)
2.outlet (outflow)
3.wall(wall)
- Solution - Solution Controls - Equations- Turbulence (Deactivate) User Scalar 0-2 (Activate)

/* Transitional k-kl-w model, based on Walters, JoFE (2008) */

#include "udf.h"
#include "turb.h"
#include "mem.h"

/*USER DEFINED SCALAR*/
#define K 0
#define KL 1
#define W 2
#define SRT_K 3
#define SRT_KL 4
#define nuTldudydvdx 5
#define nuTldwdxdudz 6
#define nuTldwdydvdz 7
#define nuTldudx 8
#define nuTldvdy 9
#define nuTldwdz 10

/*USER DEFINED MEMORY*/
/*#define K 0
#define KL 1
#define W 2*/
#define MU_T 3
#define ALPHA_T 4
#define DSRTK_DX 5
#define DSRTKL_DX 6
#define PROD_K 7
#define PROD_KL 8
#define R_BP 9
#define R_NAT 10
#define F_WALL 11
#define F_OMEGA 12
#define ROT_RATE 13
#define RE_ROT 14
#define PHI_NAT 15

/* MODEL CONSTANTS */
#define A_0 4.04
#define A_BP 0.6
#define A_NAT 200.
#define A_S 2.12
#define A_TS 200.
#define A_NU 6.75
#define C_BP_CRIT 1.2
#define C_L1 3.4e-6
#define C_L2 1.e-10
#define C_INT 0.75
#define C_NAT_CRIT 1250.0
#define C_NC 0.1
#define C_R1 0.12
#define C_R_NAT 0.02
#define C_SS 1.5
#define C_TS_CRIT 1000.0
#define C_W1 0.44
#define C_W2 0.92
#define C_W3 0.3
#define C_WR 1.5
#define C_TAU 4360.0
#define C_MU_STD 0.09
#define C_L 2.495
#define SIG_K 1.0
#define SIG_W 1.17
#define TINY 1.e-12

DEFINE_TURBULENT_VISCOSITY(mod_mu_t, c, t)
{
return C_UDMI(c,t,MU_T);
}

DEFINE_DIFFUSIVITY(user_diffusivity, c, t, eqn)
{

real diff;
switch(eqn)
{
case K:
diff = C_MU_L(c,t) + C_R(c,t)*C_UDMI(c,t,ALPHA_T)/SIG_K;
break;
case KL:
diff = C_MU_L(c,t);
break;
case W:
diff= C_MU_L(c,t) + C_R(c,t)*C_UDMI(c,t,ALPHA_T)/SIG_W;
break;
default:
diff = C_MU_L(c,t) + C_UDMI(c,t,MU_T);
}
return diff;
}

DEFINE_ADJUST(adjust_fn_kw_sst, domain)
{

Thread *t;
cell_t c;

real nu;
real l_turb, l_eff;
real rot_rate_ten;
real re_t, f_nu, f_ss, f_int, c_mu, f_omega;
real k_s, k_l, mu_t_s, mu_t_l, mu_t_l0, mu_t_l1;
real d_eff, b_ts, f_tau;
real phi_bp, beta_bp, f_nat_crit, phi_nat, beta_nat;

thread_loop_c(t,domain)
{
if (&C_UDSI_G(0,t,K)[0] != NULL) /* UDS Gradient von k */
{
begin_c_loop(c,t)
{
nu = C_MU_L(c,t)/C_R(c,t); /* Kinematic viscosity */

/* Limit turbulence variables to be non-zero for stability */
C_UDMI(c,t,K) = MAX(C_UDSI(c,t,K),0.5*(C_UDMI(c,t,K)+1.e-18));
C_UDMI(c,t,KL) = MAX(C_UDSI(c,t,KL),0.5*(C_UDMI(c,t,KL)+1.e-18));
C_UDMI(c,t,W) = MAX(C_UDSI(c,t,W),0.5*(C_UDMI(c,t,W)+1.e-18));

C_K(c,t) = C_UDMI(c,t,K);

/* Compute near wall destruction terms */
C_UDSI(c,t,SRT_K) = sqrt(C_UDMI(c,t,K));
C_UDSI(c,t,SRT_KL) = sqrt(C_UDMI(c,t,KL));
C_UDMI(c,t,DSRTK_DX) = C_UDSI_G(c,t,SRT_K)[0]*C_UDSI_G(c,t,SRT_K)[0]+C_UDSI_G(c,t,SRT_K)[1]*C_UDSI_G(c,t,SRT_K)[1];
C_UDMI(c,t,DSRTKL_DX) = C_UDSI_G(c,t,SRT_KL)[0]*C_UDSI_G(c,t,SRT_KL)[0]+C_UDSI_G(c,t,SRT_KL)[1]*C_UDSI_G(c,t,SRT_KL)[1];
/*#if RP_3D*/
C_UDMI(c,t,DSRTK_DX) += C_UDSI_G(c,t,SRT_K)[2]*C_UDSI_G(c,t,SRT_K)[2];
C_UDMI(c,t,DSRTKL_DX) += C_UDSI_G(c,t,SRT_KL)[2]*C_UDSI_G(c,t,SRT_KL)[2];
/*#endif*/

/* Turbulence length scales, effective (wall-limited) turbulence length scale and kinematic wall damping function */
l_turb = sqrt(C_UDMI(c,t,K))/C_UDMI(c,t,W);
l_eff = MIN(C_L*C_WALL_DIST(c,t),l_turb);
C_UDMI(c,t,F_WALL) = l_eff/l_turb; /*Mistake, correction Fuerst(2012): pow(l_eff/l_turb,0.666667);*/

/*Vorticity Omega: C.KATENBRING*/
rot_rate_ten=((C_DWDY(c,t)-C_DVDZ(c,t))+(C_DUDZ(c,t)-C_DWDX(c,t))+(C_DVDX(c,t)-C_DUDY(c,t))+TINY);
C_UDMI(c,t,ROT_RATE)=sqrt(pow(rot_rate_ten,2));

/* Damping functions */
re_t = pow(C_UDMI(c,t,F_WALL),2.0)*C_UDMI(c,t,K)/(nu*C_UDMI(c,t,W));
f_nu = 1.0 - exp(-sqrt(re_t)/A_NU); /* Viscous damping */
f_ss = exp(-pow(C_SS*nu*C_UDMI(c,t,ROT_RATE)/C_UDMI(c,t,K),2.0)); /* Shear sheltering*/
f_int = MIN(C_UDMI(c,t,KL)/(C_INT*(C_UDMI(c,t,K)+C_UDMI(c,t,KL))),1.0); /* Intermittency */
f_omega = 1.0 - exp(-0.41*pow(l_eff/l_turb,4));

/* Turbulent viscosity coefficient */
c_mu = 1.0/(A_0+A_S*(C_STRAIN_RATE_MAG(c,t)/C_UDMI(c,t,W)));

/* Small-scale and large-scale turbulence contributions */
k_s = f_ss*C_UDMI(c,t,F_WALL)*C_UDMI(c,t,K); /* Effective Small-Scale Turbulence*/
k_l = C_UDMI(c,t,K) - k_s;

/* Small-scale effective turbulent viscosity */
mu_t_s = C_R(c,t)*f_nu*f_int*c_mu*sqrt(k_s)*l_eff; /*contain rho already*/

/* Large-scale effective turbulent viscosity */
C_UDMI(c,t,RE_ROT) = C_WALL_DIST(c,t)*C_WALL_DIST(c,t)*C_UDMI(c,t,ROT_R ATE)/nu;
b_ts = 1.0 - exp(-pow(MAX(C_UDMI(c,t,RE_ROT) - C_TS_CRIT,0.0),2.0)/A_TS);
f_tau = 1.0 - exp(-C_TAU*k_l/(pow(l_eff*C_UDMI(c,t,ROT_RATE),2.0)+TINY));

mu_t_l0 = C_R(c,t)*C_L1*f_tau*C_UDMI(c,t,ROT_RATE)*pow(l_eff ,3.0)*sqrt(k_l)/nu;
mu_t_l1 = 0.5* C_R(c,t)*(MAX(C_UDMI(c,t,KL)+k_l,0))/(C_STRAIN_RATE_MAG(c,t)+TINY); /*, MAX and TINY for stability*/
mu_t_l = MIN(mu_t_l0+C_R(c,t)*b_ts*C_L2*C_UDMI(c,t,RE_ROT)* pow(C_WALL_DIST(c,t),2.0)*C_UDMI(c,t,ROT_RATE),mu_ t_l1);/*contains rho*/

/* Scalar diffusivity */
C_UDMI(c,t,ALPHA_T) = f_nu*C_MU_STD*sqrt(k_s)*l_eff;

/* Transition source terms */

/* Bypass */
phi_bp = MAX(C_UDMI(c,t,K)/(C_UDMI(c,t,ROT_RATE)*nu+TINY) - C_BP_CRIT,0.0);
beta_bp = 1.0 - exp(-phi_bp/A_BP);
C_UDMI(c,t,R_BP) = C_R1*beta_bp*MAX(C_UDMI(c,t,KL),0.0)*C_UDMI(c,t,W)/C_UDMI(c,t,F_WALL);/*Max - Stability*/
/* Natural */
f_nat_crit = 1.0 - exp(C_NC*sqrt(C_UDMI(c,t,KL))*C_WALL_DIST(c,t)/nu);/*Negative sign in JoF:Walters 2008*/
C_UDMI(c,t,PHI_NAT) = MAX(C_UDMI(c,t,RE_ROT)-C_NAT_CRIT/(f_nat_crit+TINY),0.0);
beta_nat = 1.0 - exp(-C_UDMI(c,t,PHI_NAT)/A_NAT);
C_UDMI(c,t,R_NAT) = C_R_NAT*beta_nat*MAX(C_UDMI(c,t,KL),0.0)*C_UDMI(c, t,ROT_RATE);/*Max - Stability*/

/* Compute turbulent viscosity */
C_UDMI(c,t,MU_T) = mu_t_s + mu_t_l;

/* Adjust derivative */
C_UDSI(c,t,nuTldudx) = C_UDMI(c,t,MU_T)*C_DUDX(c,t)*2;
C_UDSI(c,t,nuTldvdy) = C_UDMI(c,t,MU_T)*C_DVDY(c,t)*2;
C_UDSI(c,t,nuTldwdz) = C_UDMI(c,t,MU_T)*C_DWDZ(c,t)*2;
C_UDSI(c,t,nuTldudydvdx) = C_UDMI(c,t,MU_T)*(C_DUDY(c,t)+C_DVDX(c,t));
C_UDSI(c,t,nuTldwdxdudz) = C_UDMI(c,t,MU_T)*(C_DWDX(c,t)+C_DUDZ(c,t));
C_UDSI(c,t,nuTldwdydvdz) = C_UDMI(c,t,MU_T)*(C_DWDY(c,t)+C_DVDZ(c,t));

/* Compute turbulence production terms */
C_UDMI(c,t,PROD_KL) = mu_t_l*C_STRAIN_RATE_MAG(c,t)*C_STRAIN_RATE_MAG(c, t);
C_UDMI(c,t,PROD_K) = mu_t_s*C_STRAIN_RATE_MAG(c,t)*C_STRAIN_RATE_MAG(c, t);

/* Kinematic damping function */
C_UDMI(c,t,F_OMEGA) = f_omega;
}
end_c_loop(c,t)
}
else
{
begin_c_loop(c,t)
{
C_UDMI(c,t,K) = MAX(C_UDSI(c,t,K),1.e-16);
C_UDMI(c,t,KL) = MAX(C_UDSI(c,t,KL),1.e-16); /*Maybe change*/
C_UDMI(c,t,W) = MAX(C_UDSI(c,t,W),0.1*(C_MU_L(c,t)/C_R(c,t))/(C_WALL_DIST(c,t)*C_WALL_DIST(c,t)));

C_UDMI(c,t,MU_T) = C_R(c,t)*C_UDMI(c,t,K)/C_UDMI(c,t,W);
}
end_c_loop(c,t)
}
}
}

DEFINE_SOURCE(x_source,c,t,dS,eqn)
{
real x[ND_ND];
real source;
real k_TOT;

source = C_R(c,t)*((-2.0/3.0)*(C_UDSI_G(c,t,K)[0]+C_UDSI_G(c,t,KL)[0])+(C_UDSI_G(c,t,nuTldudx)[0]+C_UDSI_G(c,t,nuTldudydvdx)[1]+C_UDSI_G(c,t,nuTldwdxdudz)[2]+TINY));

return source;
}

DEFINE_SOURCE(y_source,c,t,dS,eqn)
{
real x[ND_ND];
real source;

source = C_R(c,t)*((-2.0/3.0)*(C_UDSI_G(c,t,K)[1]+C_UDSI_G(c,t,KL)[1])+(C_UDSI_G(c,t,nuTldudydvdx)[0]+C_UDSI_G(c,t,nuTldvdy)[1]+C_UDSI_G(c,t,nuTldwdydvdz)[2]+TINY));

return source;
}

DEFINE_SOURCE(z_source,c,t,dS,eqn)
{
real x[ND_ND];
real source;

source = C_R(c,t)*((-2.0/3.0)*(C_UDSI_G(c,t,K)[2]+C_UDSI_G(c,t,KL)[2])+(C_UDSI_G(c,t,nuTldwdxdudz)[0]+C_UDSI_G(c,t,nuTldwdydvdz)[1]+C_UDSI_G(c,t,nuTldwdz)[2]+TINY));

return source;
}

DEFINE_SOURCE(k_source, c, t, dS, eqn)
{

real Source;

Source = C_UDMI(c,t,PROD_K); /*Prod_K contain rho already*/
Source += C_R(c,t)*(C_UDMI(c,t,R_BP)+C_UDMI(c,t,R_NAT));
Source -= C_R(c,t)*C_UDMI(c,t,W)*C_UDSI(c,t,K);
Source -= C_R(c,t)*2.0*C_MU_L(c,t)*C_UDMI(c,t,DSRTK_DX)*(C_U DSI(c,t,K)/C_UDMI(c,t,K));/*Last fraction for derivative*/

dS[eqn] = - C_R(c,t)*C_UDMI(c,t,W) - C_R(c,t)*2.0*C_MU_L(c,t)*C_UDMI(c,t,DSRTK_DX)*(1.0/C_UDMI(c,t,K));

return Source;

}

DEFINE_SOURCE(KL_source, c, t, dS, eqn)
{

real Source;

Source = C_UDMI(c,t,PROD_KL); /*Prod_K contain rho already*/
Source -= C_R(c,t)*(C_UDMI(c,t,R_BP)+C_UDMI(c,t,R_NAT));
Source -= C_R(c,t)*2.0*C_MU_L(c,t)*C_UDMI(c,t,DSRTKL_DX)*(C_ UDSI(c,t,KL)/C_UDMI(c,t,KL));

dS[eqn] = -C_R(c,t)*2.0*C_MU_L(c,t)*C_UDMI(c,t,DSRTKL_DX)*(1. 0/C_UDMI(c,t,KL));

return Source;

}

DEFINE_SOURCE(omega_source, c, t, dS, eqn)
{

real Source;

Source = C_W1*C_UDMI(c,t,PROD_K)*C_UDSI(c,t,W)/C_UDMI(c,t,K); /*Prod_K contain rho already*/
Source += C_R(c,t)*(C_WR/C_UDMI(c,t,F_WALL)-1.0)*(C_UDMI(c,t,R_BP)+C_UDMI(c,t,R_NAT))*C_UDSI(c ,t,W)/C_UDMI(c,t,K);
Source -= C_R(c,t)*C_W2*C_UDMI(c,t,W)*C_UDSI(c,t,W);
Source += C_R(c,t)*C_W3*C_UDMI(c,t,F_OMEGA)*C_UDMI(c,t,ALPHA _T)*pow(C_UDMI(c,t,F_WALL),2)*sqrt(C_UDMI(c,t,K))/pow(C_WALL_DIST(c,t),3);

dS[eqn] = C_W1*C_UDMI(c,t,PROD_K)*1.0/C_UDMI(c,t,K) + C_R(c,t)*(C_WR/C_UDMI(c,t,F_WALL)-1.0)*(C_UDMI(c,t,R_BP)+C_UDMI(c,t,R_NAT))*1.0/C_UDMI(c,t,K) - C_R(c,t)*C_W2*C_UDMI(c,t,W);

return Source;

}

[Madreed]

I want to ask you a question about hooking the diffusivity
when I load the UDF, the default diffusion is set to (defined per uds)
I click edit and change it to user defined and choose the diffusion equation ,but unfortunately when I initialize an error occurs and the fluent closes , How could you initialize your case?
am I missing something ?