CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

UDF for statistics of dpm trapped on the wall

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By obscureed
  • 1 Post By obscureed

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 2, 2018, 22:13
Post UDF for statistics of dpm trapped on the wall
  #1
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
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!
bigfish_1023 is offline   Reply With Quote

Old   January 3, 2018, 06:41
Default
  #2
Senior Member
 
Join Date: Sep 2017
Posts: 224
Rep Power: 8
obscureed is on a distinguished road
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.
wc34071209 likes this.
obscureed is offline   Reply With Quote

Old   January 3, 2018, 08:50
Default
  #3
Senior Member
 
Join Date: Nov 2013
Posts: 1,860
Rep Power: 23
pakk will become famous soon enough
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.
pakk is offline   Reply With Quote

Old   January 3, 2018, 20:35
Default
  #4
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
Quote:
Originally Posted by obscureed View Post
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.
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);
}
}
bigfish_1023 is offline   Reply With Quote

Old   January 3, 2018, 21:46
Default
  #5
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
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?
bigfish_1023 is offline   Reply With Quote

Old   January 4, 2018, 06:48
Default
  #6
Senior Member
 
Join Date: Sep 2017
Posts: 224
Rep Power: 8
obscureed is on a distinguished road
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:
Code:
if (rp_axi_swirl)
  Message("Trapped %12.8g,%12.8g,%12.8g\n", P_POS(p)[0], R, P_FLOW_RATE(p));
obscureed is offline   Reply With Quote

Old   January 4, 2018, 07:08
Default
  #7
Senior Member
 
Join Date: Sep 2017
Posts: 224
Rep Power: 8
obscureed is on a distinguished road
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:
Code:
/* In the declarations (always at the start of the function): */
  real area,rate;
  real NV_VEC(area_vec);
  Thread *ct;
  cell_t c;
/* ... */
#define UDM_FOR_WALL_COLLISIONS 0
#define NOTEWORTHY_NORMAL_PVEL 1.0
    if(vn >= NOTEWORTHY_NORMAL_PVEL)
    {
      if (rp_axi_swirl)
        Message("%12.8g,%12.8g,%12.8g \n", P_POS(p)[0], R, P_FLOW_RATE(p));

      F_AREA(area_vec, f, ft);
      area = NV_MAG(area_vec);

      rate = P_FLOW_RATE(p) / area; /* reconsider use of P_FLOW_RATE if particles are unsteady */
      F_UDMI(f, ft, UDM_FOR_WALL_COLLISIONS) += rate;
      /* Copy value to neighbouring cell: */
      c  = F_C0(f, ft);
      ct = THREAD_T0(ft);
      if (NNULLP(ct))
        C_UDMI(c, ct, UDM_FOR_WALL_COLLISIONS) += rate;
      c  = F_C1(f, ft);
      ct = THREAD_T1(ft);
      if (NNULLP(ct))
        C_UDMI(c, ct, UDM_FOR_WALL_COLLISIONS) += rate;

      return PATH_ABORT;
    }
wc34071209 likes this.
obscureed is offline   Reply With Quote

Old   January 7, 2018, 19:27
Default
  #8
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
Quote:
Originally Posted by obscureed View Post
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:
Code:
if (rp_axi_swirl)
  Message("Trapped %12.8g,%12.8g,%12.8g\n", P_POS(p)[0], R, P_FLOW_RATE(p));
Hi, obscureed. Thanks again! It works and displays what I want.
bigfish_1023 is offline   Reply With Quote

Old   January 9, 2018, 20:55
Default
  #9
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
Quote:
Originally Posted by obscureed View Post
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:
Code:
/* In the declarations (always at the start of the function): */
  real area,rate;
  real NV_VEC(area_vec);
  Thread *ct;
  cell_t c;
/* ... */
#define UDM_FOR_WALL_COLLISIONS 0
#define NOTEWORTHY_NORMAL_PVEL 1.0
    if(vn >= NOTEWORTHY_NORMAL_PVEL)
    {
      if (rp_axi_swirl)
        Message("%12.8g,%12.8g,%12.8g \n", P_POS(p)[0], R, P_FLOW_RATE(p));

      F_AREA(area_vec, f, ft);
      area = NV_MAG(area_vec);

      rate = P_FLOW_RATE(p) / area; /* reconsider use of P_FLOW_RATE if particles are unsteady */
      F_UDMI(f, ft, UDM_FOR_WALL_COLLISIONS) += rate;
      /* Copy value to neighbouring cell: */
      c  = F_C0(f, ft);
      ct = THREAD_T0(ft);
      if (NNULLP(ct))
        C_UDMI(c, ct, UDM_FOR_WALL_COLLISIONS) += rate;
      c  = F_C1(f, ft);
      ct = THREAD_T1(ft);
      if (NNULLP(ct))
        C_UDMI(c, ct, UDM_FOR_WALL_COLLISIONS) += rate;

      return PATH_ABORT;
    }
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.
bigfish_1023 is offline   Reply With Quote

Old   January 11, 2018, 08:54
Default
  #10
Senior Member
 
Join Date: Sep 2017
Posts: 224
Rep Power: 8
obscureed is on a distinguished road
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 is offline   Reply With Quote

Old   January 15, 2018, 20:48
Default
  #11
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
Quote:
Originally Posted by obscureed View Post
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!
bigfish_1023 is offline   Reply With Quote

Old   January 16, 2018, 09:45
Default
  #12
Senior Member
 
Join Date: Nov 2013
Posts: 1,860
Rep Power: 23
pakk will become famous soon enough
Quote:
Originally Posted by bigfish_1023 View Post
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.
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).
pakk is offline   Reply With Quote

Old   January 16, 2018, 22:23
Default
  #13
New Member
 
BigZ
Join Date: Jan 2018
Posts: 11
Rep Power: 5
bigfish_1023 is on a distinguished road
Quote:
Originally Posted by pakk View Post
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?
bigfish_1023 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
UDF for suction velocity on a wall ahmed425 Fluent UDF and Scheme Programming 0 November 2, 2017 09:50
Conditional Release of DPM particles - UDF nvschandra Fluent UDF and Scheme Programming 0 December 10, 2013 11:02
DPM UDF trapped particle deposition barzin Fluent UDF and Scheme Programming 0 September 4, 2012 10:27
DPM UDF particle position using the macro P_POS(p)[i] dm2747 FLUENT 0 April 17, 2009 01:29
UDF to monitor minimum pressure on the wall Arvind Jayaprakash FLUENT 0 September 20, 2007 11:53


All times are GMT -4. The time now is 15:09.