[post_grad]

Hi everyone,
I am currently doing research on gas/solid fluidized beds. I have used gidaspow, and the emms model provided in the fluent package however there is a new model called the bubbling emms model which is supposedly a better algorithm for drag coefficient. I implemented the required equations from papers made a custom UDF and compiled it without errors. However, when I start the simulation I am getting a floating point error. So I started searching for parameters that might come to 0 and so on.
I think the problem I am facing is that my voidage for gas in the cells is coming to be 1 which is wrong i think. I am not sure why that's the case. I am attaching my code as a reference. Any help would be greatly appreciated.

the paper i referred to is as follows:http://dx.doi.org/10.1155/2016/8256816

/* Routine to calculate de drag coefficient using the EMMS */

#include "udf.h"
#include "sg_mphase.h"
#include "mem.h"

#define pi 4.*atan(1.)
#define diam2 57e-5
#define g 9.81
#define Dt 15e-2
#define U_mf 0.4
#define U_g 0.574
#define rho_p 1100


DEFINE_EXCHANGE_PROPERTY(bubbling_EMMS_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, z_vel_g, z_vel_s, u_b, u_e, u_ge, u_ge_1, u_ge_2, u_s, u_pe, u_g, U_slip_be, abs_vs, rho_g, rho_s, rho_e, mu_g, mu_e,
void_g, void_e, delta_b, C_dbo, C_db, a_b, Beta, Ym, Re_i, d_bm, d_b, A, B, C; /*add required variables*/

/* find the threads for the gas (primary) and solids (secondary phases) These phases appear in columns 2 and 1 in the Interphase panel respectively*/

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);
z_vel_g = C_W(cell, thread_g);

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

u_g = sqrt(x_vel_g*x_vel_g + y_vel_g*y_vel_g + z_vel_g*z_vel_g);
Message("The value of u_g is %g\n",u_g);
u_s = sqrt(x_vel_s*x_vel_s + y_vel_s*y_vel_s + z_vel_s*z_vel_s);
Message("The value of u_s is %g\n",u_s);

if (u_s == u_g)
{
u_s = 1e-4;
}

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

mu_g = C_MU_L(cell, thread_g);


/*Solution scheme starts here*/

/*Step 1 : calculation of voidage of emulsion*/
void_g = C_VOF(cell, thread_g);/* gas vol frac*/
void_e = 1-((0.58*pow((1-void_g),1.48))/(0.13+pow((1-void_g),1.48)));
Message("The value of void_g is %g\n",void_g);
Message("The value of void_e is %g\n",void_e);

/*Step 2 : calculation of volume fraction of bubbles*/
delta_b = (void_g - void_e)/(1 - void_e);

Message("The value of delta_b is %g\n",delta_b);

/*Calculating velocity of particles in emulsion (4)*/

u_pe = (u_s*(1-void_g))/((1-delta_b)*(1-void_e));


Message("The value of u_pe is %g\n",u_pe);

/*Calculating velocity of gas in emulsion eq(5)*/
A = (150*(pow((1-void_e),2)*mu_g))/(void_e*diam2*diam2);
B = (7*(1-void_e)*rho_g)/(4*diam2);
C = (1-void_e)*(rho_s-rho_g)*g*void_e;

u_ge = u_pe + ((sqrt(pow(A,2)+(4*B*C))-A)/(2*B));

Message("The value of u_ge is %g\n",u_ge);


/*

u_ge_1 = u_pe + ((sqrt(pow(A,2)+(4*B*C))-A)/(2*B));

u_ge_2 = u_pe + ((-sqrt(pow(A,2)+(4*B*C))-A)/(2*B));
if (isRealNumber(u_ge_1))
{
u_ge = u_ge_1;
}
else
{
u_ge = u_ge_2;
}


*/


/*Calculation of bubble velocity eq (6)*/

u_b = ((u_g*void_g) - ((1-delta_b)*void_e*u_ge))/(delta_b*1);
Message("The value of u_b is %g\n",u_b);


/*Calculation of dia of bubble eq (7)*/
Ym = 2.56 * (0.01*(sqrt(Dt/g))/U_mf);
Message("The value of Ym is %g\n",Ym);
d_bm = 2.59*pow(g,-0.2)*pow(((U_g-U_mf)*((pi*Dt*Dt)/4)),0.4);
Message("The value of d_bm is %g\n",d_bm);
d_b = pow(-Ym + sqrt(pow(Ym,2)+((4*d_bm)/Dt)),2)*(Dt/4);

Message("The value of d_b is %g\n",d_b);

/*Calculation of accn of bubble eq (8)*/
rho_e = rho_s*(1-void_e) + (rho_g*void_e);
mu_e = mu_g*(1 + 2.5*(1-void_e) + 10.05*pow((1-void_e),2) + (0.00273*exp(16.6*(1-void_e))));
u_e = ((rho_g*void_e*u_ge) + (rho_s*(1-void_e)*u_pe))/((rho_g*void_e)+(rho_s*(1-void_e)));



U_slip_be = (1-delta_b)*(u_b-u_e);
Re_i = (rho_e*d_b*U_slip_be)/mu_e;

if (Re_i > 1.8)
{
C_dbo = 2.7 + (24/Re_i);
}
else if (Re_i>0 && Re_i<= 1.8)
{
C_dbo = 38*pow(Re_i,-1.5);
}

C_db = pow((1-delta_b),-0.5)*C_dbo;
a_b = ((3/4)*C_db*(delta_b/d_b)*rho_e*pow(U_slip_be,2))/(delta_b*(rho_e-rho_g)) - g;

Message("The value of a_b is %g\n",a_b);

/*Calculation of Drag coeff eq (9)*/
Beta = (void_g*(delta_b*(rho_e-rho_g)*(g+a_b) + (1-delta_b)*(1-void_e)*(rho_s-rho_g)*g))/(u_g-u_s);
return (Beta);
}