uds in multiphase problem
 

hi,everyone.
I want to solve a scalar eqution that donot have the voidage in the whole eqution for one phase in multiphase (note that uds in multiphas fluent will auto multipy the phase's voidage in both unsteady term and convection term), so I decide to write a udf to change the default, but one problem is that the udf looks do not work fine ,because I find the max value is inceasing with time ,that should not happened! Can anybody help me to check my uds? or tell me why the max value is increasing during the iterate? thanks!
the following is my uds

DEFINE_UDS_UNSTEADY(my_uds_unsteady,c,t,i,apu,su)
{
real physical_dt, vol, rho, phi_old;
physical_dt = RP_Get_Real("physical-time-step");
vol = C_VOLUME(c,t);
rho = C_R_M1(c,t);
*apu = -rho*vol / physical_dt ;/*implicit part*/
phi_old = C_STORAGE_R(c,t,SV_UDSI_M1(i));
*su = rho*vol*phi_old/physical_dt ;/*explicit part*/
}

DEFINE_UDS_FLUX(my_uds_flux,f,t,i)
{
cell_t c0, c1 = -1;
Thread *t0, *t1 = NULL;

real NV_VEC(psi_vec), NV_VEC(A), flux = 0.0;

c0 = F_C0(f,t);
t0 = F_C0_THREAD(f,t);
F_AREA(A, f, t);

/* If face lies at domain boundary, use face values; */
/* If face lies IN the domain, use average of adjacent cells. */

if (BOUNDARY_FACE_THREAD_P(t)) /*Most face values will be available*/
{
real dens;

/* Depending on its BC, density may not be set on face thread*/
if (NNULLP(THREAD_STORAGE(t,SV_DENSITY)))
dens = F_R(f,t); /* Set dens to face value if available */
else
dens = C_R(c0,t0); /* else, set dens to cell value */

NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, dens);

flux = NV_DOT(psi_vec, A); /* flux through Face */
}
else
{
c1 = F_C1(f,t); /* Get cell on other side of face */
t1 = F_C1_THREAD(f,t);

NV_DS(psi_vec, =, C_U(c0,t0),C_V(c0,t0),C_W(c0,t0),*,C_R(c0,t0));
NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,C_R(c1,t1));

flux = NV_DOT(psi_vec, A)/2.0; /* Average flux through face */
}

return flux;
}