
[Sponsors] 
February 25, 2013, 13:40 
Define new turbulence model in Fluent

#1 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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: UDF for turbulence models And I also find some examples of kepsilon udf implementations. http://www.fhpinkafeld.ac.at/fhplus/eum/pdf/skelue.c I have some further questions by far: 1. Does anybody have examples of komega 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 Last edited by micro11sl; February 26, 2013 at 19:21. 

February 26, 2013, 00:52 

#2 
Senior Member
SSL
Join Date: Oct 2012
Posts: 227
Rep Power: 6 
It is an scalar equation. I have written so many scalars for fluent. for komega you should define at least two UDS. Then write sources and diffusion fluxes that needed for each equation.


February 26, 2013, 07:03 

#3 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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 

February 26, 2013, 07:27 

#4 
Senior Member
SSL
Join Date: Oct 2012
Posts: 227
Rep Power: 6 
If the equation you are using is general, this would also work anywhere.


February 26, 2013, 19:08 

#5 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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 komega 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" /* Userdefined constants */ #define SIG_TKE 2.0 #define SIG_OMG 2.0 /* Userdefined 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; } 

February 26, 2013, 19:11 

#6 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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_kY_k; } 

February 26, 2013, 19:12 

#7 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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_k2*C_R(c,t)*betai*fbeta*C_UDSI(c,t,OMG); return G_wY_w; } 

February 27, 2013, 14:09 

#8 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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 komega 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 komega 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 

March 2, 2013, 00:49 

#9 
Senior Member
SSL
Join Date: Oct 2012
Posts: 227
Rep Power: 6 
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. 

July 10, 2013, 05:32 

#10 
New Member
Wang
Join Date: Dec 2010
Posts: 8
Rep Power: 6 
Hi there, I am developing BSL kw 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, 10e10); } 

July 10, 2013, 06:21 

#11  
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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:


July 10, 2013, 07:45 

#12 
New Member
Wang
Join Date: Dec 2010
Posts: 8
Rep Power: 6 
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 pressurebased 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. 

July 10, 2013, 09:00 

#13  
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
Quote:


July 12, 2013, 03:08 

#14 
New Member
Wang
Join Date: Dec 2010
Posts: 8
Rep Power: 6 
I'm now trying the source linearization methods according to Menter(AIAA932906) which he think it is very robust. May be you can have a try.
for k equation: d(PkDk)/dk ~= (1/k)*Dk for omg equation: d(PwDw+Cw)/dw ~= (1/w)*(Cw+2Dw) 

January 17, 2014, 11:07 
Udf

#15 
Senior Member
rkhr
Join Date: May 2011
Posts: 226
Rep Power: 7 
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! 

January 17, 2014, 11:28 

#16 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 
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. 

January 17, 2014, 11:43 

#17  
Senior Member
rkhr
Join Date: May 2011
Posts: 226
Rep Power: 7 
Hi,
thanks for quick answer! I mean that I can not switch off default turbulence models. can you help me? thanks a lot! Quote:


January 17, 2014, 13:19 

#18 
Member
Sheng
Join Date: Jun 2011
Posts: 58
Rep Power: 6 

January 17, 2014, 13:44 

#19 
Senior Member
rkhr
Join Date: May 2011
Posts: 226
Rep Power: 7 

April 1, 2014, 14:05 
change in turbulence model, kw sst

#20 
Member
Ali.E
Join Date: Sep 2010
Location: Lisboa
Posts: 78
Rep Power: 7 
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 kequation and then, Fluent uses the modified dissipation term instead of itselfs. My turbulence model is kw 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. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
An error has occurred in cfx5solve:  volo87  CFX  5  June 14, 2013 17:44 
how to define the symmetric boundary condition for Menter's SST turbulence model?  flyingseed  Main CFD Forum  8  November 24, 2012 03:53 
What model of turbulence choose to study an external aerodynamics case  raffale  OpenFOAM  0  August 23, 2012 05:45 
Reynolds Stress Model in Fluent Vs CFX  Tim  FLUENT  0  December 6, 2005 23:03 
Sinclair Model + secondary turbulence  Yi  FLUENT  0  October 26, 2001 13:37 