UDF for statistics of dpm trapped on the wall
 

Hi, everyone. I wanted to simulate the particles which are trapped on the wall when it exceeds the critical velocity. I tried the macro of DEFINE_DPM_BC and it worked. Additionally I need to give the statistical results about the number of particles trapped at different position of the wall. Thus in the above DEFINE_DPM_BC , I created a *.dat file using fprintf and write the relevant variable into this file, but it didn’t work and reported error.

So could anyone can help me or to explain why?
Thanks!

There are several things that could be going wrong -- maybe you should post your code? First, if you are running in parallel, then writing to a single file from multiple nodes is dangerous/impossible. A simpler plan might be to write to the Fluent text window (using "Message(...)" instead of "fprintf(fp,...)") -- you can save this to a text file as a transcript.

A nicer plan might be to allocate a user-defined memory and store location-based results on the wall facets. The UDFs for this might be more lengthy, but they should not require too much extra thought to run in parallel. I try to store mesh-independent results, which in this case means "per unit face area". You can obtain total rates from this by surface integral reports.

Best regards, Ed.

I agree with obscureed. You probably made a mistake in your code, but without seeing your code it is very difficult to guess what.

And his "better plan" is exactly what I have done in the past when I wanted to do something similar.

Thanks for your quick reply! Below is my code. As you guess, i was running this 2D simulation in parallel, and some error occured. But it didnt work in serial mode either. I will try the Message macro and could you give me more details about the second way?

#include "udf.h"

DEFINE_DPM_BC(boundary_capture,p,t,f,f_normal,dim)
{
real alpha; /*angle of particle hitting the wall*/
real nor_coeff = 0.3;
real tan_coeff = 0.9;
real vn = 0;
real normal[3];
cell_t c0;
Thread *t0;
int i, idim = dim;
real NV_VEC(x);

FILE *fd;
fd = fopen("result.dat", "a");


#if RP_2D
/* dim is always 2 in 2D compilation. Need special treatment for 2d
axisymmetric and swirl flows */
if (rp_axi_swirl)
{
real R = sqrt(P_POS(p)[1]*P_POS(p)[1] +
P_POS(p)[2]*P_POS(p)[2]);
if (R > 1.e-20)
{
idim = 3;
normal[0] = f_normal[0];
normal[1] = (f_normal[1]*P_POS(p)[1])/R;
normal[2] = (f_normal[1]*P_POS(p)[2])/R;
}
else
{
for (i=0; i<idim; i++)
normal[i] = f_normal[i];
}
}
else
#endif
for (i=0; i<idim; i++)
normal[i] = f_normal[i];
if(p->type==DPM_TYPE_INERT)
{
alpha = M_PI/2. - acos(MAX(-1.,MIN(1.,NV_DOT(normal,P_VEL(p))/MAX(NV_MAG(P_VEL(p)),DPM_SMALL))));
if ((NNULLP(t)) && (THREAD_TYPE(t) == THREAD_F_WALL))
F_CENTROID(x,f,t);
for(i=0; i<idim; i++)
vn += P_VEL(p)[i]*normal[i];
if(vn >= 1.0)
{
fprintf(fd, "%f,%f,%f \n", P_POS(p)[0], P_POS(p)[1], P_FLOW_RATE(p));
return PATH_ABORT;
}
else
{
for(i=0; i<idim; i++)
P_VEL(p)[i] -= vn*normal[i];
/* Apply tangential coefficient of restitution. */
for(i=0; i<idim; i++)
P_VEL(p)[i] *= tan_coeff;
/* Add reflected normal velocity. */
for(i=0; i<idim; i++)
P_VEL(p)[i] -= nor_coeff*vn*normal[i];

/* Store new velocity in P_VEL0 of particle */
for(i=0; i<idim; i++)
P_VEL0(p)[i] = P_VEL(p)[i];
return PATH_ACTIVE;
}
fclose(fd);
}
}

Hi, obscureed. I have tried the Message macro, and it succeeds. Thank you very much.
But there are something wrong happenning. I use Message to display the information of the trapped particles, such as the position. In my 2D axisymmetric simulation, the domain of x and y are all positive, but the coordinate of y of some trapped particle using P_POS(p)[1] is negative. why?

Hi bigfish_1023. Regarding your question about surprising P_POS(p)[1] -- this is not an error. In axisymmetric-swirl 2D simulations, DPM particles have 3D coordinates and velocities. The 2D elements such as walls and fluid velocities are rotated around the x-axis to recreate the axisymmetric surroundings for the particles. You can see this if you pick through the code wrapped in "if (rp_axi_swirl)" -- a 3D copy of the face's 2D normal vector is made. So you might want something like:

