[Achini]

Hello!

I'm using DPM body force macro to model particle interactions.
According to the UDF manual, this macro should return the real value of the acceleration due to the body force (in m/s 2).

How to include the direction of this acceleration?
It gives an error when returning the acceleration as a vector.

Thank you!

[AlexanderZ]

example from Ansys Fluent Customization manual

/* UDF for computing the magnetic force on a charged particle */
#include "udf.h"
#define Q 1.0 /* particle electric charge */
#define BZ 3.0 /* z component of magnetic field */
#define TSTART 18.0 /* field applied at t = tstart */
/* Calculate magnetic force on charged particle. Magnetic */
/* force is particle charge times cross product of particle */
/* velocity with magnetic field: Fx= q*bz*Vy, Fy= -q*bz*Vx */
DEFINE_DPM_BODY_FORCE(particle_body_force,p,i)
{
real bforce=0;
if(P_TIME(p)>=TSTART)
{
if(i==0) bforce=Q*BZ*P_VEL(p)[1];
else if(i==1) bforce=-Q*BZ*P_VEL(p)[0];
}
else
bforce=0.0;
/* an acceleration should be returned */
return (bforce/P_MASS(p));
}

Where i is An index (0, 1, or 2) that identifies the Cartesian component of the body force that is to be returned by the function.
0 - x axis
1 - y axis
2 - z axis

[Achini]

Hello AlexanderZ,

Thank you for your reply.
In my application, I have a three-dimensional force and I want to return all three components (0,1,2).

Do you know how to do that?
I'd be grateful if you could help me with this.

Thank you!

[AlexanderZ]

in example all three directions are considered.
use same approach

[Achini]

I used return (bforce/m);

I have declared bforce as a vector and m is the particle mass.
It gives the following error. Could you please tell me the reason for this error?

..\..\src\newudf.c(87) : error C2296: '/' : illegal, left operand has type 'real [3]'
..\..\src\newudf.c(87) : warning C4033: 'particle_body_force' must return a value

[AlexanderZ]

i have no idea about the code you are trying to use. problem is at line 87

1. compile your code
2. show log output
3. show code

[Achini]

The code I mentioned below is to model particle wall interactions. I need to get the force exerted by the wall on each particle.

Also, the force only acts on half of the particles in the domain. I have used an integer "a" and took its modulus with a while loop to make the force apply only on half of the particles.

And for each particle loop, I have created a loop over the boundary face in order to get the force exerted by each boundary cell.

This is the error I get when compiling.

..\..\src\newudf.c(88) : error C2296: '/' : illegal, left operand has type 'real [3]'
..\..\src\newudf.c(88) : warning C4033: 'particle_body_force' must return a value
..\..\src\newudf.c(89) : fatal error C1075: end of file found before the left brace '{' at '..\..\src\newudf.c(28)' was matched

Please help me fix errors with this code.
Thank you.

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

#define K 0.000003

double bf,a,remainder,dw,x1,y1,z1,x0,y0,z0,fmag,m,acc;
real force[ND_ND];
real xw[ND_ND];
real vec[ND_ND];
real n[ND_ND];

div_t divide;


DEFINE_DPM_BODY_FORCE(particle_body_force,p,i) /* Fluent macro to define particle body force*/
{
Domain *d;
cell_t c;
Thread *t;
/*Particle *p;*/
Particle *pi;
face_t f;

a = 1;

begin_particle_cell_loop(pi, c, t)
{
bf = 0;

c = P_CELL(pi); /* Get the cell and thread that the particle is currently in */
t = P_CELL_THREAD(p);
m = P_MASS(pi);

x0 = P_POS(pi)[0];
y0 = P_POS(pi)[1];
z0 = P_POS(pi)[2];


begin_f_loop(f, t)
{

F_CENTROID(xw,f,t);
x1 = xw[0];
y1 = xw[1];
z1 = xw[2];

/* dw = sqrt(SQR(x1-x0) + SQR(y1-y0) + SQR(z1-z0)); wall to particle distance calculation*/

divide = div(a,2);

remainder = divide.rem;

while (remainder == 1)
{

if (THREAD_TYPE(t)==THREAD_F_WALL)
{

vec[0] = x1-x0; /* vector in the direction of the force */
vec[1] = y1-y0;
vec[2] = z1-z0;

dw = NV_MAG(vec); /* Magnitude of the vector in the direction of the force*/

bf = K/(dw*dw*dw*dw*dw*dw*dw); /* calculation of the resulting force */


n[0]=vec[0]/dw; /* unit vector in the diretion of the force*/
n[1]=vec[1]/dw;
n[2]=vec[2]/dw;


NV_V_VS(force, =, force, +, force, *, bf); /* calculates the sum of all interaction forces*/

}
else

bf=0;
}

end_f_loop(f, t)
}
a = a+1;
}

return (force/m);
}

