particles wall interaction
 

hello guys i created this UDF that works without error but i see that when the particle in injected it is not affected by the lateral force and i do not know why.
I post the UDF.
Thanks in advance


/* UDF adding an external body force to the particles */
#include "udf.h"
#include "dpm.h"
#include "dpm_mem.h"
#include "surf.h"


DEFINE_DPM_BODY_FORCE(particle_body_force, p, i)
{

/* declaration of variables */
double w, Dh, Ufx, Ufy, Ufz, Gx, rho, a, F_L, c_height, c_length,
c_volume, side, H, W, f_height_total, Renx, Reny, Renz, crDh, mu, aUfx, aUfy, aUfz, C1, C2, C3, C4, C_Ly, C_Lz, z, y, distancewall, cosine;



real xw[ND_ND];
real vec[ND_ND];
Domain* d;
cell_t c;
Thread* t;
face_t f; //f is the face index that is inside the face thread t(our wall). You will extract the face area vector from f index
d = Get_Domain(1);
c = P_CELL(p); //cell in witch the particle is currently in
t = Lookup_Thread(d, 9); //9 is the ID number of the upper wall. t is the thread index(pointer) of the wall boundary surface 9



distancewall = 100; //i'm setting an generic initial value
begin_f_loop(f, t)
{
F_CENTROID(xw, f, t); //center of each face
NV_VV(vec, =, xw, -, P_POS(p)); //vector connecting the particle to the face center that is placed at the wall. The NV_VV do the operation between 2 vectors :(vec=xw-P_POS(p))




if (distancewall >= NV_MAG(vec))
{
distancewall = NV_MAG(vec);
}
}
end_f_loop(f, t)




//i'm doing a double cicling of all faces of the thread. I start with distancewall=100 and then i compare it with the distance of the particle respect to a certain face centroid\\
//if distancewall is bigger than that value, you fix it equal to them and then you do again a comparin but with the next face centroid. if is not more bigger you close the loop
//and you get the minimum distance.


Thread* tt = P_CELL_THREAD(p);



H = 0.00005; /*height of the cross section of the channel */
W = 0.0001; /*width of the cross section of the channel */
Dh = 0.00006667; /*hydraulic diameter of the channel */
mu = 0.001003; /*viscosity of the fluid */
rho = C_R(c, tt);
Ufx = C_U(c, tt);
Ufy = C_V(c, tt);
Ufz = C_W(c, tt);
aUfx = fabs(Ufx); /* absolute value of Ufx */
aUfy = fabs(Ufy); /* absolute value of Ufy */
aUfz = fabs(Ufz); /* absolute value of Ufz */


/* local Reynolds number calculation*/

Renx = (rho * aUfx * Dh) / mu;
Reny = (rho * aUfy * Dh) / mu;
Renz = (rho * aUfz * Dh) / mu;


Gx = (Ufx) / H; /*max shear rate of the channel*/

a = P_DIAM(p); /* particle diameter */

y = fabs(((H / 2) - distancewall) / (H / 2));

C1 = -3.0374;
C2 = -0.3046;
C3 = 0.9790;
C4 = -0.0016;

C_Ly = C1 * pow(y, 3) + C2 * pow(y, 2) + C3 * pow(y, 1) + C4 * pow(y, 0); //valido nel primo e secondo quadrante

if (i == 1)
{
if (distancewall <= H / 2)
{
F_L = rho * pow(Gx, 2) * pow(a, 4) * C_Ly;
}
else
{
F_L = -rho * pow(Gx, 2) * pow(a, 4) * C_Ly;
}
}












t = Lookup_Thread(d, 8);



distancewall = 100; //i'm setting an generic initial value
begin_f_loop(f, t)
{
F_CENTROID(xw, f, t); //center of each face
NV_VV(vec, =, xw, -, P_POS(p)); //vector connecting the particle to the face center that is placed at the wall. The NV_VV do the operation between 2 vectors :(vec=xw-P_POS(p))





if (distancewall >= NV_MAG(vec))
{
distancewall = NV_MAG(vec);
}
}
end_f_loop(f, t)



cosine = vec[i] / (NV_MAG(vec)); //this is the cosine director because vec[i] is the i component of the vector vec




Thread* ttt = P_CELL_THREAD(p);




rho = C_R(c, ttt);
Ufx = C_U(c, ttt);
Ufy = C_V(c, ttt);
Ufz = C_W(c, ttt);
aUfx = fabs(Ufx); /* absolute value of Ufx */
aUfy = fabs(Ufy); /* absolute value of Ufy */
aUfz = fabs(Ufz); /* absolute value of Ufz */


/* local Reynolds number calculation*/

Renx = (rho * aUfx * Dh) / mu;
Reny = (rho * aUfy * Dh) / mu;
Renz = (rho * aUfz * Dh) / mu;


Gx = (Ufx) / H; /*max shear rate of the channel*/

a = P_DIAM(p); /* particle diameter */

z = fabs(((W / 2) - distancewall) / (H / 2));



if (z >= 0 && z <= 1.2)
{
C_Lz = ((0.047 - y * (0.077 / 0.6)) / 1.2) * z;
}
else if (z > 1.2)
{

C_Lz = (7.3046 * ((y * (0.077 / 0.6)) - 0.127) * pow(z, 2) - 17.531 * ((y * (0.077 / 0.6)) - 0.127) * z + 9.51863 * ((y * (0.077 / 0.6)) - 0.135405)); /*value of C_L for 2z/H */
}





if (i == 2)
{
if (distancewall <= W / 2)
{
F_L = (rho * pow(Gx, 2) * pow(a, 4) * C_Lz)*cosine;
}
else
{
F_L = -(rho * pow(Gx, 2) * pow(a, 4) * C_Lz)*cosine;
}
}


if (i == 0)
{
if (distancewall <= W / 2)
{
F_L = (rho * pow(Gx, 2) * pow(a, 4) * C_Lz)*cosine;
}
else
{
F_L = -(rho * pow(Gx, 2) * pow(a, 4) * C_Lz)*cosine;
}
}



Message("x position:%g\n", P_POS(p)[0]);
Message("y position:%g\n", P_POS(p)[1]);
Message("z position:%g\n", P_POS(p)[2]);


return (F_L / P_MASS(p));




}

__________