CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (http://www.cfd-online.com/Forums/fluent-udf/)
-   -   Gidaspow udf (http://www.cfd-online.com/Forums/fluent-udf/116318-gidaspow-udf.html)

ravindra April 17, 2013 01:26

Gidaspow udf
 
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.:)

ghost82 April 17, 2013 03:02

Maybe it could be of great help if you attach your udf.

Daniele

ravindra April 17, 2013 11:18

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

ghost82 April 18, 2013 04:31

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 April 18, 2013 04:41

In addition I think it is better to add a fabs:

Code:

abs_v = sqrt(fabs(slip_x*slip_x + slip_y*slip_y));

ravindra April 19, 2013 01:26

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 April 25, 2013 05:07

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

mohsen0488 October 14, 2013 03:27

Hi
i want to write a udf about "gidaspow drag model with switch function" can you help me please


thanks a lot in advance.

Musa October 15, 2013 06:19

@ 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.

mohsen0488 October 15, 2013 06:38

Thanks a lot Mr Musa for your reply.
I will try and if I had a problem I ask you.

Thanks again.

Musa October 15, 2013 06:53

@ 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.

mohsen0488 October 15, 2013 07:14

Ok.Thanks.

mohsen0488 October 17, 2013 11:48

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

Musa October 17, 2013 19:21

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

Musa October 17, 2013 19:23

Try my code and let me know if you get the desired results. All the best.

mohsen0488 October 18, 2013 06:47

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?

Musa October 18, 2013 07:52

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!

mohsen0488 October 18, 2013 10:52

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

Musa October 18, 2013 17:46

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.

mohsen0488 October 18, 2013 18:20

Musa I tried with that code and the results were similar results of paper
so do you think results of this paper are incorrect?


All times are GMT -4. The time now is 08:06.