[sjbub]

I'm trying to have the particles "stick" or abort their paths at the wall with these particular boundary conditions but they are being labeled as escaped by fluent during the dpm. Anyone else seen this?

PHP Code:

/**********************************************************************/
/*UDF Particle Sticking Criteria*/
/*funz Theory*/
/**********************************************************************/

#include "udf.h"
#include "dpm.h"
DEFINE_DPM_BC(my_trapbc,p,t,f,f_normal,dim)
{
    /*Variable declarations*/
    real capture_vel;
    real diameter;
    real vmag;  /*particle velocity magnitude*/
    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];
    real NV_VEC(x);
    int i, idim = dim; /*working counter for script*/


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

    /*Calculate capture velocity*/
    capture_vel=1000;

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

    /*if vmag is greater than capvel then particle is fast and bounces*/
    if(vmag > capture_vel)
    {
        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 if (vmag < capture_vel)
        {
        for(i=0; i<dim; i++)
            P_VEL(p)[i] = 0;
        }
    /*Stop tracking the particle*/
    return PATH_ABORT;    
    }
}
/**********************************************************************/ 


[sjbub]

I believe it may be my else statement at the end.

[`e`]

If they're being reflected then that suggests that 'vmag' is always greater than 'capture_vel'. A velocity of 1000 m/s is about Mach 3 in air... what context has particles travelling at this speed?

[sjbub]

Originally Posted by `e`

If they're being reflected then that suggests that 'vmag' is always greater than 'capture_vel'. A velocity of 1000 m/s is about Mach 3 in air... what context has particles travelling at this speed?

It doesn't represent anything real, I wanted to try to force everything to bounce off the wall and exit. My previous udf had only about 2% of the mass exiting the domain and I was thinking that my udf wasn't telling the particles to return to their path if they didn't meet the capture velocity criteria.

[`e`]

If the capture velocity is less than the particle velocity at the wall then they should stick (according to PATH_ABORT). Try printing a message when the code reaches the line before aborting the particle path to check whether or not this line of code is being executed.

/*Stop tracking the particle*/
Message("A particle is about to be aborted...\n");
return PATH_ABORT;

Another possible cause is if you've forgotten to enable the user-defined boundary condition for your wall (because the default is reflect on walls and it sounds like that's happening).

[pakk]

Just a minor point: when your particles are smaller than the capture velocity, your code effectively does this:

for(i=0; i<dim; i++)
  P_VEL(p)[i] = 0;
}
return PATH_ABORT;

So first you set the velocity of the particle to zero, and then you tell Fluent to stop following the particle. You don't have to set the velocity of the particle to zero, you can just return PATH_ABORT and the particle will be stuck.

[mimi0201]

you just used exactly same code as i used, in determining the particle impacting on the wall!!! vmag cannot exceed 1000m/s right???

[sjbub]

Yes, true the velocity cannot be that high. I was testing my udf to force everything to stick by putting a velocity that I know would not be possible.

I had my original UDF set up incorrectly. I had the velocity capture part set up right but I forgot to tell the particle what to do (continue on its path) if it didnt meet the criteria.

[sjbub]

Originally Posted by pakk

Just a minor point: when your particles are smaller than the capture velocity, your code effectively does this:

for(i=0; i<dim; i++)
  P_VEL(p)[i] = 0;
}
return PATH_ABORT;

So first you set the velocity of the particle to zero, and then you tell Fluent to stop following the particle. You don't have to set the velocity of the particle to zero, you can just return PATH_ABORT and the particle will be stuck.

oh okay thanks pakk

[hwet]

I was going over the same example in the UDF manual today, could you tell what is the point of adding and then subtracting the same normal component from the particle velocity.
Also how does the above code actually make the particle reflect, i think it is just making the particle hit the wall again with a slightly smaller velocity again and again?

[`e`]

They're not adding and removing the normal velocity from the particle. They're evaluating the normal velocity of the particle, storing this value in the variable 'vn', and then removing this value from the particle (therefore the particle is motionless). Next, they apply a new reflected velocity to the particle based on 'vn' but scaled back by a coefficient of restitution.

[hwet]

Shouldn't the new reflected velocity have a negative component to change its direction (reflected), what they have done appears to be just reducing the velocity of the particle by a certain coefficient, that i think will not make it change its direction to be reflected, it will still be moving towards the wall just with a changed velocity and a different angle maybe?

[`e`]

Ok, let's work through this reflection code with an example.

Suppose our particle exists in 2-D and has a velocity vector p = (3,-2) m/s. The wall is parallel to the x-axis and therefore has a normal unit vector of n = (0,1). Units are absent to keep simplicity in the following text.

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

vn = p dot n = (3,-2) dot (0,1) = -2 (the same speed our particle was travelling towards the wall, "normal component")

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

p = p - vn*n = (3,-2) - (-2)*(0,1) = (3,0) (now the normal component for our particle is zero)

/* Apply tangential coefficient of restitution. */
for(i=0; i<idim; i++)
P_VEL(p)[i] *= tan_coeff;

Suppose the tangential coefficient is 0.5:
p = p*tangential_coefficient = (3,0)*0.5 = (1.5,0) (the tangential speed is reduced but the direction remains the same)

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

Suppose the normal coefficient is 0.5:
p = p - normal_coefficient*vn*n = (1.5,0) - 0.5*(-2)*(0,1) = (1.5,1) (the normal speed is reduced and changes direction)

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

Some DPM solvers may use the velocity of the particle as it enters the current cell for calculating the next velocity or position.

As we expect, our reflected velocity, (1.5,1) m/s, has been reduced and directed away from the wall. This code works.

[hwet]

Best explanation ever! Thanks

[hwet]

Could you tell why in the code originally posted above, the angle is calculated, i dont see it being use anywhere?

[`e`]

The code in the original post of this thread is copied from an example in the UDF Manual. Perhaps the developers writing the manual wanted to provide this code for calculating the angle in case a user wanted this value.

[hwet]

Hmm that makes sense.

Was just wondering if the angle of incidence is used in anyway in the reflected velocity calculation for this particular example, apparently it is calculated but not used.

Actually that example is for reflected particles as well, and explaining the procedure it states 'First, the angle of incidence is calculated. Next, the normal particle velocity, with respect to wall...' So i thought I was just missing out on something. Thanks

[`e`]

What might have happened, is developer A wrote the code using the angle of incidence (component-wise method). Then, a second developer B was tasked with optimising or reviewing this code and changed to using the dot product and normal vector but forgot to edit the comments and remove the angle calculation. Who knows, there's a few things in Fluent I don't understand their motivations for...

[mimi0201]

Im just wondering what does this sentence do in this code? "if ((NNULLP(t)) && (THREAD_TYPE(t) == THREAD_F_WALL))
F_CENTROID(x,f,t);"

[`e`]

F_CENTROID(x,f,t);

This macro returns the coordinates of the centroid of the face f belonging to thread t.

if ((NNULLP(t)) && (THREAD_TYPE(t) == THREAD_F_WALL))

This if statement restricts when the above macro is called. Specifically if this thread exists (has storage allocated, NNULLP(t) returns true) and if this thread is a wall boundary type (THREAD_TYPE(t) == THREAD_F_WALL).