I don't think/know the production term and dissipation term can be replaced individually. This is also what I want to do before. I think this is one short coming of Fluent, i.e. it lets you add source term, but it doesn't allow you modify source term that have already been there. You probably need to define a new turbulence model. This may introduce a lot of work, although we just want to have a small change. If Fluent can hear me, I strongly suggest Fluent make it easier for the user to implement turbulence models. I have seen some researcher publish their work using their own turbulence model in Fluent, but I still don't manage to do it easily. When I get all the equations right in the code, the solver get divergence and I get frustrated and then give up. I may try it out later after my graduation, but at that time I may switch to other open source code as full control can be available from there. Good luck with your work!

Thank you very much for your answer. You are right, we can not play with source term that have already been in Fluent and it is a gap that could be covered in next version of the software.
Any way, you said that your implement model has not been converged. I do not know you change the under relaxation of your new trasport equation or not. You may control it by decreasing this value and it directly has affected on your convergence, too.
Moreover, I read your code in this thread, I think that it is beter to use DEFINE_TURBULENT_VISCOSITY instead of DEFINE_ADJUST for introducing eddy viscosity.
It would be appricated if you tell me more about your experiances. Now I just want to write k-w SST model in UDF and after that I got the same results with fluent, then I will change it.

hi,
did you manage to get results from this UDF code?
particularly I am interested in tensor definitions:
for example:"S[m][i]=0.5*(puixj[m][i]+puixj[i][m]); //S is the stress tensor”
is it getting the right values of derivatives if you define like puixj[m][i] and puixj[i][m] derivatives they should give different values…in my case they are giving the same values…
thanks for help in advance!

Hi Kanarya,
Sorry I didn't understand. Can you illustrate better?

Regards,

Hi,
thanks for the quick answer!
did you manage to get results from your code?
and it seems I am getting same values for puixj[i][m] and puixj[m][i] although they are with opposite indices?
do you have any idea why?
thanks!

Yes and no. I manage to make it running. But I can't get it converged. Later I give up.
For that puixj[i][m], following the definition, they are derivatives so they are not necessarily the same? For example, partial u / partial y and partial v / partial x, they don't need to be the same?
I don't think you need to worry about this, as the code gets derivatives from FLUENT. If there is something wrong, it is due to FLUENT, rather than the code. And if that is the case, you need to calculate the derivative yourself instead of using FLUENT, i.e., write a part of code instead of using the built-in FLUENT macro.

Regards,

hi,
thanks!
I have convergence problem as well....do you thing i should initialize the tensors?
thanks!

Hi Kanarya,
I have got a lot of work to do now. I can't provide a clear answer to you at this moment. I also tried to seek solutions for that convergence problem. What I can suggest is, if you have sufficient time, consider other opensource code you have access to. I am not saying FLUENT is not a good code. It is good but we know nothing about it, so when we have problems we have limited chance to know what is actually going on. On the contrary, with a code either your research group used for a long time or else, you have full control on any part of the code, so that investigating convergence and underlying mistakes in your implementation becomes possible. If you don't have time, one thing you can do might define a turbulence model that is available inside FLUENT first, so then you can at least compare the one you define and the one already in to see what is suspecting to be wrong. Boundary condition and other stuff, everything matters!

Regards,

Sst udf
 

Hello all,
Does anyone write an UDF for SST model? Did you get the same results that Fluent gives us?

Actually, I am writing an UDF for a new turbulence model that is very close to the SST model. Therefore, I decided that write SST model by UDF, firstly. To make sure the UDF working correctly, I want to compare the result of UDF and Fluent.

Hi ,
If you use the fluent theory guide, you should get similar results...

Hello,
Thank you very much for your answer. Did you try to write an UDF for SST? Which part of Fluent theory guide is your aim?

Actually, I am writting an UDF for SST model, my problem is the boundary condition and the value of them is not the same with Fluent one.

Hi,
I did several turbulence models but not specifically SST and I think it does not matter that much!
I think you should describe your problem more in detail...
And then you could receive some help...
thanks!

Thank you very much for your answer. My goal is to write an UDF for SST model and solve it for a 2D flat plate case.
I run it but it gives me a wrong answer and the result is not the same with Fluent.

this is my code:
/**This UDF is k-w SST model without Low Reynolds correction and any other options,just SST equation is applied**/

#include "udf.h"
#include "mem.h"
#include "math.h"
#include "sg_udms.h"
#include "global.h"

/* User-defined constants */
#define SIG_TKE1 1.176
#define SIG_OMG1 2.0
#define SIG_TKE2 1.0
#define SIG_OMG2 1.168
#define ALPHA_1 0.31
#define ALPHA_infstar 1.0
#define ALPHA_star 1.0
#define ALPHA_0 1/9
#define BETA_i1 0.075
#define BETA_i2 0.0828
#define BETA_infstar 0.09 // the compressibility correction is not enabled for k
#define R_BETA 8.0
#define R_TKE 6.0
#define R_OMG 2.95
#define ZETA_star 1.5
#define M_t0 0.25
#define K_karman 0.41

