uds's flux term in multiphase flow
 

dear all
I tried to use DEFINE_UDS_FLUX to return the default flux term (UA) in multiphase(Eulerian model) flow, but the calculated scalar is out of range, i.e. greadter than 1 even if the inlet value is 1.　I have tried first-order upwind and central-differencing scheme but the scalars were still out of range.
so can anyone help? Thanks at first!

Below are my codes.

#include "udf.h"

DEFINE_UDS_FLUX(cd, face, thread_g, i)
{
cell_t c0, c1 = -1;
Thread *tg0, *tg1 = NULL;

real NV_VEC(A), flux = 0.0;

real flux0 = 0.0, flux1 = 0.0, NV_VEC(psi_vec0), NV_VEC(psi_vec1);

real rho_g0, vof_g0, vof_g1;


c0 = F_C0(face, thread_g);
tg0 = F_C0_THREAD(face, thread_g);
F_AREA(A, face, thread_g);

if (BOUNDARY_FACE_THREAD_P(thread_g))
{
if (NNULLP(THREAD_STORAGE(thread_g, SV_DENSITY)))
{
rho_g0 = F_R(face, thread_g);
}
else
{
rho_g0 = C_R(c0, tg0);
}

NV_DS(psi_vec0, =, F_U(face, thread_g), F_V(face, thread_g), F_W(face, thread_g), *, rho_g0);
flux = NV_DOT(psi_vec0, A) * F_VOF(face, thread_g);
}
else
{
c1 = F_C1(face, thread_g);
tg1 = F_C1_THREAD(face, thread_g);

vof_g0 = C_VOF(c0, tg0);
vof_g1 = C_VOF(c1, tg1);

NV_DS(psi_vec0, =, C_U(c0, tg0), C_V(c0, tg0), C_W(c0, tg0), *, C_R(c0, tg0));
NV_DS(psi_vec1, =, C_U(c1, tg1), C_V(c1, tg1), C_W(c1, tg1), *, C_R(c1, tg1));

flux0 = NV_DOT(psi_vec0, A) * vof_g0;
flux1 = NV_DOT(psi_vec1, A) * vof_g1;
flux = 0.5 *(flux0 +flux1);
}
return flux ;
}