Learning Forum

[PrantikDas]

I have written the UDF for Burns et. al Model: -

# include "udf.h"

DEFINE_VECTOR_EXCHANGE_PROPERTY(burnetal,cell,mixture,liquid,solid,Ftdvector)

{

    Thread *thread_l, *thread_s;

    real rho_l,rho_s,Sigmapq,alpha_l,alpha_s,term,Ftd,Ctd,Kpq;

    real px,py,pz,sx,sy,sz;

    real x_vel_l,y_vel_l,z_vel_l,x_vel_s,y_vel_s,z_vel_s,slip_x,slip_y,slip_z,slip;

    thread_l=THREAD_SUB_THREAD(mixture,liquid);

    thread_s=THREAD_SUB_THREAD(mixture,solid);

    /* Calculation of Gradients */

    alpha_l=C_VOF(cell, thread_l); /* liquid phase volume fraction */

    alpha_s=C_VOF(cell, thread_s); /* Solid Phase Volume fraction */

    rho_l=C_R(cell,thread_l);

    rho_s=C_R(cell,thread_s);

    Ctd=1;

    Sigmapq=.9;

    Kpq=alpha_l*rho_l*C_K(cell, thread_l)+alpha_s*rho_s*C_K(cell, thread_s);

    term=Ctd*Kpq*C_MU_T(cell,thread_l)/(rho_l*Sigmapq);

    px=C_VOF_G(cell,thread_l)[0]/alpha_l;

    py=C_VOF_G(cell,thread_l)[1]/alpha_l;

    pz=C_VOF_G(cell,thread_l)[2]/alpha_l;

    sx=C_VOF_G(cell,thread_s)[0]/alpha_s;

    sy=C_VOF_G(cell,thread_s)[1]/alpha_s;

    sz=C_VOF_G(cell,thread_s)[2]/alpha_s;

    /* Calculation of the Final Expression */

    Ftdvector[0]=term*(sx-px);

    Ftdvector[1]=term*(sy-py);

    Ftdvector[2]=term*(sz-pz);

}   

But I am getting a floating-point error after certain iterations and also unable to know why the formula used by ANSYS THEORY GUIDE for the Burns et. al model (Turbulent Dispersion) is different from that used in the original paper.

[Rahul Kumar]

Hello,
Are you sure that the floating point error is because of the UDF? How do the monitors look like before the error?

[PrantikDas]

#include"udf.h"
DEFINE_VECTOR_EXCHANGE_PROPERTY(burnetal,cell,mixture,liquid,solid,Ftdvector)
{
Thread*thread_l,*thread_s;
realx_vel_l,y_vel_l,z_vel_l,x_vel_s,y_vel_s,z_vel_s,slip_x,slip_y,slip_z,slip,rho_l,Cd,Sigmapq,diam,alpha_l,alpha_s,term,rey;
realpx,py,pz,sx,sy,sz,mu_e,tau_p,D,K,mut,f,mu_l,k,e;
thread_l=THREAD_SUB_THREAD(mixture,liquid);
thread_s=THREAD_SUB_THREAD(mixture,solid);
/*Definingtheslipvelocity*/
x_vel_l=C_U(cell,thread_l);
y_vel_l=C_V(cell,thread_l);
z_vel_l=C_W(cell,thread_l);
x_vel_s=C_U(cell,thread_s);
y_vel_s=C_V(cell,thread_s);
z_vel_s=C_W(cell,thread_s);
slip_x=x_vel_l-x_vel_s;
slip_y=y_vel_l-y_vel_s;
slip_z=z_vel_l-z_vel_s;
slip=sqrt(pow(slip_x,2)+pow(slip_y,2)+pow(slip_z,2));
rho_l=C_R(cell,thread_l);
Sigmapq=0.9;
k=C_K(cell,thread_l);
e=C_D(cell,thread_l);
mut=.09*pow(k,2)/e;
mu_l=C_MU_L(cell,thread_l);
diam=C_PHASE_DIAMETER(cell,thread_s);
alpha_l=C_VOF(cell,thread_l);
alpha_s=C_VOF(cell,thread_s);
mu_e=C_MU_EFF(cell,thread_l);
/*DefiningReynoldsNumber*/
rey=(rho_l*slip*diam)/mu_l;
/*CalculationofCd*/
if(rey==0)
{
Cd=0;
}
elseif(rey<1000)
{
Cd=((24/(rey*alpha_l))*(1+(0.15*pow((alpha_l*rey),0.687))));
}
else
{
Cd=0.44;
}
/*CalculationofGradients*/
D=mut/(rho_l);
f=Cd*rey/24;
tau_p=C_R(cell,thread_s)*pow(diam,2)/(18*mu_l);
K=(C_R(cell,thread_s)*f*alpha_s*alpha_l)/(tau_p);
term=K*1*D/Sigmapq;
px=C_VOF_G(cell,thread_l)[0]/alpha_l;
py=C_VOF_G(cell,thread_l)[1]/alpha_l;
pz=C_VOF_G(cell,thread_l)[2]/alpha_l;
sx=C_VOF_G(cell,thread_s)[0]/alpha_s;
sy=C_VOF_G(cell,thread_s)[1]/alpha_s;
sz=C_VOF_G(cell,thread_s)[2]/alpha_s;
/*CalculationoftheFinalExpression*/
Ftdvector[0]=term*(sx-px);
Ftdvector[1]=term*(sy-py);
Ftdvector[2]=term*(sz-pz);
printf("%f\t%f\n",k,e);
}

