[Achini]

Hello!

The following code is for interaction forces between particle in a rectangular channel.
The expression I'm using for particle wall interactions is, F = K/(dw^2), which is the force exerted on the particle by a wall.
so, when the particle is placed in the channel, there will be 4 forces exerted by the 4 channel walls. And the net force can be calculated.

I have considered x,y,z components separately, and created 2 face loops to loop over faces parallel to xy plane and yz plane separately.
And I have used a face loop to get the minimum distance between the particle and the wall.

The code runs without any error.
But the particle tracks look different.

Can someone help me with this?
Thank you!

#include "udf.h"
#include "dpm.h"
#include "dpm_mem.h"
#include "surf.h"

int particleid;

DEFINE_DPM_BODY_FORCE(particle_body_force,p,i) /* Fluent macro to define particle body force*/
{
double bf,dw,area,facearea,e,K;
real xw[ND_ND];
real bforce[ND_ND];
real vec[ND_ND];

Domain*d;
cell_t c;
Thread *t;
face_t f;

d = Get_Domain(1);

K = 4.00365*(10^(-24));


particleid = (int)p->part_id;

if (particleid % 2 == 0) return(0); /*This part makes it such that only half of the particles get the force. */

c = P_CELL(p); /*Get the cell and thread that the particle is currently in */

/* t = P_CELL_THREAD(p);*/

if (i==0) /* to consier the direction */
{
t = Lookup_Thread(d, 7); /* channel wall on yz plane */

NV_S(bforce, =, 0.0); /* reset force*/

dw=200;
begin_f_loop(f, t)
{


F_CENTROID(xw,f,t); /* Center of the face*/


NV_VV(vec, =, xw, -, P_POS(p)); /* vector conecting the particlea nd the face center*/


if (dw>NV_MAG(vec))
{
dw = NV_MAG(vec);

/*Cos(Q) = vec[i]/dw;*/
}



end_f_loop(f, t)

}

bf = (K/(vec[i]* vec[i]))-(K/((150-vec[i])*(150-vec[i])));
/* Force from the infinitesimal face area*/
bforce[i] = bf; /* force in the relavent cartesian direction; vec[i]/dw is the cos(theta)*/
}

if (i==1) /* to consier the direction */
{
t = Lookup_Thread(d, 8); /* channel wall on xy plane*/

NV_S(bforce, =, 0.0); /* reset force*/

dw=200;
begin_f_loop(f, t)
{


F_CENTROID(xw,f,t); /* Center of the face*/


NV_VV(vec, =, xw, -, P_POS(p)); /* vector conecting the particlea nd the face center*/


if (dw>NV_MAG(vec))
{
dw = NV_MAG(vec);

}



end_f_loop(f, t)

}

bf = (K/(vec[i])* (vec[i]))-(K/((150-vec[i])*(150-vec[i]))); /* calculation of the net force by attractions from walls on yz plane*/

bforce[i] = bf;
}

if (i==2)
{
bforce[i]=0;
}

return (bforce[i]/P_MASS(p));


}