Oh, and the idea about saving facet-based data to User-Defined Memory: this is slightly less attractive for 2D simulations, because the walls are only lines and so contour plots on the walls don't look good. In 3D, you can build up records of where noteworthy collisions have occurred, and display the results as contour plots. The steps are roughly as follows:
(1) Request one User-Defined Memory (UDM).
(2) Set up an on-demand UDF to set the UDM to zero.
(3) Inside the DEFINE_DPM_BC, add particle-flowrate-per-unit-area to the UDM for each noteworthy collision.
(4) Display contours of UDM-0, or report totals by surface integrals.

The code in step (3) is something like this:

Hi, obscureed. Thanks again! It works and displays what I want.

Obscureed. sorry to bother you again. i have some other questions about the code you mentioned above. i have reported the totals by surface integrals of UDM-0， but the obtained value was not the same as the sum total of P_FLOW_RATE(p). In my opinion, these two ways should lead to the same result. Additionally, both the surface integrals of UDM-0 and the sum of P_FLOW_RATE(p) were much larger than the inlet mass flow of particle. this is unreasonable. I have checked the code carefully and still dont konw why.

On your observation that reported wall-collision particle mass flowrate is greater than released particle mass flowrate: some possibilities leap to mind:
(1) Perhaps you released the injection multiple times (for example, during two-way interaction, every N iterations) without resetting the UDM.
(2) Perhaps the particles do not die as intended when you have noted their effects. Hopefully you can display some particle tracks and actually see them not surviving. A good way to debug this would be to lower your criterion (NOTEWORTHY_NORMAL_PVEL) to zero, temporarily, so that all particles should die at their first wall collision. One way to kill a particle is "MARK_TP(tp,P_FL_REMOVED);" or possibly "MARK_PARTICLE(p,P_FL_REMOVED);" or both.
(2a) Really confusing situations can arise if you change the DPM numerics to Runge-Kutta. R-K integration involves making some trial steps before the final step. If one of these trial steps triggers the DEFINE_DPM_BC, and then the final step also triggers it, you might get some double-counting. I can't suggest a fix for this apart from avoiding R-K.
(3) Perhaps you have multiple stochastic tries in the Turbulent Dispersion settings for DPM. If you have 10 stochastic tries per particle, then they each get (1/10)th of the effect (in momentum sources, etc), but I think they carry the full P_FLOWRATE(p) or TP_FLOWRATE(tp). Another potentially confusing area is multiple size bins in a size distribution (Rosin-Rammler etc), but from memory I think that these carry only their fraction of the mass flowrate. The best way to be sure it to check.
(4) Oh, you're doing 2D axisymmetric swirl. There are sometimes some weird and fiddly factors of (2*M_PI) flying around. It's horrible. Make sure that the UDMs are reset to zero, kill your entire injection and see if the totals add up.

But you said the UDF flowrates are "much larger" than the released flowrates -- most of these explanations will only be a few times larger. Dividing by area twice by mistake? -- that would do it. Please tell us what you find. If you still can't see what's going on, you might want to print out the P_FLOWRATE of every particle that is released -- a new UDF, either DEFINE_DPM_INJECTION_INIT or DEFINE_SCALAR_UPDATE with the initialize flag. Failing that, please post your code and see who else has any ideas.

Obscureed. Thanks for your detailed explanation. Based on your suggestions, I have found the reason. In my simulation, the particles experienced reaction, and the P_FLOW_RATE(p) gives the initial mass flow rate of a stream, not the description that it returns the strength multiplied by P_MASS(p) at the current particle position. On the other hand, as you pointed out, I have multiple stochastic tries. Thus the deposition rate of trapped particles on the specific position is P_FLOW_RATE(p)*P_MASS(p)/P_INIT_MASS(p)/(number of tries). And the totals by surface integrals of UDM-0 is just nearly (2*M_PI) times larger than the sum of deposition rate calculated by the above equation. Please tell me the reason.
Moreover, the usage of on demand to reset UDMs to zero is not preferred as the DPM iterations are going on. Is there another way to reset the UDMs to zero automatically after the complement of each DPM iteration.

Thanks again!

The reason is that Fluent in axisymmetric mode calculates everything per radian, whereas in the surface integral mode it is calculated for the complete geometry (=2*M_PI radians).

Hi， pakk. Thank you very much! Does that mean the calculation of face area from NV_MAG(area_vec) in axisymmetric mode should be multiplied by 2*M_PI?