[sjbub]

I am modeling some inert particles flow through a pipe and want to make a udf to define sticking criteria for the particles. The plan is have the ufd check the particle diamter and velcity at impact, calculate the capture velocity and then abort for sticking or track for the particle to bounce off the wall and continue on.

Anyyyyways my question is that the user manual says that if you are going to use this DEFINE_DPM_BC function you must compute the new velocity of the particle after it hits the wall within your udf. I am not altering the coefficients of restitution (my setting is 0.99 - very small particles and assuming no deformation on impact) so do I still need to include this calculation in my udf?

Sincerely-
trying-not-to-pull-my-hair-out

[`e`]

In your case, the particles have two options (use an if statement) once they hit the wall:
(1) particle is captured by the wall ("stick" == trapped; using a return of PATH_ABORT)
(2) particle continues (using a return of PATH_ACTIVE) with reflected velocities. Either assume perfectly elastic conditions or use a relevant equation to calculate the new velocities.

If you simply return a particle from the boundary with the same velocities it had, then I suspect the particle would collide with the same boundary again, creating an infinite loop (try this!).

[sjbub]

Originally Posted by `e`

In your case, the particles have two options (use an if statement) once they hit the wall:
(1) particle is captured by the wall ("stick" == trapped; using a return of PATH_ABORT)
(2) particle continues (using a return of PATH_ACTIVE) with reflected velocities. Either assume perfectly elastic conditions or use a relevant equation to calculate the new velocities.

If you simply return a particle from the boundary with the same velocities it had, then I suspect the particle would collide with the same boundary again, creating an infinite loop (try this!).

thanks e for the kind and quick response

I have decided I lack the programming experience to get this udf together properly so I will not be able to let you know how it goes. But hopefully soon I will have a not have a not-as-sophisticated method but a workable one

[`e`]

The concept for this UDF isn't too complicated but I can understand your hesitation without programming experience, but surely it'll be fun!

Here's a start anyway, you'll need to add in the capture velocity equation and how you'll determine the reflected particle velocities.

#include "udf.h"

DEFINE_DPM_BC(my_dpm_bc,p,t,f,f_normal,dim)
{
	/* Variable declarations */
	real u, v, w, diameter, capture_velocity, critical_capture_velocity=0.;

	/* Variable evaluations */
	u=P_VEL(p)[0];
	v=P_VEL(p)[1];
	w=P_VEL(p)[2];
	diameter=P_DIAM(p);

	// Calculate the capture velocity with u, v, w and diameter here, then:
	if (capture_velocity > critical_capture_velocity) // particle travels fast
	{
		// Update particle velocities for reflection case
		return PATH_ACTIVE;
	}
	else // particle travels slow and is stuck
	{
		return PATH_ABORT;
	}
}

[sjbub]

You are too kind and thank you for the encouragement.

So here is my first attempt. I don't need to calculate the particle velocity so I deleted the u, v and w velocity components. I grabbed an example from the user manual to calc the reflection velocity.

I am a bit confused on why in the variable declarations the critical_capture_velocity is equal to zero? (I did try to search this). Is this initializing the variable maybe?

I think I'll be up for a while working on this. Let you know how it goes.

PHP Code:

/***********************************************************************
UDF Particle Sticking Criteria
If a particle impacts the wall and its velocity is less than the 
capture velocity for its diameter it will stick. if not it will reflect.
************************************************************************/

#include "udf.h"

