Define new turbulence model in Fluent
 

Note:
26/2/13: The essential UDF codes other than the wall boundary condition are uploaded from post #5 to #7. Any suggestions or comments are welcome.

Original post in 25/2/13:
Hi all,
I'm going to define a new turbulence model in Fluent instead of using the default ones.
I have done some search and found a related post is:
http://www.cfd-online.com/Forums/flu...ce-models.html
And I also find some examples of k-epsilon udf implementations.
http://www.fh-pinkafeld.ac.at/fhplus/eum/pdf/ske-lue.c

I have some further questions by far:
1. Does anybody have examples of k-omega udf implementations? I will eventually implement the correlation based γ-Reθ SST model.
2. As a general strategy for defining new turbulence model in Fluent, is it always necessary to select the default model first , add user defined UDS, and disable the solution for the default transport equations for turbulence? Any alternatives?

Any replies are welcomed. I plan to keep updating the status of defining the turbulence model with udfs.

Regards,
Sheng

It is an scalar equation. I have written so many scalars for fluent. for k-omega you should define at least two UDS. Then write sources and diffusion fluxes that needed for each equation.

A quick question:
I am going to consider the equations are still valid for compressible flows, so I need to retain density in the equations. Can I use this udf valid for compressible flows for incompressible flows? Do I need to write udfs for incompressible and compressible flows separately?

Cheers,
Sheng

If the equation you are using is general, this would also work anywhere.

I have finished essential parts except the wall boundary condition today. I am uploading in separated posts. I haven't debug it yet. I'll do it in following days and keep updating.

The k-omega model formulation follows the one presented in the Fluent theory guide.
The first part is as follows. It deals with eddy viscosity and diffusivity.

#include "udf.h" // This header is required.
#include "mem.h" // To access density, etc.
#include "math.h"
/* User-defined constants */
#define SIG_TKE 2.0
#define SIG_OMG 2.0
/* User-defined scalars */
enum
{
TKE,
OMG,
N_REQUIRED_UDS
};

DEFINE_ADJUST(eddy_viscosity, domain) /* Consider using DEFINE_TURBULENT_VISCOSITY?*/
{
Thread *t;
cell_t c;
real alpha_star=1.0;
/* Set the turbulent viscosity, looping over all threads and then cells */
thread_loop_c (t, domain)
{
if (FLUID_THREAD_P(t))
{
begin_c_loop(c,t)
{
C_MU_T(c,t) =alpha_star*C_R(c,t)*C_UDSI(c,t,TKE)/C_UDSI(c,t,OMG);
}
end_c_loop(c,t)
}
}
}

DEFINE_DIFFUSIVITY(kw_diff, c, t, eqn)
{
real diff; // define the diffusion coeffcient
real mu = C_MU_L(c,t);
real mu_t = C_MU_T(c,t);
switch (eqn)
{
case TKE:
diff = mu_t/SIG_TKE + mu;
break;
case OMG:
diff = mu_t/SIG_OMG + mu;
break;
default:
diff = mu_t + mu;
}
return diff;
}

The second part deals with the source term of k: production - dissipation

DEFINE_SOURCE(k_source, c, t, dS, eqn)
{
int i,j,m;
real G_k; // refers to the production of k
real Y_k; // refers to the dissipation of k
real beta_star=0.09; // the compressibility correction is not enabled for k
real betai=0.072; // the compressibility correction is not enabled for w
real Xk, gk, gw;

G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t)); // Bounsinesq hypothesis, production of k
gdk = C_UDSI_G(c,t,TKE); // retrieve the gradient of k
gdw = C_UDSI_G(c,t,OMG); // retrieve the gradient of w
Xk = pow(C_UDSI(c,t,OMG),-3.0)*(gdk[0]*gdw[0]+gdk[1]*gdw[1]+gdk[2]*gdw[2]);
if(Xk>0)
{
fbeta_star=(1+680*pow(Xk,2))/(1+400*pow(Xk,2));
}
else
{
fbeta_star=0;
}
Y_k = C_R(c,t)*beta_star*fbeta_star*C_UDSI(c,t,TKE)*C_UD SI(c,t,OMG);
dS[eqn] = -C_R(c,t)*beta_star*fbeta_star*C_UDSI(c,t,OMG);
return G_k-Y_k;
}