#define MAX(x,y) (((x) > (y)) ? (x) : (y))
#define MIN(x,y) (((x) < (y)) ? (x) : (y))

/* User-defined scalars */
enum
{
TKE,
OMG,
N_REQUIRED_UDS
};

/*==========================define dissipation term of k (Y_K)============================*/
real Y_k (cell_t c,Thread *t)
{return C_R(c,t)*BETA_infstar*C_UDSI(c,t,TKE)*C_UDSI(c,t,O MG);}

/*======================define production of k(G_Kbar) for each cell=====================*/
real G_kbar(cell_t c,Thread *t)
{
real G_k,G_bar,s1; // refers to the production of k
G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t)); // Bounsinesq hypothesis, production of k
s1 = 10*Y_k(c,t);

G_bar = MIN(G_k,s1);
return G_bar;
}
/**********************************Source term of TKE************************************/
DEFINE_SOURCE(k_source, c, t, dS, eqn)
{
real dGK_bar,dY_k;

dY_k = C_R(c,t)*BETA_infstar*C_UDSI(c,t,OMG);
dGK_bar = MIN((10*dY_k),0);

/* dS[eqn] = dGK_bar-dY_k; */
dS[eqn] = -C_R(c,t)*BETA_infstar*C_UDSI(c,t,OMG);
return G_kbar(c,t)-Y_k(c,t);
}
/*==========================================F1===== =====================================*/
real F1(cell_t c,Thread *t)
{
real s1,s2,s3,s4,cdkw1,arg1,cdkw,cdkw2=pow(10,-10);

s1 = sqrt(C_UDSI(c,t,TKE))/(0.09*C_WALL_DIST(c,t)*MAX(C_UDSI(c,t,OMG),1.e-10));
s2 = 500*C_MU_L(c,t)/(C_R(c,t)*SQR(C_WALL_DIST(c,t))*MAX(C_UDSI(c,t,OMG ),1.e-10));
s3 = MAX(s1,s2);

/* compute CDkw */
cdkw1 = (2*C_R(c,t)*NV_DOT(C_UDSI_G(c,t,TKE),C_UDSI_G(c,t, OMG)))/(SIG_OMG2*MAX(C_UDSI(c,t,OMG),1.e-10));
cdkw = MAX(cdkw1,cdkw2);

s4 = 4*C_R(c,t)*C_UDSI(c,t,TKE)/(SIG_OMG2*cdkw*SQR(C_WALL_DIST(c,t)));
arg1 = MIN(s4, s3);
return tanh(arg1*arg1*arg1*arg1);
}
/**********************************Source term of OMG************************************/
DEFINE_SOURCE(omg_source, c, t, dS, eqn)
{
real ALPHA_inf,ALPHA_inf1,ALPHA_inf2,BETA,dG_w,dY_w;

ALPHA_inf1 = (BETA_i1/BETA_infstar)-(SQR(K_karman)/(SIG_OMG1*sqrt(BETA_infstar)));
ALPHA_inf2 = (BETA_i2/BETA_infstar)-(SQR(K_karman)/(SIG_OMG2*sqrt(BETA_infstar)));
/* ALPHA_inf = F1(c,t)*ALPHA_inf1+(1-F1(c,t))*ALPHA_inf2; */
ALPHA_inf = 0.52;
BETA = F1(c,t)*BETA_i1+(1-F1(c,t))*BETA_i2;

dY_w = 10*C_R(c,t)*BETA_infstar*C_UDSI(c,t,TKE);
dG_w = MIN(dY_w,0);

dS[eqn]=-(fabs(2*(1-F1(c,t))*C_R(c,t)*NV_DOT(C_UDSI_G(c,t,TKE),C_UDSI_ G(c,t,OMG))/(MAX(C_UDSI(c,t,OMG),1.e-10)*SIG_OMG2))+2*C_R(c,t)*BETA*C_UDSI(c,t,OMG))/MAX(C_UDSI(c,t,OMG),1.e-10);
/* ALPHA_inf*dG_w*C_R(c,t)/C_MU_T(c,t); */

return ALPHA_inf*G_kbar(c,t)*C_R(c,t)/C_MU_T(c,t)-C_R(c,t)*BETA*SQR(C_UDSI(c,t,OMG))+2*(1-F1(c,t))*C_R(c,t)*NV_DOT(C_UDSI_G(c,t,TKE),C_UDSI_ G(c,t,OMG))/(MAX(C_UDSI(c,t,OMG),1.e-10)*SIG_OMG2);
}
/*==========================================F2===== =====================================*/
real F2(cell_t c,Thread *t)
{
real s1,s2,arg2;

s1 = 2*sqrt(C_UDSI(c,t,TKE))/(0.09*C_WALL_DIST(c,t)*MAX(C_UDSI(c,t,OMG),1.e-10));
s2 = 500*C_MU_L(c,t)/(C_R(c,t)*SQR(C_WALL_DIST(c,t))*MAX(C_UDSI(c,t,OMG ),1.e-10));
arg2 = MAX(s1,s2);

return tanh(arg2*arg2);
}
/*************************************eddy viscosity**************************************/
DEFINE_TURBULENT_VISCOSITY(user_mu_t,c,t)
{
real mu_t,s1,s2;

s1=Strainrate_Mag(c,t)*F2(c,t)/(ALPHA_1*MAX(C_UDSI(c,t,OMG),1.e-10));
s2=1/ALPHA_star;

mu_t=C_R(c,t)*C_UDSI(c,t,TKE)/(MAX(s1,s2)*MAX(C_UDSI(c,t,OMG),1.e-10));
return mu_t;
}
/*================================Prandtl number for TKE==================================*/
real SIG_TKE(cell_t c,Thread *t)
{ return 1/((F1(c,t)/SIG_TKE1)+((1-F1(c,t))/SIG_TKE2));}

