CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (https://www.cfd-online.com/Forums/fluent-udf/)
-   -   Define new turbulence model in Fluent (https://www.cfd-online.com/Forums/fluent-udf/113720-define-new-turbulence-model-fluent.html)

micro11sl February 25, 2013 12:40

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

msaeedsadeghi February 25, 2013 23:52

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 February 26, 2013 06:03

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 February 26, 2013 06:27

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

micro11sl February 26, 2013 18:08

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 February 26, 2013 18:11

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 February 26, 2013 18:12

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 February 27, 2013 13:09

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 March 1, 2013 23:49

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 July 10, 2013 05:32

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 July 10, 2013 06:21

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



Quote:

Originally Posted by tmac1kobe8 (Post 438880)
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 July 10, 2013 07:45

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 July 10, 2013 09:00

Quote:

Originally Posted by tmac1kobe8 (Post 438921)
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 July 12, 2013 03:08

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 January 17, 2014 10:07

Udf
 
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 January 17, 2014 10:28

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 January 17, 2014 10:43

Hi,

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

Quote:

Originally Posted by micro11sl (Post 470523)
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 January 17, 2014 12:19

Quote:

Originally Posted by Kanarya (Post 470526)
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 January 17, 2014 12:44

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!
Quote:

Originally Posted by micro11sl (Post 470541)
solution control -> equations -> uncheck turbulence


behest April 1, 2014 14:05

change in turbulence model, k-w sst
 
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.

micro11sl April 1, 2014 15:27

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!

behest April 2, 2014 06:53

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.


Quote:

Originally Posted by micro11sl (Post 483296)
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!


Kanarya August 29, 2014 17:20

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!
Quote:

Originally Posted by micro11sl (Post 410306)
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 August 30, 2014 06:34

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

Regards,

Quote:

Originally Posted by Kanarya (Post 508407)
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!


Kanarya August 30, 2014 06:46

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!
Quote:

Originally Posted by micro11sl (Post 508434)
Hi Kanarya,
Sorry I didn't understand. Can you illustrate better?

Regards,


micro11sl September 1, 2014 04:47

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,

Quote:

Originally Posted by Kanarya (Post 508435)
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!


Kanarya September 1, 2014 06:40

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

Originally Posted by micro11sl (Post 508593)
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,


micro11sl September 1, 2014 12:28

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,

Quote:

Originally Posted by Kanarya (Post 508618)
hi,
thanks!
I have convergence problem as well....do you thing i should initialize the tensors?
thanks!


behest November 12, 2014 09:28

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.

Kanarya November 12, 2014 10:04

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

Originally Posted by behest (Post 518754)
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.


behest November 12, 2014 11:06

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.

Quote:

Originally Posted by Kanarya (Post 518767)
Hi ,
If you use the fluent theory guide, you should get similar results...


Kanarya November 12, 2014 11:34

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!
Quote:

Originally Posted by behest (Post 518784)
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.


behest November 13, 2014 05:46

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.
Quote:

Originally Posted by Kanarya (Post 518790)
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!


Kanarya November 13, 2014 06:27

Hi,

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

Originally Posted by behest (Post 518952)
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.


behest November 13, 2014 09:41

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?

Kanarya November 13, 2014 10:00

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...
Quote:

Originally Posted by behest (Post 519000)
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?


behest November 13, 2014 18:18

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


Quote:

Originally Posted by Kanarya (Post 519008)
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...


Kanarya November 14, 2014 12:52

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!
Quote:

Originally Posted by behest (Post 519085)
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


behest November 14, 2014 13:01

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

Quote:

Originally Posted by Kanarya (Post 519205)
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!


behest November 20, 2014 07:18

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.

Quote:

Originally Posted by Kanarya (Post 519205)
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!



All times are GMT -4. The time now is 14:11.