[pakk]

return(force[i]/m)

[Achini]

return(force[i]/m)

This worked. Thank you!

but it gives the following error even though I have added brackets.

..\..\src\newudf.c(89) : fatal error C1075: end of file found before the left brace '{' at '..\..\src\newudf.c(28)' was matched

Below mentioned is my code. I want to get the body force on each particle. therefore, should I put the return bforce; command inside the particle loop?


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

#define K 0.000003

double bf,a,remainder,dw,x1,y1,z1,x0,y0,z0,fmag,m,acc;
real force[ND_ND];
real xw[ND_ND];
real vec[ND_ND];
real n[ND_ND];

div_t divide;


DEFINE_DPM_BODY_FORCE(particle_body_force,p,i) /* Fluent macro to define particle body force*/
{
Domain *d;
cell_t c;
Thread *t;
/*Particle *p;*/
Particle *pi;
face_t f;

a = 1;

begin_particle_cell_loop(pi, c, t)
{
bf = 0;

c = P_CELL(pi); /* Get the cell and thread that the particle is currently in */
t = P_CELL_THREAD(p);
m = P_MASS(pi);

x0 = P_POS(pi)[0];
y0 = P_POS(pi)[1];
z0 = P_POS(pi)[2];


begin_f_loop(f, t)
{

F_CENTROID(xw,f,t);
x1 = xw[0];
y1 = xw[1];
z1 = xw[2];

/* dw = sqrt(SQR(x1-x0) + SQR(y1-y0) + SQR(z1-z0)); wall to particle distance calculation*/

divide = div(a,2);

remainder = divide.rem;

while (remainder == 1)
{

if (THREAD_TYPE(t)==THREAD_F_WALL)
{

vec[0] = x1-x0; /* vector in the direction of the force */
vec[1] = y1-y0;
vec[2] = z1-z0;

dw = NV_MAG(vec); /* Magnitude of the vector in the direction of the force*/

bf = K/(dw*dw*dw*dw*dw*dw*dw); /* calculation of the resulting force */


n[0]=vec[0]/dw; /* unit vector in the diretion of the force*/
n[1]=vec[1]/dw;
n[2]=vec[2]/dw;


NV_V_VS(force, =, force, +, force, *, bf); /* calculates the sum of all interaction forces*/

}
else

bf=0;
}
return (force[i]/m);
end_f_loop(f, t)
}
a = a+1;
}

/*return (force[i]/m);*/
}


Thank you!

[pakk]

Well, you are doing a completely useless cell-loop... You are asked to give the acceleration on particle p, there is no need to loop through all particles anymore...

I cleaned up your code, removed most of the unneccessary lines, simplified things, and made comments where you have to change things but I don't know what because I don't know which problem you are solving.

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

#define K 0.000003

DEFINE_DPM_BODY_FORCE(particle_body_force,p,i) /* Fluent macro to define particle body force*/
{
 double bf,dw2;
 real xw[ND_ND];
 real force[ND_ND];
 real vec[ND_ND];
 
 cell_t c;
 Thread *t;
 face_t f;
 
 int 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 
   ** YOU DON'T WANT THIS: YOU WANT THE THREAD OF THE WALL THAT YOU ARE INTERACTING WITH**
   the particle most likely is in a gas thread, so you will never see interaction in this way. */
 t = P_CELL_THREAD(p);
 
 NV_S(force, =, 0.0);
 begin_f_loop(f, t)
 {  
  if (THREAD_TYPE(t)==THREAD_F_WALL)
  {
   F_CENTROID(xw,f,t);
   NV_VV(vec, =, xw, -, P_POS(p));/* vector in the direction of the force */
   dw2 = NV_MAG2(vec); /* Magnitude of the vector in the direction of the force*/
   bf = K/(dw2*dw2*dw2*dw2); /* calculation of the resulting force divided by distance*/
   NV_V_VS(force, =, force, +, vec,*, bf); /* calculates the sum of all interaction forces*/
   /* As I told you in another question: the line above makes no sense physically. Currently, your 
      result depends on the mesh size. If you have a smaller mesh in your wall, you will get a
      bigger force. This makes no sense. I have no idea which physics you want to simulate, so
      I don't know how to fix this. My two guesses were that you want to calculate the
      maximum force (or the minimum distance), or that you want to integrate over the surface,
      but in the latter case the face size should be taken into account. */ 
   end_f_loop(f, t) 
  }
 }
 return (force[i]/P_MASS(p));
}

