[ravindra]

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]

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

Daniele

[ravindra]

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]

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]

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

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

[ravindra]

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]

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]

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]

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

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

Thanks again.

[Musa]

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

Ok.Thanks.

[mohsen0488]

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]

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]

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

[mohsen0488]

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]

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]

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]

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]

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