DEFINE_DPM_BC(my_dpm_bc,p,t,f,f_normal,dim)
{
    /* Variable declarations */
    real u, v, w, diameter, critical_capture_velocity=0.;
    real alpha; /* angle of particle path with face normal */
    real vn=0.;
       real nor_coeff = 0.99.;
       real tan_coeff = 0.99;
       real normal[3];
       int i, idim = dim;
    real NV_VEC(x);

    /* Variable evaluations */
    diameter=P_DIAM(p);

    // Calculate the capture velocity with u, v, w and diameter here, then:
       critical_capture_velocity = (5.15e-07)/diamter;
    if (P_VEL(p) > critical_capture_velocity) // particle travels fast
    {
        // Update particle velocities for reflection case (grabbed from fluent example and updated coeff)

#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);

         /* calculate the normal component, rescale its magnitude by
           the coefficient of restitution and subtract the change */

         /* Compute normal velocity. */
         for(i=0; i<idim; i++)
           vn += P_VEL(p)[i]*normal[i];

         /* Subtract off normal velocity. */
           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;
    }
    else // particle travels slow and is stuck
    {
        return PATH_ABORT;
    } 


[`e`]

No worries, syntax and rules of coding languages can be easily found on the internet but it's a bit trickier starting a UDF from scratch without the big picture.

You don't need to initialise the 'critical_capture_velocity' variable; you could either hard code a constant in this initialisation (instead of having zero) or calculate a value as you have now done.

Which velocity are you intending to compare with the critical capture velocity, a single component (perpendicular, etc) or magnitude? You're currently comparing an array (P_VEL(p)) with a scalar (critical_capture_velocity).

Also: 'diamter' -> 'diameter'

[sjbub]

I would like to compare the velocity magnitude of the particle to the critical capture velocity.

I thought that P_VEL(p) was the velocity magnitude of the particle but I'm guessing if its an array it holds the velocity info for all the particles.

Let the google searching commence!

[sjbub]

Ok I am seeing why you had u, v and w originally in the variable evaluations. I need to calculate the particles velocity magnitude.

(slowly but surely)

[sjbub]

I wanted to say thanks again for giving me a place to start from. I actually do enjoy doing this sort of thing its just such a jump from beginning to anything else. A long time ago I took a C class but it was kind of a disacter since most all except myself and another person had programming experience so it was taught to that intermediate level.

Anyhoo since you were so nice I wanted to let you know I am working through getting all the errors I am getting when compiling in fluent. Hopefully soon I'll have this firgured out.

PHP Code:

/***********************************************************************
UDF Particle Sticking Criteria
If a particle impacts the wall and its velocity is less than the 
capture velocity for its diameter it will stick. if not it will reflect
************************************************************************
#include "udf.h"
DEFINE_DPM_BC(my_dpm_bc,p,t,f,f_normal,dim)
{
/* Variable declarations */
real u, v, w, diameter, particle_velocity, critical_capture_velocity;
/* Variable evaluations */
diameter=P_DIAM(p);
u=P_VEL(p)[0];
v=P_VEL(p)[1];
w=P_VEL(p)[2];
/* Calculate the particle velocity with u, v, w and critical capture velocity with the diameter here */ 
particle_velocity = sqrt (pow (u,2) + pow (v,2) + pow (w,2));
critical_capture_velocity = (0.000000515)/diameter;
/* particle travels fast */
if (particle_velocity > critical_capture_velocity)
/* Update particle velocities for reflection case (grabbed from fluent example and updated coeff */
{
#include "udf.h"
DEFINE_DPM_BC(bc_reflect,p,t,f,f_normal,dim)
{
real alpha; /* angle of particle path with face normal */
real vn=0.;
real nor_coeff = 1.;
real tan_coeff = 0.3;
real normal[3];
int i, idim = dim;
real NV_VEC(x);
#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);
/* Calculate the normal component, rescale its magnitude by
the coefficient of restitution and subtract the change */
/* Compute normal velocity. */
for(i=0; i<idim; i++)
vn += P_VEL(p)[i]*normal[i];
/* Subtract off normal velocity. */
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;

}
}
}
/* particle travels slow and is stuck*/
else
return PATH_ABORT;
} 


[sjbub]

The latest excitement...