I rewrote the code as previously, I considered Kpq as Mixture Turbulent kinetic energy but it is the momentum exchange term. The turbulent viscosity is causing the divergence. I am getting floating point error due to that.

The residuals are given above. I am continuously getting the error :-#I have a message "Turbulent viscosity limited to viscosity ratio of 1.000000e+05 in xxx cells. How can I fix this error message in Ansys fluent? So I think the problem is with the formula. The formula given in the original paper is quite different from what is given in Ansys Theory Guide for Burns et al Model.
By the way, I am coding according to ANSYS THEORY GUIDE.

[DrAmine]

Why are you coding what is already available in the code? Does it run with the built-in Burns et al?

[PrantikDas]

I am using the formula given in ANSYS THEORY GUIDE. I don't know whether it is the formula used in the built-in Burns et al Model. I am coding it so that I can compare the order of magnitude for lift, drag, and turbulent dispersion force for local cells.

[DrAmine]

Okay I understand. What you can do is avoid this kind of unbounded division (division by alpha).
I have the formula in the original paper and is the same as used in our codes as the author is an Ansys colleague.
I think the vector which should be returned is relative to the primary phase like in the example of Lopez de Bertodano Case.

But before you do that: have you verified running with the Burns. et al. without UDF?

[PrantikDas]

Yes, and it is working fine. Can you provide me the correct UDF?

[DrAmine]

Unfortunately, I cannot provide you a correct UDF.

[PrantikDas]

Can you provide me the correct formula so that I can incorporate it in my UDF? or any advice to make my code workable?

[Rob]

What error is in the Theory Guide?
General advice is to comment code and indent to make it easier to see what section is linked to and/or does what. An editor like Notepad++ or the like that recognises C code is also helpful as it'll tend to colour text based on what it's doing, eg DEFINE lines may be coloured differently to REAL etc.

[PrantikDas]

