Modified k-e turbulence model UDF
 

Hi,

I'm trying to implement a modified k-epsilon model written by Heschl, Sanz, and Klanatsky. My experience with UDFs is limited to a simple parabolic velocity inlet (the one in the FLUENT UDF User Guide,) so this complex turbulence model UDF is above my head. Can anyone tell me the process of implementing the UDF?

I can get as far as compiling the UDF. Beyond that, I'm lost.

Here's the entire UDF code:

/* turbulence model: nonlinear k-epsilon turbulence model by Ehrhard (1999) */ /* UDF (User Defined Function) Routine for FLUENT */ /* see also the paper: */ /* Implementation and comparison of different turbulence models for three dimensional wall jets with Fluent */ /* FLUENT CFD Forum 2005 Bad Nauheim, Germany */ /* Heschl Ch*., Sanz W.**, Klanatsky P.* */ /* * FachhochschulstudiengçÏge Burgenland GmbH, SteinamangerstraŽße 21, A-7423 Pinkafeld */ /* ** TU Graz, Institut of Thermal Turbomachinery and Machine Dynamics, Inffeldgasse 25, A-8010 Graz */ /* Contact: christian.heschl@fh-burgenland.at, wolfgang.sanz@tugraz.at, peter.klanatsky@fh-burgenland.at */

/* Die Autoren ?ernhemen keine Verantwortung oder Haftung f? den Inhalt der unten angef?rten UDF-Routinen */

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

/* Turbulence model constants */ const real C_1=-0.2; const real C_2=0.4; const real C_5=0.0; const real CE_1=1.44; const real CE_2=1.92; const real SIG_K=1.0; const real SIG_E=1.3; const real C_MU_ke=0.09; const real kappa=0.42; const real E_log=9.793;

/* User-defined scalars */ enum { UU, VV, WW, UV, UW, VW };

DEFINE_SOURCE(u_source, c, t, dS, eqn) { real source;

dS[eqn]=0.0; source= - C_R(c,t) * ( C_UDSI_G(c,t,UU)[0] + C_UDSI_G(c,t,UV)[1] + C_UDSI_G(c,t,UW)[2] );

return source; }

DEFINE_SOURCE(v_source, c, t, dS, eqn) { real source;

dS[eqn]= 0.0; source= - C_R(c,t) * ( C_UDSI_G(c,t,UV)[0] + C_UDSI_G(c,t,VV)[1] + C_UDSI_G(c,t,VW)[2] );

return source; }

DEFINE_SOURCE(w_source, c, t, dS, eqn) { real source;

dS[eqn]= 0.0; source= - C_R(c,t) * ( C_UDSI_G(c,t,UW)[0] + C_UDSI_G(c,t,VW)[1] + C_UDSI_G(c,t,WW)[2] );

return source; }

