DPM Body Force Needs Help.
 

Hello Everyone,

I try to define the body force under DPM model, however, even though I can complie the code, it seems like that my body force always equal to zero. The goal for this udf code is try to calcuate boday force only act on two particles.

Here is my UDF:

#include "udf.h"
#include "mem.h"
#include "dpm.h"
#include "surf.h"
#include "stdio.h"
#define A 3.64e-20 /*Hamaker constant*/
#define KD 6.89e+12 /*Debye-Huckel Parameter*/
#define R 8.314 /*Idea gas constant*/
#define Ve 2.30e+15 /*Valance of the electrolyte*/
#define F 96485.3365 /*Faraday constant*/
#define Z1 -0.0517 /*Zeta potential of particle*/
#define Z2 -0.02 /*Zeta potential of substrate*/
#define T 300 /*Thermaldynamic temperature at Kelvin*/
#define er 85.0 /*Dielectric constant of the medium*/
#define eo 8.85e-12 /*Permittivity of vacuum*/
#define M_PI 3.14159265358979323846

DEFINE_DPM_BODY_FORCE(particle_body_force,p,i)
{
Domain *d;
Thread *t;
cell_t c;
real bforce;
real vdwF;
real eleF;
real d1;
real d2;
real m;
real s;
real h;
real uex;
real uey;
real OP;
real op;

Particle *pi; /*loop over all particles in all cells and find their position and diameter*/
int ii;
real p1[2], p2[2];

d=Get_Domain(1);
for(ii = 0; ii<2; ii++)
{
p1[ii]= p->state.pos[ii];
{
d1 = P_DIAM(p);
m = P_MASS(p);
thread_loop_c(t,d) /*loops over all cell threads in domain*/
{
begin_c_loop(c,t) /*loops over cells in a cell thread*/
}

begin_particle_cell_loop(pi,c,t)
{
for (ii=0;ii<2;ii++)
{
p2[ii]= p->state.pos[ii];
}
d2 = pi->state.diam;
s = sqrt(pow(p1[0]-p2[0],2)+pow(p1[1]-p2[1],2)); /*particle center to center distance*/
h = s-((d1+d2)/2); /*particle surface separation distance*/
OP = pow((s*s-d1*d1-2*d1*d2-d2*d2),2);
op = pow((s*s-d1*d1-2*d1*d2-d2*d2+4*d1*d2),2);
if (s>=(d1/2+d2/2)||h>=0)
{

uex = ((p1[0]-p2[0])/s);
uey = ((p1[1]-p2[1])/s);
vdwF = (-1)*A*64*pow((d1*d2),3)*s/(6*OP*op);
eleF = 64*M_PI*er*eo*KD*pow((R*T/(Ve*F)),2)*tanh(Ve*F*Z1/(4*R*T))*tanh(Ve*F*Z2/(4*R*T))*(d1*d2)*exp((-1)*KD*h)/(d1+d2);
if (i==0)
bforce=(vdwF*uex+eleF*uex) ;
else if(i==1)
bforce=(vdwF*uey+eleF*uey) ;

}

}
end_particle_cell_loop(pi,c,t)
}

end_c_loop(c,t)

} /*end of thread_loop_c(t,d)*/

return (bforce/m);
}


I think the problem is i failed to update the particle information at each loop.

So, please help me out and thank you in advance.!!

Best Regards