The third part deals with the source term of omega, I use the denotation of w in somewhere.

DEFINE_SOURCE(w_source, c, t, dS, eqn)
{
int i,j,m;
real G_k; // refers to the production of k
real G_w; // refers to the production of w
real Y_w; // refers to the dissipation of w
real alpha=1.0;
real beta_star=0.09; // the compressibility correction is not enabled for k
real betai=0.072; // the compressibility correction is not enabled for w

real Xw=0;
real Og[ND_ND][ND_ND], S[ND_ND][ND_ND], puixj[ND_ND][ND_ND];

puixj[0][0]= C_DUDX(c,t);
puixj[0][1]= C_DVDX(c,t);
puixj[1][0]= C_DUDY(c,t);
puixj[1][1]= C_DVDY(c,t);
if(ND_ND==3)
{
puixj[0][2]= C_DWDX(c,t);
puixj[1][2]= C_DWDY(c,t);
puixj[2][0]= C_DUDZ(c,t);
puixj[2][1]= C_DVDZ(c,t);
puixj[2][2]= C_DWDZ(c,t);
}
G_k = C_MU_T(c,t)*SQR(Strainrate_Mag(c,t)); // Bounsinesq hypothesis
G_w = alpha*C_UDSI(c,t,OMG)/C_UDSI(c,t,TKE)*G_k; // production of w, requiring G_k
for(m=0;m<ND_ND;m++)
{
for(j=0:j<ND_ND;j++)
{
for(i=0;i<ND_ND;i++)
{
Og[i][j]=0.5*(puixj[i][j]-puixj[j][i]); //Og is the rotation tensor
Og[j][m]=0.5*(puixj[j][m]-puixj[m][j]);
S[m][i]=0.5*(puixj[m][i]+puixj[i][m]); //S is the stress tensor
Xw = Og[i][j]*Og[j][m]*S[m][i]+ Xw;
}
}
}
Xw=abs(Xw/pow(beta_star*C_UDSI(c,t,OMG),3));
fbeta=(1+70*Xw)/(1+80*Xw);
Y_w = C_R(c,t)*betai*fbeta*pow(C_UDSI(c,t,OMG),2);
dS[eqn] = alpha/C_UDSI(c,t,TKE)*G_k-2*C_R(c,t)*betai*fbeta*C_UDSI(c,t,OMG);
return G_w-Y_w;
}

Hi all,
I get stuck today. I find there's more work that I expected before.

Question 1:
Because I am implementing a compressible turbulence model, does this imply that I need to define an energy equation as well? Because turbulent kinetic energy should be presented in the energy equation.

Question 2:
After I load the UDF for user defined scalars, when setting the boundary condition tab, is there any difference between the "specific value" and "specific flux" if I am going to use my own boundary condition profile?

Question 3:
How to get the "turbulence intensity" and "viscosity ratio" I set for k-omega based model before for my newly defined model? Do I need to calculate the value of k and omega explicitly and assign them to my UDS? A related question is do I need to write up an udf for the initialization of my UDS? Also udfs for postprocessing?

Question 4:
As I read from many other threads in this forum, I guess some modification should be done before the udf can be run in parallel. Is this guess correct?

Considering from Question 1 to 4, I feel there's lots of work to finish. It seems I can't finish shortly. To simplify, I have a very rough idea. Because the original k-omega model will be selected but not solved (the UDS equations are solved instead), I might transfer the value of the initialized k and omega to my UDS in the beginning, and transfer my UDS back to the original k and omega scalars at the end of one computation. In this way, there's no need to write up udfs to initialize and postprocessing my UDS. But I don't know it's possible or not that Fluent will let me assign values back and forth between UDS and original k and omega transport equations.

Are there any simple approaches? Any comments?

Regards,
Sheng

I have written so User Defined Scalars before.
There is difference between flux and value. It's clear.
There were no need to change UDF for parallel in my previous UDFs.