CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Fluent UDF and Scheme Programming

Gidaspow udf

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree4Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   April 17, 2013, 01:26
Default Gidaspow udf
  #1
New Member
 
Join Date: Jul 2012
Posts: 19
Rep Power: 5
ravindra is on a distinguished road
hey hi friends,
i have a little problem but i am unable to figure it out
i am using this gidaspow drag model for my fluidised bed simulation but i am getting different results when i am incorporating this drag law through inbuilt gidaspow law and by writing a UDF for it.
i have seen from fluent manual that the equations which i have used are same in the udf as given in manual.

Can anybody please tell me what may be the reason behind this?

Thank you in advance.
ravindra is offline   Reply With Quote

Old   April 17, 2013, 03:02
Default
  #2
Senior Member
 
ghost82's Avatar
 
Daniele
Join Date: Oct 2010
Location: Italy
Posts: 918
Rep Power: 15
ghost82 will become famous soon enough
Maybe it could be of great help if you attach your udf.

Daniele
ghost82 is offline   Reply With Quote

Old   April 17, 2013, 11:18
Default
  #3
New Member
 
Join Date: Jul 2012
Posts: 19
Rep Power: 5
ravindra is on a distinguished road
this is my udf for gidaspow model
#include "udf.h"

#define diam2 875.e-7


DEFINE_EXCHANGE_PROPERTY(custom_drag,cell,mix_thre ad,s_col,f_col)
{
Thread *thread_g, *thread_s;
real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y,rho_g, rho_s, mu_g, reyp;
real void_g,void_s,fdrgs,k_g_s;

/* find the threads for the gas (primary) */
/* and solids (secondary phases) */

thread_g = THREAD_SUB_THREAD(mix_thread, s_col);/* gas phase */
thread_s = THREAD_SUB_THREAD(mix_thread, f_col);/* solid phase*/

/* find phase velocities and properties*/

x_vel_g = C_U(cell, thread_g);
y_vel_g = C_V(cell, thread_g);

x_vel_s = C_U(cell, thread_s);
y_vel_s = C_V(cell, thread_s);

slip_x = x_vel_g - x_vel_s;
slip_y = y_vel_g - y_vel_s;

rho_g = C_R(cell, thread_g);
rho_s = C_R(cell, thread_s);

mu_g = C_MU_L(cell, thread_g);

/*compute slip*/
abs_v = sqrt(slip_x*slip_x + slip_y*slip_y);

/*compute Reynold's number*/
reyp = rho_g*abs_v*diam2/mu_g;

void_g = C_VOF(cell, thread_g);/* gas vol frac*/
void_s=1.-void_g;

if(void_g>0.8)
{
if(reyp<=1000)
{
fdrgs=(24/(void_g*reyp))*(1+0.15*(pow((void_g*reyp),0.687))) ;
}
else
{
fdrgs=0.44;
}

k_g_s=(3/4)*(fdrgs*void_s*void_g*rho_g*abs_v)*pow(void_g,-2.65)/diam2;
}
else
{
k_g_s=150*void_s*void_s*mu_g/(void_g*diam2*diam2)+1.75*void_s*rho_g*abs_v/diam2;
}
return k_g_s;
}
ravindra is offline   Reply With Quote

Old   April 18, 2013, 04:31
Default
  #4
Senior Member
 
ghost82's Avatar
 