/*================================Prandtl number for OMG==================================*/
real SIG_OMG(cell_t c,Thread *t)
{ return 1/((F1(c,t)/SIG_OMG1)+((1-F1(c,t))/SIG_OMG2));}

/***************************Diffusivity term of TKE and OMG********************************/
DEFINE_DIFFUSIVITY(kw_diff, c, t, eqn)
{
real diff; // define the diffusion coeffcient
switch (eqn)
{
case TKE:
diff=C_MU_T(c,t)/SIG_TKE(c,t)+C_MU_L(c,t);
break;
case OMG:
diff=C_MU_T(c,t)/SIG_OMG(c,t)+C_MU_L(c,t);
break;
default:
diff=C_MU_T(c,t)+C_MU_L(c,t);
}
return diff;
}
/*===================Wall boundary=======================*/
DEFINE_PROFILE(wall_d_bc,t,i)
{
Thread *t0;
face_t f;
cell_t c, c0;
real yplus,wshear;

begin_f_loop(f,t)
{
c0 = F_C0(f,t);
t0 = THREAD_T0(t);

yplus = C_STORAGE_R(f,t,SV_WALL_YPLUS_UTAU); /* Y+*/
wshear = C_MU_L(c0,t0)*C_U(c0,t0)/C_WALL_DIST(c0,t0); /* wall shear stress */
F_PROFILE(f,t,i) = 3.*wshear/(BETA_infstar*C_MU_L(c0,t0)*SQR(yplus));
}
end_f_loop(f,t)
}

You can use this code for other 2D problems and check that it does not give us the correct results.

Hi,

I can not debug your code but I can give suggestions...do you have some results which you can show here...
thanks!

Actually, I have not gotten a good results yet and must of time the simulation would be divege after some itteration. It means that some thing is wrong; therefore, boundary condition is one of them.
Did you write an UDF for any turbulence model and get a good result? or did you see any simulation with different model by UDF and reach the same result that Fluent obtains?
I read all threads and nobody said that he gets the correct answer by UDF. may be Fluent performs some extra work on the turbulence model that is not written in the theory guide. What do you think?

Hi,
I thought so, did the simulations crashing immidiately or they run for a while?
yes , I got reasonable results from my code for turbulence!
the problem can be the marco like C_WALL_DIST(c0,t0), if the simulations crash immidiatly, this macro you can use only with enhanced wall function in standard k-epsilon model or so on... but I do not know whether it is ok for your case, you should check it...

Actually, my simulation is crash after some iteration.
I have a question about wall treatment, did you check that your code gives the same value on the wall surface?
Do you know how I hook "DEFINE_TURBULENT_VISCOSITY" in Fluent? How did you define wall function?

I have some macro like C_WALL_DIST(c0,t0), too. I check them. Thank you

hi,
to hook DEFINE_TURBULENT_VISCOSITY:
go to Models-->click on Viscous Model-->then you should see bottom right hand side User-Defined Function Turbulent Viscosity. you should hook there...
did you deactivate the turbulence from Equation panel...?
best!

Thank you for your answer. I did do, deactivate the turbulence model in equation panel. Moreover, I hook turbulent viscosity correctly.

Hello,
May I ask you about that how you consider y+ and wall shear stress near the wall in your UDF? Did you get them from Fluent by these expressions? C_UDMI(c0,t0,0)=C_STORAGE_R(f,t,SV_WALL_YPLUS_UTAU ); /* Y+*/
C_UDMI(c0,t0,2)=C_STORAGE_R(f,t,SV_WALL_SHEAR); /* wall shear */ C_UDMI(c0,t0,3)=C_STORAGE_R(f,t,SV_WALL_YPLUS); /* Ystar */

and did you save them in C_UDMI or in F_UDMI?

thanks alot for your consideration.