#include"udf.h"
DEFINE_VECTOR_EXCHANGE_PROPERTY(burnetal,cell,mixture,liquid,solid,Ftdvector)
{
Thread*thread_l,*thread_s;
realx_vel_l,y_vel_l,z_vel_l,x_vel_s,y_vel_s,z_vel_s,slip_x,slip_y,slip_z,slip,rho_l,Cd,Sigmapq,diam,alpha_l;
real alpha_s,term,rey;
realpx,py,pz,sx,sy,sz,tau_p,D,K,mu_l,Kpq;
realCtd=1;
thread_l=THREAD_SUB_THREAD(mixture,liquid);
thread_s=THREAD_SUB_THREAD(mixture,solid);
/*Definingtheslipvelocity*/
x_vel_l=C_U(cell,thread_l);
y_vel_l=C_V(cell,thread_l);
z_vel_l=C_W(cell,thread_l);
x_vel_s=C_U(cell,thread_s);
y_vel_s=C_V(cell,thread_s);
z_vel_s=C_W(cell,thread_s);
slip_x=x_vel_l-x_vel_s;
slip_y=y_vel_l-y_vel_s;
slip_z=z_vel_l-z_vel_s;
slip=sqrt(pow(slip_x,2)+pow(slip_y,2)+pow(slip_z,2));
/*Densityoftheliquid*/
rho_l=C_R(cell,thread_l);
/*DefiningsigmausedinBurnsetalformulainANSYSTHEORYGUIDE*/
Sigmapq=0.9;
/*Defininglaminarviscosity,diameterandvolumefraction*/
mu_l=C_MU_L(cell,thread_l);
diam=C_PHASE_DIAMETER(cell,thread_s);
alpha_l=C_VOF(cell,thread_l);
alpha_s=C_VOF(cell,thread_s);
/*DefiningReynoldsNumber*/
rey=(rho_l*slip*diam)/mu_l;
/*CalculationofCdbyusinggidaspowmodel*/
/*DeterminationofmomentumexchangecoefficientbygidaspowmodelforsolidliquidmomentumexchangeandtheformulahasbeentakenfromANSYSTHEORYGUIDE*/
if(alpha_l>.8)
{
if(rey==0)
{
Cd=0;
}
else
{
Cd=(24*pow((rey*alpha_l),-1)*(1+(0.15*pow((alpha_l*rey),0.687))));
}
Kpq=(.75*Cd*alpha_l*alpha_s*rho_l*slip*pow(alpha_l,-2.65))*pow((diam*4),-1);
}
else
{
Kpq=(150*(alpha_s*(1-alpha_l)*mu_l)*pow(alpha_l*pow(diam,2),-1))+(1.75*rho_l*alpha_s*slip*pow(diam,-1));
}
/*Definingthedispersionscalar*/
D=C_MU_T(cell,thread_l)*pow(rho_l,-1);
/*Terminfrontofgradient*/
term=Ctd*Kpq*D*pow(Sigmapq,-1);
/*CalculationofGradients*/
px=C_VOF_G(cell,thread_l)[0]/alpha_l;
py=C_VOF_G(cell,thread_l)[1]/alpha_l;
pz=C_VOF_G(cell,thread_l)[2]/alpha_l;
sx=C_VOF_G(cell,thread_s)[0]/alpha_s;
sy=C_VOF_G(cell,thread_s)[1]/alpha_s;
sz=C_VOF_G(cell,thread_s)[2]/alpha_s;
/*CalculationoftheFinalExpressionandtheformulahasbeentakenfromANSYSTHEORYGUIDE*/
Ftdvector[0]=term*(sx-px);
Ftdvector[1]=term*(sy-py);
Ftdvector[2]=term*(sz-pz);
}
I have commented on the code though I am getting the same error.

[Rob]

What is the flow field doing a few iterations or time steps before the error?

[DrAmine]

The formula is the one mentioned in the Theory Guide. UDF Manual describes an example for Lopez de Bertodano Dispersion.

[PrantikDas]

I have used UDMI to store the values:-
C_UDMI(cell,mixture,0)=D;
C_UDMI(cell,mixture,1)=term;
C_UDMI(cell,mixture,2)=Ftdvector[0];
C_UDMI(cell,mixture,3)=Ftdvector[1];
C_UDMI(cell,mixture,4)=Ftdvector[2];
and I have multiplied the term by cell volume as given in Lopez Model by using the command:-
term*=C_VOLUME(cell,mixture);
I am getting a floating-point error at the 32nd iteration and the contours at the inlet at the 31st iteration are shown above. These are the contours of various variables at the inlet of my pipe (The variables which are stored in my UDMI)

[PrantikDas]

Another question I have
in Lopez model sample UDF given in ANSYS UDF MANUAL what does this lines mean ?
real turb_k = (mp_ke_type==TURB_MP_KE_MIXTURE ||
mp_rsm_type==TURB_MP_RSM_MIXTURE)?
C_K(c,t) : C_K(c,pt[P_PHASE]);
2.if(NNULLP(THREAD_STORAGE(tj,SV_MP_DRIFT_S_P_COEFF)))
P_DRIFT_COEFF(c,tj) = 0.0;
3.if(NNULLP(THREAD_STORAGE(tj,SV_MP_DRIFT_COEFF)))
S_DRIFT_COEFF(c,tj) = term;

[DrAmine]

The documentation goes in some aspects into what you are asking.

1/Getting TKE using some sanity checks. C programming + some Fluent Predefined variables.
2/The two others terms are not documented but what I can understand is that the second one requires your return value for the drift coefficient and the other ones is not required as you are not adding term into the continuity equation.

I can double check the meaning of those terms. I will write back if I have something different to share but I am confident what I wrote is correct.

[PrantikDas]

Any suggestions regarding my UDF code and contour?

[DrAmine]

No I do not have suggestions for the UDF apart from what I have shared already and to rely on the semantic from the manual (force to be returned is the one exerted on the primary phase)

[Rob]

Now look at the flow field, this requires more than looking at just the values on the inlet/outlet/walls. Something happens in the domain, and this may help guide you to a solution.

[DrAmine]

P_DRIFT_COEFF is the coefficient of the volume fraction gradient of the primary phase (zero for Lopez et. al.) S_DRIFT_COEFF the one for dispersed volume fraction gradient.

[Sofen Kumar Jena]

What is the final expression for slip velocity for these models ? How the volumetric force is converted into slip velocity term ?