CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

Laminar Kinetic Energy Model (Walters 2008)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 2 Post By logoswort

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 8, 2017, 10:33
Default Laminar Kinetic Energy Model (Walters 2008)
  #1
New Member
 
Join Date: May 2017
Posts: 2
Rep Power: 0
logoswort is on a distinguished road
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.
Attached Files
File Type: c k-kl-w.c (8.9 KB, 32 views)
grossz and dehai like this.
logoswort is offline   Reply With Quote

Old   May 9, 2017, 10:51
Default
  #2
New Member
 
Join Date: May 2017
Posts: 2
Rep Power: 0
logoswort is on a distinguished road
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;

}
logoswort is offline   Reply With Quote

Old   May 19, 2017, 20:41
Default
  #3
New Member
 
Mohammed
Join Date: Jan 2013
Posts: 7
Rep Power: 13
Madreed is on a distinguished road
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 ?
Madreed is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
question about turbulent kinetic energy junker4236 Main CFD Forum 19 April 19, 2017 05:46
Turbulent kinetic energy in k-epsilon model khunyeu FLUENT 0 October 7, 2013 02:31
Question about grids for laminar model and turbulent model Anna Tian Main CFD Forum 0 March 3, 2013 20:44
K - epsilon VS SST turbulence model Maicol Main CFD Forum 0 November 30, 2012 17:25
diff bet total energy model and thermal energy model?? vijeshjoshi23 Main CFD Forum 0 October 8, 2009 03:29


All times are GMT -4. The time now is 12:09.