[Luca_Panariello]

Hi everybody
I am trying to implement the srivastava and sundaresan model for friction pressure and viscosity via UDF in fluent. My problem consist in a fluidized bed where air is injected through a jet. Johnson&Jsckson and Schaeffer models, available in fluent, failed to predict the bubbling of the bed. I am now trying to implement this model, but the UDF I write causes the divergence of the solver. You can find the UDF code below. Is there anyone who knows what's wrong?

DEFINE_PROPERTY(frictional_pressure, cell_s, thread_s)
{

real es, dudx, dudy, dvdx, dvdy, div, S11, S12, S21, S22, SS, theta, pc, pf, n, den, exponent, Fr;


/* find velocity gradient for solid phase, solid volume fraction and granular temperature */

es = C_VOF(cell_s, thread_s); /*solid volume fraction */

dudx = C_DUDX(cell_s, thread_s); /* velocity gradient components */

dudy = C_DUDY(cell_s, thread_s); /* velocity gradient components */

dvdx = C_DVDX(cell_s, thread_s); /* velocity gradient components */

dvdy = C_DVDY(cell_s, thread_s); /* velocity gradient components */

theta = C_GT(cell_s, thread_s); /* granular temperature 8/

/* define divergence and compute S:S */

div = dudx+dvdy; /* divergence of solid velocity */

S11 = 2/3*dudx-1/3*dvdy; /* component 11 of the strain rate tensor */

S12 = 0.5*(dudy+dvdx); /* component 12 of the strain rate tensor */

S21 = S12; /* component 21 of the strain rate tensor */

S22 = 2/3*dvdy-1/3*dudx; /* component 22 of the strain rate tensor */

SS = abs(pow(S11,2)+2*S12+pow(S22,2)); /* S:S */

/* compute pc, Johnson&Jacskon model */

Fr = 0.1*es;

if(es > esmin)

pc = Fr*(pow((es-esmin),r)/pow((esmax-es),s));

else
pc = 0;

/* compute n */

if(div >= 0)

n = sqrt(3)/(2*sin(fi));

else
n = 1.03;

/* compute and return pf */

den = n*sqrt(2)*sin(fi)*sqrt(SS+theta/(pow(ds,2)));

exponent = 1/(n-1);

if (den == 0)

pf = pc;

else

pf = pc*(pow((1-div/den),exponent));

return pf;
}

DEFINE_PROPERTY(frictional_viscosity, cell_s, thread_s)
{

real es, dudx, dudy, dvdx, dvdy, div, S11, S12, S21, S22, SS, theta, pc, pf, n, den1, den2, exponent, Fr, muf;


/* find velocity gradient for solid phase, solid volume fraction and granular temperature */

es = C_VOF(cell_s, thread_s); /*solid volume fraction */

dudx = C_DUDX(cell_s, thread_s); /* velocity gradient components */

dudy = C_DUDY(cell_s, thread_s); /* velocity gradient components */

dvdx = C_DVDX(cell_s, thread_s); /* velocity gradient components */

dvdy = C_DVDY(cell_s, thread_s); /* velocity gradient components */

theta = C_GT(cell_s, thread_s); /* granular temperature 8/

/* define divergence and compute S:S */

div = dudx+dvdy; /* divergence of solid velocity */

S11 = 2/3*dudx-1/3*dvdy; /* component 11 of the strain rate tensor */

S12 = 0.5*(dudy+dvdx); /* component 12 of the strain rate tensor */

S21 = S12; /* component 21 of the strain rate tensor */

S22 = 2/3*dvdy-1/3*dudx; /* component 22 of the strain rate tensor */

SS = abs(pow(S11,2)+2*S12+pow(S22,2)); /* S:S */

/* compute pc, Johnson&Jacskon model */

Fr = 0.1*es;

if(es > esmin)

pc = Fr*(pow((es-esmin),r)/pow((esmax-es),s));
else
pc = 0;

/* compute n */

if(div >= 0)

n = sqrt(3)/(2*sin(fi));
else
n = 1.03;

/* compute pf */

den1 = n*sqrt(2)*sin(fi)*sqrt(SS+theta/(pow(ds,2)));

exponent = 1/(n-1);

if (den1 == 0)

pf = pc;

else

pf = pc*(pow((1-div/den1),exponent));

/* compute and return muf */

den2 = sqrt(SS+theta/(pow(ds,2)));

if((den2 == 0) && (pc != 0))

muf = pf*sin(fi);

if (pc == 0)

muf = 0;

else

muf = sqrt(2)*pf*sin(fi)/den2*(n-(n-1)*pow(pf/pc, exponent));

return muf;

}

[tarzy]

Shouldn't be the in the S:S calculation:
pow(S11,2)+2*pow(S12,2)+pow(S22,2)

instead of
pow(S11,2)+2*S12+pow(S22,2) ?