Daniele
Join Date: Oct 2010
Location: Italy
Posts: 918
Rep Power: 15
ghost82 will become famous soon enough
Hi,
I read the fluent user guide; when you compute fdrgs you are discriminating For the particle Reynolds number, but as you can see here (https://www.sharcnet.ca/Software/Flu...mp-cd-gidaspow) the Gidaspow model calculates fdrgs with one formula:

fdrgs=(24/(void_g*reyp))*(1+0.15*(pow((void_g*reyp),0.687)))

Accordingly to the user guide your calculations of Cd refer to the "symmetric" or to the "Shiller and Naumann" models.
Maybe this is the "error" in your udf.

Daniele
ghost82 is offline   Reply With Quote

Old   April 18, 2013, 04:41
Default
  #5
Senior Member
 
ghost82's Avatar
 
Daniele
Join Date: Oct 2010
Location: Italy
Posts: 918
Rep Power: 15
ghost82 will become famous soon enough
In addition I think it is better to add a fabs:

Code:
abs_v = sqrt(fabs(slip_x*slip_x + slip_y*slip_y));
ghost82 is offline   Reply With Quote

Old   April 19, 2013, 01:26
Default
  #6
New Member
 
Join Date: Jul 2012
Posts: 19
Rep Power: 5
ravindra is on a distinguished road
hey Daniele,
i have also tried in which i didnot discriminate for reynolds number...but still the results are different for cases with udf and without udf.

i will do what you have suggested for using fabs() function and get back to you...

thanks
ravindra is offline   Reply With Quote

Old   April 25, 2013, 05:07
Default
  #7
New Member
 
Join Date: Jul 2012
Posts: 19
Rep Power: 5
ravindra is on a distinguished road
hey daniele,
i am using this udf but results are still different with those obtained without udf.
can you suggest what is wrong.

#include "udf.h"

#define diam2 875.e-7


DEFINE_EXCHANGE_PROPERTY(custom_drag,cell,mix_thre ad,s_col,f_col)
{
Thread *thread_g, *thread_s;
real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y,rho_g, rho_s, mu_g, reyp;
real void_g,void_s,fdrgs,k_g_s;

/* find the threads for the gas (primary) */
/* and solids (secondary phases) */

thread_g = THREAD_SUB_THREAD(mix_thread, s_col);/* gas phase */
thread_s = THREAD_SUB_THREAD(mix_thread, f_col);/* solid phase*/

/* find phase velocities and properties*/

x_vel_g = C_U(cell, thread_g);
y_vel_g = C_V(cell, thread_g);

x_vel_s = C_U(cell, thread_s);
y_vel_s = C_V(cell, thread_s);

slip_x = x_vel_g - x_vel_s;
slip_y = y_vel_g - y_vel_s;

rho_g = C_R(cell, thread_g);
rho_s = C_R(cell, thread_s);

mu_g = C_MU_L(cell, thread_g);

/*compute slip*/
abs_v = sqrt(fabs(slip_x*slip_x + slip_y*slip_y));

/*compute Reynold's number*/
reyp = rho_g*abs_v*diam2/mu_g;

void_g = C_VOF(cell, thread_g);/* gas vol frac*/
void_s=1.-void_g;

if(void_g>0.8)
{

fdrgs=(24/(void_g*reyp))*(1+0.15*(pow((void_g*reyp),0.687))) ;

k_g_s=(3/4)*(fdrgs*void_s*void_g*rho_g*abs_v)*pow(void_g,-2.65)/diam2;
}
else
{
k_g_s=150*void_s*void_s*mu_g/(void_g*diam2*diam2)+1.75*void_s*rho_g*abs_v/diam2;
}
return k_g_s;
}
ravindra is offline   Reply With Quote

Old   October 14, 2013, 03:27
Default
  #8
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
Hi
i want to write a udf about "gidaspow drag model with switch function" can you help me please


thanks a lot in advance.
mohsen0488 is offline   Reply With Quote

Old   October 15, 2013, 06:19
Default
  #9
Member
 
Musango Lungu
Join Date: Dec 2011
Location: China
Posts: 73
Rep Power: 5
Musa is on a distinguished road
@ Mohsen it's very easy to do that. Bear in mind that there are several stitching functions have been proposed for the Gidaspow model with the most commonly used being that due to Hulin and Gidaspow. The model is actually available in ANSYS FLUENT 14. Alternatively you can browse the source codes for the drag models in MFIX code and get an idea how to code it for FLUENT. If you still have trouble or doubts then feel free to contact me and I will you the specific details.
Musa is offline   Reply With Quote

Old   October 15, 2013, 06:38
Default
  #10
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
Thanks a lot Mr Musa for your reply.
I will try and if I had a problem I ask you.

Thanks again.
mohsen0488 is offline   Reply With Quote

Old   October 15, 2013, 06:53
Default
  #11
Member
 
Musango Lungu
Join Date: Dec 2011
Location: China
Posts: 73
Rep Power: 5
Musa is on a distinguished road
@ Mohsen you are welcome. Look at this paper Leboreiro, J., Joseph, G. G., & Hrenya, C. M. Revisiting the standard drag law for bubbling, gas-fluidized beds. Powder Technol 2008, 183, 385-400. for the discussion regarding the different stitching functions to address the discontinuity of the Gidaspow model. For browsing the MFIX code you have to be a member to have access. It's easy and free to join. Membership will usually take a day or two to be activated uopn application. Though the source code is in FORTRAN it can be easily inferred for use in C.
Musa is offline   Reply With Quote

Old   October 15, 2013, 07:14
Default
  #12
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
Ok.Thanks.
mohsen0488 is offline   Reply With Quote

Old   October 17, 2013, 11:48
Default
  #13
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
Dear Musa
I read article that said you, and I wrote the code as follows:
But my results are different from other authors , and i don't know what's the problem can you help me please? Thanks

#include "udf.h"
#define pi 4.*atan(1.)
#define diam2 1.41e-3
DEFINE_EXCHANGE_PROPERTY(custom_drag, cell, mix_thread, s_col, f_col)
{
Thread *thread_g, *thread_s;
real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y,
rho_g, rho_s, mu_g, reyp, cd,
void_g,void_s, vfac, fi_gs, k_g_s, k_g_s_eg, k_g_s_wy;
/* find the threads for the gas (primary) */
/* and solids (secondary phases) */
thread_g = THREAD_SUB_THREAD(mix_thread, s_col);/* gas phase */
thread_s = THREAD_SUB_THREAD(mix_thread, f_col);/* solid phase*/
/* find phase velocities and properties*/
x_vel_g = C_U(cell, thread_g);
y_vel_g = C_V(cell, thread_g);
x_vel_s = C_U(cell, thread_s);
y_vel_s = C_V(cell, thread_s);
slip_x = x_vel_g - x_vel_s;
slip_y = y_vel_g - y_vel_s;
rho_g = C_R(cell, thread_g);
rho_s = C_R(cell, thread_s);
mu_g = C_MU_L(cell, thread_g);
/*compute slip*/
abs_v = sqrt(slip_x*slip_x + slip_y*slip_y);
void_g = C_VOF(cell, thread_g);/* gas vol frac*/
void_s = C_VOF(cell, thread_s);/* solid vol frac*/
/*compute reynolds number*/
reyp = rho_g*abs_v*diam2/mu_g;
/*compute Cd*/
if(reyp<=1000)
cd=24*(1+0.15*pow(void_g*reyp,0.687))/void_g*reyp;
else
cd=0.44;
fi_gs=atan(150*1.75*(void_g-0.8))/pi+0.5;
/*compute drag and return drag coeff, k_g_s*/
k_g_s_eg=(150.*(pow(void_s, 2)*mu_g)/(void_g*pow(diam2, 2)))+(1.75*void_s*rho_g*abs_v)/diam2;
k_g_s_wy=((3./4.)*cd*void_s*void_g*rho_g*abs_v*pow(void_g, -2.65))/diam2;
k_g_s=fi_gs*k_g_s_wy+(1-fi_gs)*k_g_s_eg;
return k_g_s;
}
mohsen0488 is offline   Reply With Quote

Old   October 17, 2013, 19:21
Default
  #14
Member
 
Musango Lungu
Join Date: Dec 2011
Location: China
Posts: 73
Rep Power: 5
Musa is on a distinguished road
Hi Mohsen you forgot the void fraction condition that switches from Ergun to Wen and Yu. Here is my code below. Also notice how I define the particle diameter inside the code to avoid changing the particle diameter every time I have a new particle size.

#include "udf.h"
#define pi 4.*atan(1.)


DEFINE_EXCHANGE_PROPERTY(gidaspow_blend,cell,mix_t hread,s_col,f_col)
{
Thread *thread_g, *thread_s;
real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y,
rho_g, rho_s, mu_g, reyp, Res, afac, Cd,
vof_g, vof_s, k_g_s,Ergun,WenYu,phi_gs, diam2;

/* find the threads for the gas (primary) */
/* and solids (secondary phases) */

thread_g = THREAD_SUB_THREAD(mix_thread, s_col);/* gas phase */
thread_s = THREAD_SUB_THREAD(mix_thread, f_col);/* solid phase*/

/* find phase velocities and properties*/

x_vel_g = C_U(cell, thread_g);
y_vel_g = C_V(cell, thread_g);

x_vel_s = C_U(cell, thread_s);
y_vel_s = C_V(cell, thread_s);

slip_x = x_vel_g - x_vel_s;
slip_y = y_vel_g - y_vel_s;

rho_g = C_R(cell, thread_g);
rho_s = C_R(cell, thread_s);

mu_g = C_MU_L(cell, thread_g);

diam2 = C_PHASE_DIAMETER(cell, thread_s);

/*compute slip*/
abs_v = sqrt(slip_x*slip_x + slip_y*slip_y);

/*compute Reynold's number*/

reyp = rho_g*abs_v*diam2/mu_g;

/* compute particle relaxation time */

vof_g = C_VOF(cell, thread_g);/* gas vol frac*/
vof_s = C_VOF(cell, thread_s);/* particle vol frac*/

Res = vof_g * reyp;

if(Res>=1000)
Cd=0.44;
else
Cd=24/Res*(1+0.15*pow(Res,0.687));


if(vof_g <=0.8){
/* Ergun Drag Model */
Ergun = 150*vof_s*(1-vof_g)/vof_g*mu_g/diam2/diam2 + 1.75*vof_s*rho_g/diam2*abs_v;
}


else
{
/* Wen & Yu Model */
WenYu = 3/4*Cd*vof_g*vof_s*rho_g*abs_v/diam2*pow(vof_g,-2.65);
}

phi_gs = atan(150*1.75*(0.2-vof_s))/pi +0.5;

k_g_s = (1-phi_gs)*Ergun + phi_gs*WenYu;

return k_g_s;
}
coolgame2004 and esmaeil like this.
Musa is offline   Reply With Quote

Old   October 17, 2013, 19:23
Default
  #15
Member
 
Musango Lungu
Join Date: Dec 2011
Location: China
Posts: 73
Rep Power: 5
Musa is on a distinguished road
Try my code and let me know if you get the desired results. All the best.
esmaeil likes this.
Musa is offline   Reply With Quote

Old   October 18, 2013, 06:47
Default
  #16
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
Dear Musa
Thank you so much
I tried with your code , but still the results are different from other authors
Do you think the code has still a problem?
esmaeil likes this.
mohsen0488 is offline   Reply With Quote

Old   October 18, 2013, 07:52
Default
  #17
Member
 
Musango Lungu
Join Date: Dec 2011
Location: China
Posts: 73
Rep Power: 5
Musa is on a distinguished road
Hi Mohsen the code does not have a problem. I have coded it as given in several references. I also made reference to the MFIX source code of this drag model. Results might be different due to other reasons perhaps. You can go through the authors complete CFD code and see if matches with yours. Ensure that the boundary conditions and other parameters are the same. Check your Grid also. If you still don't get the same results then perhaps you can contact the authors of the work you are trying to mimic and ask them for the specifics. Good luck!
Musa is offline   Reply With Quote

Old   October 18, 2013, 10:52
Default
  #18
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
hi musa
this paper " CFD modeling of the gas–particle flow behavior in spouted beds , Wu Zhonghua ⁎, Arun S. Mujumdar, " used this code as below: i tried by this code and my results are same.but i think it has a problem in solid reynolds definition do you think too?
code as follows:
#include "udf.h"
#define pi 4.*atan(1.)
#define diam2 1.41e-3
DEFINE_EXCHANGE_PROPERTY(custom_drag, cell, mix_thread, s_col, f_col)
{
Thread *thread_g, *thread_s;
real x_vel_g, x_vel_s, y_vel_g, y_vel_s, abs_v, slip_x, slip_y,
rho_g, rho_s, mu_g, reyp, cd,
void_g, vfac, fi_gs, k_g_s, k_g_s_eg, k_g_s_wy;
/* find the threads for the gas (primary) */
/* and solids (secondary phases) */
thread_g = THREAD_SUB_THREAD(mix_thread, s_col);/* gas phase */
thread_s = THREAD_SUB_THREAD(mix_thread, f_col);/* solid phase*/
/* find phase velocities and properties*/
x_vel_g = C_U(cell, thread_g);
y_vel_g = C_V(cell, thread_g);
x_vel_s = C_U(cell, thread_s);
y_vel_s = C_V(cell, thread_s);
slip_x = x_vel_g - x_vel_s;
slip_y = y_vel_g - y_vel_s;
rho_g = C_R(cell, thread_g);
rho_s = C_R(cell, thread_s);
mu_g = C_MU_L(cell, thread_g);
/*compute slip*/
abs_v = sqrt(slip_x*slip_x + slip_y*slip_y);
void_g = C_VOF(cell, thread_g);/* gas vol frac*/
/*compute reynolds number*/
reyp = rho_g*abs_v*diam2/mu_g;
/*compute Cd*/
if(reyp<=1000)
cd=24*(1+0.15*pow(reyp,0.687))/reyp;
else
cd=0.44;
fi_gs=atan(150*1.75*(void_g-0.8))/pi+0.5;
/*compute drag and return drag coeff, k_g_s*/
k_g_s_eg=(150*(1-void_g)+1.75*reyp)*(1-void_g)*mu_g/void_g/diam2/diam2;
k_g_s_wy=0.75*cd*(1-void_g)*mu_g*reyp*pow(void_g,-2.65)/diam2/diam2;
k_g_s=fi_gs*k_g_s_wy+(1-fi_gs)*k_g_s_eg;
return k_g_s;
}
mohsen0488 is offline   Reply With Quote

Old   October 18, 2013, 17:46
Default
  #19
Member
 
Musango Lungu
Join Date: Dec 2011
Location: China
Posts: 73
Rep Power: 5
Musa is on a distinguished road
Mohsen you are absolutely right the reyp has to be multiplied by void_g in the equation containing Cd but this has not been done in the code although it has been shown in the paper . Secondly there is no condition for switching from one regime to the other i.e from Ergun for the dense bed (void_g > 0.8) to WenYu for ( void_g<= 0.8). The stitching function simply ensures smooth transition from one regime to the other meaning that the condition of regime change must still be given.
Musa is offline   Reply With Quote

Old   October 18, 2013, 18:20
Default
  #20
Member
 
mohsen
Join Date: Sep 2013
Posts: 34
Rep Power: 3
mohsen0488 is on a distinguished road
Musa I tried with that code and the results were similar results of paper
so do you think results of this paper are incorrect?
mohsen0488 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
UDF parallel error: chip-exec: function not found????? shankara.2 Fluent UDF and Scheme Programming 1 January 16, 2012 23:14
How to add a UDF to a compiled UDF library kim FLUENT 3 October 26, 2011 21:38
Dynamic Mesh UDF Qureshi FLUENT 1 December 2, 2009 01:27
UDF...UDF...UDF...UDF Luc SEMINEL FLUENT 0 November 25, 2002 05:03
UDF, UDF, UDF, UDF Luc SEMINEL Main CFD Forum 0 November 25, 2002 05:01


All times are GMT -4. The time now is 15:35.