..\..\src\my_dpm_bc.c(12) : error C2061: syntax error : identifier 'u'
..\..\src\my_dpm_bc.c(12) : error C2059: syntax error : ';'
..\..\src\my_dpm_bc.c(13) : error C2061: syntax error : identifier 'v'
..\..\src\my_dpm_bc.c(13) : error C2059: syntax error : ';'
..\..\src\my_dpm_bc.c(14) : error C2061: syntax error : identifier 'w'
..\..\src\my_dpm_bc.c(14) : error C2059: syntax error : ';'
..\..\src\my_dpm_bc.c(15) : error C2061: syntax error : identifier 'diameter'
..\..\src\my_dpm_bc.c(15) : error C2059: syntax error : ';'
..\..\src\my_dpm_bc.c(16) : error C2061: syntax error : identifier 'particle_velocity'
..\..\src\my_dpm_bc.c(16) : error C2059: syntax error : ';'
..\..\src\my_dpm_bc.c(17) : error C2061: syntax error : identifier 'critical_capture_velocity'
..\..\src\my_dpm_bc.c(17) : error C2059: syntax error : ';'
..\..\src\my_dpm_bc.c(20) : error C2065: 'p' : undeclared identifier
..\..\src\my_dpm_bc.c(20) : error C2099: initializer is not a constant
..\..\src\my_dpm_bc.c(21) : error C2109: subscript requires array or pointer type
..\..\src\my_dpm_bc.c(22) : error C2109: subscript requires array or pointer type
..\..\src\my_dpm_bc.c(23) : error C2109: subscript requires array or pointer type
..\..\src\my_dpm_bc.c(27) : error C2099: initializer is not a constant
..\..\src\my_dpm_bc.c(28) : error C2099: initializer is not a constant
..\..\src\my_dpm_bc.c(32) : error C2059: syntax error : 'if'

PHP Code:

/***********************************************************************
UDF Particle Sticking Criteria
If a particle impacts the wall and its velocity is less than the 
capture velocity for its diameter it will stick. if not it will reflect
************************************************************************
#include "udf.h"
DEFINE_DPM_BC(my_dpm_bc,p,t,f,f_normal,dim)
{
/* Variable declarations */
real u;
real v;
real w;
real diameter;
real particle_velocity; 
real critical_capture_velocity;
/* Variable evaluations */
diameter=P_DIAM(p);
u=P_VEL(p)[0];
v=P_VEL(p)[1];
w=P_VEL(p)[2];
/* Calculate the particle velocity with u, v, w and critical capture velocity with the diameter here */ 
particle_velocity = sqrt (pow (u,2) + pow (v,2) + pow (w,2));
critical_capture_velocity = (0.000000515)/diameter;
/* particle travels fast */
if (particle_velocity > critical_capture_velocity)
/* Update particle velocities for reflection case (grabbed from fluent example and updated coeff */
{
#include "udf.h"
DEFINE_DPM_BC(bc_reflect,p,t,f,f_normal,dim)
{
real alpha; /* angle of particle path with face normal */
real vn=0.;
real nor_coeff = 1.;
real tan_coeff = 0.3;
real normal[3];
int i, idim = dim;
real NV_VEC(x);
#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);
/* Calculate the normal component, rescale its magnitude by
the coefficient of restitution and subtract the change */
/* Compute normal velocity. */
for(i=0; i<idim; i++)
vn += P_VEL(p)[i]*normal[i];
/* Subtract off normal velocity. */
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;
}
/* particle travels slow and is stuck*/
else
return PATH_ABORT;
}
} 


[`e`]

No worries, I find problem solving enjoyable.

Why is there DEFINE_DPM_BC(bc_reflect,p,t,f,f_normal,dim) nestled inside DEFINE_DPM_BC(my_dpm_bc,p,t,f,f_normal,dim), and the commented title does not have a trailing */?

[sjbub]

cause I'm a spaz and not paying attention...

I got this to compile on fluent with out any errors. So its definitely an understatement to say I'm excited. I did some more research for udf's online and spent a lot of time messing around with them. Lots of guess and check since like I said before I'm very green with this stuff.

PHP Code:

/**********************************************************************/

#include "udf.h"
#include "dpm.h"
DEFINE_DPM_BC(best_dpmbc,p,t,f,f_normal,dim)
{
    /*Variable declarations*/
    real capture_vel;
    real diameter;
    real vmag;  /*particle velocity magnitude*/
    int i;  /*working counter for script*/

    /*Varaible evaluations*/
    diameter=P_DIAM(p);

    /*Calculate capture velocity*/
    capture_vel=(0.000000515/diameter);

    /*Calculate velocity magnitude of particle*/
    for(i=0; i<dim; i++)
    {
        vmag += P_VEL(p)[i]*P_VEL(p)[i];
    }
    vmag = sqrt(vmag);

    /*If particle velocity is less than capture velocity, set to zero*/
    if(vmag < capture_vel)
    {
        for(i=0; i<dim; i++)
            P_VEL(p)[i] = 0;
    }
    /*Stop tracking the particle*/
    return PATH_ABORT;
}
/********************************************************************** 


/

[`e`]

Awesome, glad to hear you've got a working UDF now!

[why?]

Will the above posted UDF be OK to run in parallel to without any chances made to it? How can I tell if I need to make changes to a UDF to run it in parallel?

[`e`]

Your DPM BC should run fine in parallel because only a single process calls this macro when a particle collides with the relevant wall (particles are only tracked on a single process at a time from what I can determine).

[hwet]

Your DPM BC should run fine in parallel because only a single process calls this macro when a particle collides with the relevant wall (particles are only tracked on a single process at a time from what I can determine).

Is it possible that at some time more than one processes call this UDF. Since for me if I run the simulation for a high number of time steps the fluent particle data file resets.

This does not happen if there is no UDF or the number of time steps is small.

[`e`]

Originally Posted by hwet

Is it possible that at some time more than one processes call this UDF. Since for me if I run the simulation for a high number of time steps the fluent particle data file resets.

This does not happen if there is no UDF or the number of time steps is small.

Perhaps if multiple particles (on different processes) collide with this boundary at the same time. I've not experienced any issues with using DPM in parallel. However, I've not used the "particle data file" (but expect the developers to have coded for parallel).

[yingwen]

Originally Posted by `e`

In your case, the particles have two options (use an if statement) once they hit the wall:
(1) particle is captured by the wall ("stick" == trapped; using a return of PATH_ABORT)
(2) particle continues (using a return of PATH_ACTIVE) with reflected velocities. Either assume perfectly elastic conditions or use a relevant equation to calculate the new velocities.

If you simply return a particle from the boundary with the same velocities it had, then I suspect the particle would collide with the same boundary again, creating an infinite loop (try this!).

Excuse me, this is my first time to ask for help on CFD Online, I'm trouble with the problem for several days. I try to use the macro "DPM_BC", the particles can bounce off the wall or be trapped by the wall, but I don't know how to make them pass through the wall just like the "interior" boundary condition.

I tried to use "return PATH_ACTIVE", but the particles continued to hit the wall, and it seems that I can't change their position simply by changing the value of P_POS(p)[i].


I wonder that if u can give me some advice! thank u very much!

[pakk]

How do you imagine particles to 'pass through' your wall? What is on the other side of the wall?


Normally speaking, a wall is the end of your calculation domain. So a particle can not go beyond that wall, it would be outside of the calculation domain. Do you have a very specific situation in which it makes sense to let a particle go through a wall?

[yingwen]

Originally Posted by pakk

How do you imagine particles to 'pass through' your wall? What is on the other side of the wall?


Normally speaking, a wall is the end of your calculation domain. So a particle can not go beyond that wall, it would be outside of the calculation domain. Do you have a very specific situation in which it makes sense to let a particle go through a wall?

The wall is an internal wall, both side of the wall is fluid domain. You can imagine the gas with particles flows in a tube, and there is a net which makes the big particles stop and small ones pass through.

Thank u for your reply!