I did not test this code, so typos can still be there.

[Achini]

Thank you very much.

I will create an integration over the surface to get the force.

[Achini]

Hi Pakk,

I edited the UDF as below.
I took the sum of force multiplied with cell face area, calculated its average to get the net force.
It compiles without any error.

But I could see particles getting attracted to the wall.
Will you be able to help me figure out why?
I am much grateful for your help.

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

#define K 0.000003

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

cell_t c;
Thread *t;
face_t f;

int 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); */

/*t = P_CELL_THREAD(p);*/

NV_S(force, =, 0.0);
area=0;
begin_f_loop(f, t)
{
if (THREAD_TYPE(t)==THREAD_F_WALL)
{
F_CENTROID(xw,f,t);
F_AREA(A,f,t);
facearea = sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); /* to cell the cell area*/
area = area+facearea; /* addition of all cell feca area to get the total area*/

NV_VV(vec, =, xw, -, P_POS(p));/* vector in the direction of the force */

dw2 = NV_MAG2(vec); /* Magnitude of the vector in the direction of the force*/
bf = K/(dw2*dw2*dw2*dw2); /* calculation of the resulting force divided by distance*/
e=bf*facearea;/* force on face area*/
NV_V_VS(force, =, force, +, vec,*, e); /* calculates the sum of all interaction forces*/


end_f_loop(f, t)
}

bforce[i]=force[i]/area;
}
return (force[i]/P_MASS(p));
}

[pakk]

Well, that's what happens when you add a force that attracts particles to the wall...

If you want the force in the other direction, you should add a minus sign.

[Achini]

Hi Pakk,

I'm getting an error running this. DO you know the reason behind this?

Error: received a fatal signal (Segmentation fault).
Error Object: #f

Also, how should I add the minus sign? it also gives errors.

[pakk]

You forgot to assign t to a thread.

You want your particle to interact with a wall. Tell fluent which wall you are taking about! GET_THREAD should do the trick, if I remember correctly.

[Achini]

Thank you Pakk for your support.
I used "t = Lookup_Thread(d, 5);" to assign the thread with boundary wall id.
Please have a look at my UDF and comment if there's anything wrong.
Thanks again!

#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,dw2,area,facearea,e,K;
real xw[ND_ND];
real force[ND_ND];
real bforce[ND_ND];
real vec[ND_ND];
real A[ND_ND];

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

d = Get_Domain(1);

K = 3.3*(10^(-34));


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);*/

t = Lookup_Thread(d, 5);

NV_S(force, =, 0.0); /* reset force*/
area=0;

begin_f_loop(f, t)
{

/*if (THREAD_TYPE(t)==THREAD_F_WALL)*/

F_CENTROID(xw,f,t);
F_AREA(A,f,t);

facearea = sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); /* to cell the cell area*/
area = area+facearea; /* addition of all cell feca area to get the total area*/

NV_VV(vec, =, xw, -, P_POS(p)); /* vector in the direction of the force */

dw2 = NV_MAG2(vec); /* Magnitude of the vector in the direction of the force*/

bf = K/(dw2*dw2*dw2*dw2*dw2*dw2*dw2); /* calculation of the resulting force divided by distance*/

e=bf*facearea; /* force on cell face area*/

NV_V_VS(force, =, force, +, vec,*, e); /* calculates the sum of all interaction forces*/


end_f_loop(f, t)

}

bforce[i]=(-(force[i])/area);

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

[Anshs]

Originally Posted by Achini

Thank you Pakk for your support.
I used "t = Lookup_Thread(d, 5);" to assign the thread with boundary wall id.

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

}

Hi Achini!

How to compute body force between any 2 DPM particles at any instance. Let's say I have 10 DPM partciles and I want to assign a attraction force between all these particles. Thanks

[Jacky Yang]

/* dw = sqrt(SQR(x1-x0) + SQR(y1-y0) + SQR(z1-z0)); wall to particle distance calculation*/

Genius bro! I was inspired by your wall distance calculating method.