[moloykb]

Dear all,
I am trying to simulate a problem where the transport of drug has to determined within blood using some external magnet, which will guide the drug towards the target location. In order to solve the issue, there are two different forces that are present one is the magnetic force (Fmag) and other is the drag force (Fdrag). In order to implement the effect of Fmag, I include two different source term for both x and y momentum equation. But for the UDS equation I have to define a velocity say, vp=vblood+Fmag/(Coeff of Drag), with ihis vp I have to solve the convective term of the UDS equation. Following is the UDF that I have written for this, It compiles perfectly well but once I try to solve it says "divergence detected in the solver: uds0"
bellow is what I have written:

#include "udf.h"

#define rp 1.25e-7
#define eta 3.45e-3
#define dens 1.056e3
DEFINE_UDS_FLUX(drug_flux,f,t,i)
{
real muzero00,Mag_const00,muzero01,Mag_const01;
real x,y;
real Numer00,Numer01,Denom01,Denom11,Denom02,Denom12,De nom00;
real bforcex, bforcey;
real vmagx, vmagy;
real xc[ND_ND];
real vmag[ND_ND], vmag1[ND_ND];
real vb[ND_ND], vb0[ND_ND], vb1[ND_ND];

cell_t c0,c1 = -1;
Thread *t0,*t1 = NULL;

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

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


C_CENTROID(xc,c0,t0);
x=xc[0];
y=xc[1];

muzero00=4.*PI*(pow(10.,(-7.)));
Mag_const00=-2*pow(p,2)*muzero00*Kieff;
Numer00=x-xmag;
Denom01=(pow(y-ymag,(2.)));
Denom02=(pow(x-xmag,(2.)));
Denom00=(pow(Denom01+Denom02,(3.)));
bforcex=Mag_const00*Numer00/Denom00;

muzero01=4.*PI*(pow(10.,(-7.)));
Mag_const01=-2*pow(p,2)*muzero01*Kieff;
Numer01=y-ymag;
Denom11=(pow(y-ymag,(2.)));
Denom12=(pow(x-xmag,(2.)));
Denom00=(pow(Denom11+Denom12,(3.)));
bforcey=Mag_const01*Numer01/Denom01;

vmagx=bforcex/(6*PI*eta*rp);
vmagy=bforcey/(6*PI*eta*rp);


vmag[0]=vmagx;
vmag[1]=vmagy;
vmag[2]=0.0;

if(BOUNDARY_FACE_THREAD_P(t))
{
//real dens;
//if(NULLP(THREAD_STORAGE(t,SV_DENSITY)))
// dens=F_R(f,t);
//else
//dens=C_R(c0,t0);

vb[0]=F_U(f,t);
vb[1]=F_V(f,t);
vb[2]=0.0;

NV_VV(psi1_vec, =, vmag, +, vb);
NV_S(psi1_vec, *=, dens);

flux = NV_DOT(psi1_vec, A);
}
else
{
c1=F_C1(f,t);
t1=F_C1_THREAD(f,t);

C_CENTROID(xc,c1,t1);
x=xc[0];
y=xc[1];

muzero00=4.*PI*(pow(10.,(-7.)));
Mag_const00=-2*pow(p,2)*muzero00*Kieff;
Numer00=x-xmag;
Denom01=(pow(y-ymag,(2.)));
Denom02=(pow(x-xmag,(2.)));
Denom00=(pow(Denom01+Denom02,(3.)));
bforcex=Mag_const00*Numer00/Denom00;

muzero01=4.*PI*(pow(10.,(-7.)));
Mag_const01=-2*pow(p,2)*muzero01*Kieff;
Numer01=y-ymag;
Denom11=(pow(y-ymag,(2.)));
Denom12=(pow(x-xmag,(2.)));
Denom00=(pow(Denom11+Denom12,(3.)));
bforcey=Mag_const01*Numer01/Denom01;

vmagx=bforcex/(6*PI*eta*rp);
vmagy=bforcey/(6*PI*eta*rp);

vmag1[0]=vmagx;
vmag1[1]=vmagy;

vb0[0]=C_U(c0,t0);
vb0[1]=C_V(c0,t0);

vb1[0]=C_U(c1,t1);
vb1[1]=C_V(c1,t1);

NV_VV(psi1_vec, =,vmag, +, vb0);
NV_VV(psi2_vec, =,vmag1, +, vb1);
NV_VS_VS(psi_vec, =, psi1_vec, *, dens, +, psi2_vec, *, dens);

flux = NV_DOT(psi_vec, A)/2.0;
}

return flux;
please help me to solve this issue...

[Rahul123]

could you please provide the equation for the Fmag that you are considering as the magnetophoretic body force