|
[Sponsors] |
![]() |
![]() |
#1 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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! |
|
![]() |
![]() |
![]() |
![]() |
#2 |
Senior Member
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34 ![]() ![]() |
example from Ansys Fluent Customization manual
Code:
/* 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)); } 0 - x axis 1 - y axis 2 - z axis
__________________
best regards ****************************** press LIKE if this message was helpful |
|
![]() |
![]() |
![]() |
![]() |
#3 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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! |
|
![]() |
![]() |
![]() |
![]() |
#4 |
Senior Member
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34 ![]() ![]() |
in example all three directions are considered.
use same approach
__________________
best regards ****************************** press LIKE if this message was helpful Last edited by AlexanderZ; November 3, 2020 at 04:13. |
|
![]() |
![]() |
![]() |
![]() |
#5 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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 |
|
![]() |
![]() |
![]() |
![]() |
#6 |
Senior Member
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34 ![]() ![]() |
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
__________________
best regards ****************************** press LIKE if this message was helpful |
|
![]() |
![]() |
![]() |
![]() |
#7 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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); } |
|
![]() |
![]() |
![]() |
![]() |
#8 |
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 ![]() |
return(force[i]/m)
|
|
![]() |
![]() |
![]() |
![]() |
#9 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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! |
|
![]() |
![]() |
![]() |
![]() |
#10 |
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 ![]() |
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. Code:
#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)); } |
|
![]() |
![]() |
![]() |
![]() |
#11 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
Thank you very much.
I will create an integration over the surface to get the force. |
|
![]() |
![]() |
![]() |
![]() |
#12 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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)); } |
|
![]() |
![]() |
![]() |
![]() |
#13 |
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#14 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#15 |
Senior Member
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26 ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#16 |
Member
L.A.Isanka
Join Date: Jul 2020
Posts: 55
Rep Power: 6 ![]() |
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)); } |
|
![]() |
![]() |
![]() |
![]() |
#17 | |
Member
Anshuman Sinha
Join Date: Oct 2018
Posts: 70
Rep Power: 7 ![]() |
Quote:
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 |
||
![]() |
![]() |
![]() |
![]() |
#18 |
New Member
Jacky Yang
Join Date: Jun 2023
Posts: 1
Rep Power: 0 ![]() |
/* 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. |
|
![]() |
![]() |
![]() |
Tags |
dpm bodyforce, interactions, vector |
|
|
![]() |
||||
Thread | Thread Starter | Forum | Replies | Last Post |
No matching function error: Phase change source term added to interMixingFoam | wavefunction | OpenFOAM Programming & Development | 2 | February 4, 2022 07:46 |
Creating a new field from terms of the turbulence model | HaZe | OpenFOAM Programming & Development | 15 | November 24, 2014 13:51 |
DPM UDF particle position using the macro P_POS(p)[i] | dm2747 | FLUENT | 0 | April 17, 2009 01:29 |
Missing math.h header | Travis | FLUENT | 4 | January 15, 2009 11:48 |
URGENT - DPM Macro for custom bc | cfd-user | FLUENT | 2 | July 21, 2008 11:26 |