DEFINE_TURBULENT_VISCOSITY(user_mu_t,c,t) { real S,W,c_mu,T,mu_t; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33;

S11 = 0.5*( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = 0.5*( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = 0.5*( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = 0.5*( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = 0.5*( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = 0.5*( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = 0.5*( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = 0.5*( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = 0.5*( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = 0.5*( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = 0.5*( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = 0.5*( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = 0.5*( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = 0.5*( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = 0.5*( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = 0.5*( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = 0.5*( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = 0.5*( C_DWDZ(c,t)-C_DWDZ(c,t) );

T=C_K(c,t)/C_D(c,t);

S=T*sqrt( 2*(S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23 +S31*S31+S32*S32+S33*S33) ); W=T*sqrt( 2*(W11*W11+W12*W12+W13*W13+W21*W21+W22*W22+W23*W23 +W31*W31+W32*W32+W33*W33) );

c_mu=MIN( 1./(0.9*pow(S,1.4)+0.4*pow(W,1.4)+3.5) , 0.15 );

mu_t=c_mu*T*C_K(c,t);

return mu_t; }

DEFINE_ADJUST(rsm_adjust,domain) {

Thread *t; cell_t c; real T; real mu_t, c_mu; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real S, W, Sij, Wij; real C_3, C_4, C_6, C_7; real P11, P22, P33, Pkk; real tau_w, u_tauw, u_star, y_star, Gk, u_mag;

/* Set the turbulent viscosity */ thread_loop_c(t,domain) if (FLUID_THREAD_P(t)) { begin_c_loop(c,t) {

S11 = 0.5*( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = 0.5*( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = 0.5*( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = 0.5*( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = 0.5*( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = 0.5*( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = 0.5*( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = 0.5*( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = 0.5*( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = 0.5*( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = 0.5*( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = 0.5*( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = 0.5*( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = 0.5*( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = 0.5*( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = 0.5*( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = 0.5*( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = 0.5*( C_DWDZ(c,t)-C_DWDZ(c,t) );

T=C_K(c,t)/C_D(c,t);

S=T*sqrt( 2*(S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23 +S31*S31+S32*S32+S33*S33) ); W=T*sqrt( 2*(W11*W11+W12*W12+W13*W13+W21*W21+W22*W22+W23*W23 +W31*W31+W32*W32+W33*W33) );

Sij=S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S2 3+S31*S31+S32*S32+S33*S33; Wij=W11*W11+W12*W12+W13*W13+W21*W21+W22*W22+W23*W2 3+W31*W31+W32*W32+W33*W33;

c_mu=MIN( 1./(0.9*pow(S,1.4)+0.4*pow(W,1.4)+3.5) , 0.15 );

mu_t=c_mu*T*C_K(c,t);

C_3=2.0-exp(-SQR(S-W)); C_4=-32.0*SQR(c_mu); C_6=-16.0*SQR(c_mu); C_7=16.0*SQR(c_mu);

C_UDSI(c,t,UU)= C_1*mu_t*T* (S11*S11+S12*S12+S13*S13-1./3.*Sij) + 2.*C_2*mu_t*T*(W12*S12+W13*S13) + C_3*mu_t*T*(W12*W12+W13*W13-1./3.*Wij) - 2.*C_4*mu_t*T*T * (W12*(S11*S12+S12*S22+S13*S23)+W13*(S11*S13+S12*S2 3+S13*S33)) + C_6*mu_t*T*T* (S11*Sij) + C_7*mu_t*T*T* (S11*Wij);

C_UDSI(c,t,VV)= C_1*mu_t*T* (S12*S12+S22*S22+S23*S23-1./3.*Sij) + 2.*C_2*mu_t*T*(W23*S23-W12*S12) + C_3*mu_t*T*(W12*W12+W23*W23-1./3.*Wij) + 2.*C_4*mu_t*T*T * (W12*(S11*S12+S12*S22+S13*S23)-W23*(S12*S13+S22*S23+S23*S33)) + C_6*mu_t*T*T* (S22*Sij) + C_7*mu_t*T*T* (S22*Wij);

C_UDSI(c,t,WW)= C_1*mu_t*T* (S13*S13+S23*S23+S33*S33-1./3.*Sij) - 2.*C_2*mu_t*T*(W13*S13+W23*S23) + C_3*mu_t*T*(W13*W13+W23*W23-1./3.*Wij) + 2.*C_4*mu_t*T*T * (W13*(S11*S13+S12*S23+S13*S33)+W23*(S12*S13+S22*S2 3+S23*S33)) + C_6*mu_t*T*T* (S33*Sij) + C_7*mu_t*T*T* (S33*Wij);

C_UDSI(c,t,UV)= C_1*mu_t*T* (S11*S12+S12*S22+S13*S23) + C_2*mu_t*T*(W12*(S22-S11)+W13*S23+W23*S13) + C_3*mu_t*T*(W13*W23) + C_4*mu_t*T*T* (W12*(S11*S11+S13*S13-S22*S22-S23*S23)-W13*(S12*S13+S22*S23+S23*S33)-W23*(S11*S13+S12*S23+S13*S33)) + C_6*mu_t*T*T* (S12*Sij) + C_7*mu_t*T*T* (S12*Wij);

C_UDSI(c,t,UW)= C_1*mu_t*T* (S11*S13+S12*S23+S13*S33) + C_2*mu_t*T*(W13*(S33-S11)+W12*S23-W23*S12) - C_3*mu_t*T*(W12*W23) + C_4*mu_t*T*T* (W13*(S11*S11+S12*S12-S23*S23-S33*S33)-W12*(S11*S13+S22*S23+S23*S33)+W23*(S11*S12+S12*S22 +S13*S23)) + C_6*mu_t*T*T* (S13*Sij) + C_7*mu_t*T*T* (S13*Wij);

C_UDSI(c,t,VW)= C_1*mu_t*T* (S12*S13+S22*S23+S23*S33) + C_2*mu_t*T*(W23*(S33-S22)-W12*S13-W13*S12) + C_3*mu_t*T*(W12*W13) + C_4*mu_t*T*T* (W23*(S12*S12+S22*S22-S13*S13-S33*S33)+W12*(S11*S13+S12*S23+S13*S33)+W13*(S11*S1 2+S12*S22+S13*S23)) + C_6*mu_t*T*T* (S23*Sij) + C_7*mu_t*T*T* (S23*Wij);

C_UDMI(c,t,0)= -2.*mu_t*S11 + 2./3.*C_K(c,t) + C_UDSI(c,t,UU); C_UDMI(c,t,1)= -2.*mu_t*S22 + 2./3.*C_K(c,t) + C_UDSI(c,t,VV); C_UDMI(c,t,2)= -2.*mu_t*S33 + 2./3.*C_K(c,t) + C_UDSI(c,t,WW); C_UDMI(c,t,3)= -2.*mu_t*S12 + C_UDSI(c,t,UV); C_UDMI(c,t,4)= -2.*mu_t*S13 + C_UDSI(c,t,UW); C_UDMI(c,t,5)= -2.*mu_t*S23 + C_UDSI(c,t,VW);

} end_c_loop(c,t) } }

Re: Modified k-e turbulence model UDF - help!
 

Hi, As I may underestand your question "you want to implement this source code to FLUENT", After compiling you should go to "Define\BC\fluid and set the source terms of Reynolds stresses " and add respect source to transport equation. Hope it could help.

Re: Modified k-e turbulence model UDF - help!
 

I've managed to compile the BCs, hook the ADAPT function to the model, define the turbulent viscosity via DEFINE-MODELS-VISCOUS, and define the momentum terms via DEFINE-BC-FLUID.

3 questions: first, the u_source, v_source, and w_source terms are the momentum terms, right? Looks like it to me but I wanted to be sure.

Second, the code is written for a 3d case. Can it be applied to a 2d case, or does the code have to be modified to remove the third dimension?

Third, I've applied the BCs i listed above, and I get a SEGMENTATION ERROR when I try to initialize. Or, if I use an already converged solution, I get the same SEGMENTATION ERROR when I try to iterate. Any thoughts on why?

Thanks! I know that's a lot of questions.

Re: Modified k-e turbulence model UDF - help!
 

Hi, First: They are the Reynolds Stress Source Terms. Seconde: Yes I think you could, just apply dw/dz=0 and w=0; Third: Second will solve this problem.

Hope this helps!

Hi,
Can you successfully implement this turbulence model in Fluent?

Hi Travis,

did you manage to run this code successfully?
Because I have the same problem
can you help me?
thanks in advance!
best!

hi
 

does anybody run this model ? is this model working fine. please share the updates. i want to use this code with small modifications.

thanks

Hi,
I'v just define a low Reynolds turbulence model with UDS, but it didn't work well. Have you ever run the udf above?