[micro11sl]

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

[msaeedsadeghi]

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.

[micro11sl]

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

[msaeedsadeghi]

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

[micro11sl]

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;
}

[micro11sl]

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;
}

[micro11sl]

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;
}

[micro11sl]

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

[msaeedsadeghi]

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.

[tmac1kobe8]

Hi there, I am developing BSL k-w model by using uds recently. I met some problem that rather disturb me. As following was positive cross diffusion term to calculate the coefficient blending function F1. But the gradient term such as C_K_G and C_O_G is not exist at the beginning.
It always tell error. Would anyone give me some help.


real CD_kw_positive(cell_t c, Thread *t)
{
real a,b;
/*a = NV_DOT(C_UDSI_G(c,t,TKE), C_UDSI_G(c,t,OMG)); */
a = NV_DOT(C_K_G(c,t), C_O_G(c,t) /*calculate (dk/dx_j)*(dw/dx_j)*/
b = 2.0*C_R(c,t) / (SIG_OMG_2*C_UDSI(c,t,OMG)*a);

return MAX(b, 10e-10);
}

[micro11sl]

You have to let Fluent store those values in memory during calculation. This can be done by typing these in the command line:

/solve/set/expert

and then answering "yes" for the question "Keep temporary solver memory from being freed?"

Try this to see if it helps.

Regards,
Sheng

Originally Posted by tmac1kobe8

Hi there, I am developing BSL k-w model by using uds recently. I met some problem that rather disturb me. As following was positive cross diffusion term to calculate the coefficient blending function F1. But the gradient term such as C_K_G and C_O_G is not exist at the beginning.
It always tell error. Would anyone give me some help.


real CD_kw_positive(cell_t c, Thread *t)
{
real a,b;
/*a = NV_DOT(C_UDSI_G(c,t,TKE), C_UDSI_G(c,t,OMG)); */
a = NV_DOT(C_K_G(c,t), C_O_G(c,t) /*calculate (dk/dx_j)*(dw/dx_j)*/
b = 2.0*C_R(c,t) / (SIG_OMG_2*C_UDSI(c,t,OMG)*a);

return MAX(b, 10e-10);
}

[tmac1kobe8]

Yes, I've found that. So I have to iterate some steps to make Fluent remember the gradient values, then I put the udfs in.

But the next problem I met is that I have to use pressure-based solver since I've used macros such as F_CENTROID, otherwise it will ended with error.
The initial procedure is always turn me down, the max number of iteration steps is 7. Do you have some good suggestions?

Best wishes.
WANG.

[micro11sl]

Originally Posted by tmac1kobe8

Yes, I've found that. So I have to iterate some steps to make Fluent remember the gradient values, then I put the udfs in.

But the next problem I met is that I have to use pressure-based solver since I've used macros such as F_CENTROID, otherwise it will ended with error.
The initial procedure is always turn me down, the max number of iteration steps is 7. Do you have some good suggestions?

Best wishes.
WANG.

I have startup problem too. And I've yet managed to sort it. Source term linearisation may help. But I am not sure how to do it effectively.

[tmac1kobe8]

I'm now trying the source linearization methods according to Menter(AIAA-93-2906) which he think it is very robust. May be you can have a try.

for k equation: d(Pk-Dk)/dk ~= -(1/k)*Dk
for omg equation: d(Pw-Dw+Cw)/dw ~= -(1/w)*(|Cw|+2Dw)

[Kanarya]

Hi All,

I have written udf code for turbulence model but I couldnt deactivate the existing ones in fluent and it gives me access violation error.
In order to use udf turb. model should I deactivate the turbulence models?
is there anyone who can help me to add my code to the fluent properly?

thanks in advance!

[micro11sl]

Hi Kanaya,
What do you mean that "can't deactive existing turbulence models"? When you use your own turbulence models, you have to possibly define necessary UDM when necessary. This is one reason what I know about "access violation". You also need to think about the procedure when solving it, if some quantities are not available while the equation needs them, such errors will happen. They way to solve it is to compute those quantities first in your code. Other stuff include boundary condition, switching off solving defalt turbulence model by uncheck the equation for solving.

I didn't continue my effort on implementing a turbulence model into Fluent for about a year. I later have found some errors in the code I posted here, so DO NOT USE THIS ONE. Divergence always happen when I run the simulation with my UDF turbulence models. I need more time in the future to get it sorted. One difficulty I have found is that these turbulene equations are coupled. In k equations, omega needs to be known and vice versa. However, we can't control which equation to be solved first.

[Kanarya]

Hi,

thanks for quick answer!
I mean that I can not switch off default turbulence models.
can you help me?
thanks a lot!

Originally Posted by micro11sl

Hi Kanaya,
What do you mean that "can't deactive existing turbulence models"? When you use your own turbulence models, you have to possibly define necessary UDM when necessary. This is one reason what I know about "access violation". You also need to think about the procedure when solving it, if some quantities are not available while the equation needs them, such errors will happen. They way to solve it is to compute those quantities first in your code. Other stuff include boundary condition, switching off solving defalt turbulence model by uncheck the equation for solving.

I didn't continue my effort on implementing a turbulence model into Fluent for about a year. I later have found some errors in the code I posted here, so DO NOT USE THIS ONE. Divergence always happen when I run the simulation with my UDF turbulence models. I need more time in the future to get it sorted. One difficulty I have found is that these turbulene equations are coupled. In k equations, omega needs to be known and vice versa. However, we can't control which equation to be solved first.

[micro11sl]

Originally Posted by Kanarya

Hi,

thanks for quick answer!
I mean that I can not switch off default turbulence models.
can you help me?
thanks a lot!

solution control -> equations -> uncheck turbulence

[Kanarya]

thanks a lot! still Have the same problem
if you have email address I can send you the code and case!
if you have time
thanks in advance!

Originally Posted by micro11sl

solution control -> equations -> uncheck turbulence

[behest]

Hello all,
I want to modify the dissipation rate of kinetic energy equation by UDF. Which macro can be used? Is it possible to use DEFINE_SOURCE for definning new dissipation rate?

Actually, I must change the formulation of dissipation rate of k-equation and then, Fluent uses the modified dissipation term instead of itselfs. My turbulence model is k-w SST and I just manipulate Y_k (dissipation rate of k).
Original: Y_k=ro*beta_s*k*w;
Modified: Y_k=ro*beta_s*k*w/delta; (delta is wall distance)

It would be appriciated if anyone makes a comment to how I do this task.