[ustbdynamic]

Hi friends.
I am trying to use two UDSs solve the following equations. (uds0 for eq1 and uds1 for eq2, UDS2 and UDS3 are used to storage the uds0’s Gradient in x direction and y direction respectively ).
It seems that the problem is very easy, but I haven’t got it long time. When the UDS1’s source is added, the solution cannot converge. While without it, the simulation works well. I found that the used UDS’s Gradient is very high in the boundary zone. It may be the cause.
Can somebody please help me.

thank you very much.

en.jpg

#include "udf.h"
#define B_ION 1.5e-4 /* Ion mobility - positive ion estimate */
#define EPS_0 8.85418717e-12

DEFINE_ADJUST(Calc_E_xy,d)
{
Thread *t;
cell_t c;
face_t f;
Thread *t0;
cell_t c0;
real dr0[ND_ND],dr1[ND_ND], xf[ND_ND], dy;
/* Do nothing if gradient isn't allocated yet. */
if (!Data_Valid_P()) return;
thread_loop_c(t,d)
{
if (NULL != THREAD_STORAGE(t,SV_UDS_I(0)) && NULL != T_STORAGE_R_NV(t,SV_UDSI_G(0)))
{
begin_c_loop(c,t)
{
C_UDSI(c,t,2) = -1.*C_UDSI_G(c,t,0)[0];
C_UDSI(c,t,3) = -1.*C_UDSI_G(c,t,0)[1];
}
end_c_loop(c,t)
}
}
thread_loop_f(t,d)
{
if (NULL != THREAD_STORAGE(t,SV_UDS_I(3)) && NULL != T_STORAGE_R_NV(t->t0,SV_UDSI_G(3)))
{
if (BOUNDARY_FACE_THREAD_P(t))
{
begin_f_loop(f,t)
{
c0 = F_C0(f,t);
t0 = THREAD_T0(t);
F_CENTROID(xf,f,t);
C_CENTROID(dr0,c0,t0);
NV_VV(dr1,=,xf,-,dr0);
F_UDSI(f,t,2) = C_UDSI(c0,t0,2) + (C_UDSI_G(c0,t0,2)[0]*dr1[0]);
F_UDSI(f,t,3) = C_UDSI(c0,t0,3) + (C_UDSI_G(c0,t0,3)[1]*dr1[1]);
}
end_f_loop(f,t)
}
}
}
}

/************************************************** *****************/
/* UDF_source_phi_rho_11.c for specifying the phi and rho source terms for uds-0 */
/************************************************** *****************/
DEFINE_SOURCE(phi_source_temp,c,t,dS,eqn)
{
real phi_source;
dS[eqn] = 0.0;
phi_source = C_UDSI(c,t,1)/EPS_0;
return phi_source;
}
DEFINE_SOURCE(rho_source_temp,c,t,dS,eqn)
{
real rho_source;
real xc[ND_ND];
real dens;
//dS[eqn] = 0;
C_CENTROID(xc,c,t);
dens=C_R(c,t);
rho_source = C_UDSI(c,t,1)*B_ION*(0 + C_UDSI_G(c,t,3)[1]);//C_UDSI_G(c,t,2)[0]
rho_source = -rho_source;
return rho_source;
}

DEFINE_UDS_FLUX(rho_ion_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 (BOUNDARY_FACE_THREAD_P(t))
{

NV_DS(psi_vec, =, F_U(f,t), F_V(f,t), F_W(f,t), *, 1);
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),*,1.0);
NV_DS(psi_vec, +=, C_U(c1,t1),C_V(c1,t1),C_W(c1,t1),*,1.0);
flux = NV_DOT(psi_vec, A)/2.0; /* Average flux through face */
}
return flux;
}