Learning Forum

[kiran.purushothamakeshavan]

Hello,

I have written a UDF where I count the number of particles that hit a wall and gets stuck to it under certain conditions. 

I see that the number of particles hitting this wall is 22k while I only count 4.5k particles at the inlet. Why do you think this might be happening?

/* reflect boundary condition for inert particles */

#include

#include "udf.h"

#include "dpm.h"

#include "mem.h"

#include "sg.h"

#include "stdio.h"

#include "math.h"

#define g -9.81

#define mu 0.5

#define NUM_UDM 7

#define PI 3.141592654

FILE *fp1;

FILE *fp2;

FILE *fp3;

real FD(real diam, real velocidad)

{

  return ((1/2) * 1.225 * 0.44 * M_PI*(((diam)*(diam))/4) * (velocidad*velocidad));

}

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

  real posx, posy, posz;

  // real NV_VEC(x);

  real x[ND_ND];

  real A[ND_ND];

  Domain* d;

  real Area;

  real Flux_density, tem_MFR;

  real tem_MFR_2, Flux_density_2, Area_2;

  // FILE *fp1;

  real w, drag_force, F, nf, ff, ff123;

  Thread *tcell = P_CELL_THREAD(p);

  cell_t c = P_CELL(p);

  real relVel[ND_3], dvtdn, fluidVel[ND_3];

  // mass of particles

  real tem_particle_dia = 0;

  NV_D(fluidVel, =, C_U(c, tcell), C_V(c, tcell), C_W(c, tcell));     //assigning fluid velocity components to vector

  NV_VV(relVel, =, fluidVel, -, P_VEL(p));            //calculating relative velocity of flow on particle by subtracting absolute particle velocity from fluid velocity

  dvtdn = NV_MAG(relVel);                 //calculation of relative velocity magnitude

  // P_USER_REAL(p, 0) = FD(drag_force, C_MU_L(c, tcell), P_RHO(p), P_DIAM(p), P_MASS(p), dvtdn); /*Drag force*/

  P_USER_REAL(p, 0) = FD(P_DIAM(p), dvtdn); /*Drag force*/

  P_USER_REAL(p, 1) = dvtdn;

 

  // fprintf( fp1 , "%e %f %e %f %f %e %e\n", P_USER_REAL(p, 0), P_USER_REAL(p, 1), P_USER_REAL(p, 2), drag_force, P_DIAM(p), F, P_USER_REAL(p, 3));

  real diameter_squared = P_DIAM(p) * P_DIAM(p);

  real area_term = (diameter_squared) / 4;

  real velocity_term = dvtdn * dvtdn;

  real constant_factor = (1.0 / 2.0) * 1.225 * 0.3 * PI;

  F = constant_factor * area_term * velocity_term * (10000000);

  P_USER_REAL(p, 2) = F;

  nf = C_P(c, tcell) * (PI / 4.0) * diameter_squared *10000000;

  ff = mu * nf;

  ff123 = ff;

  P_USER_REAL(p, 3) = ff123;

  fp1=fopen ("Impact-2-drag-data-mod.txt" , "a");

  fprintf(fp1, "Area Term: %e, Velocity Term: %e, F (final): %e, Constant factor: %e, friction force: %e, Normal force: %e\n", area_term, velocity_term, P_USER_REAL(p, 2), constant_factor, P_USER_REAL(p, 3), nf);

  fprintf(fp1, "Pressure: %e,  ff123: %e\n", C_P(c, tcell), ff123);

  fclose(fp1);

  Message("Inside the first section");

 #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

            normal[i] = f_normal[i];

        }

    }

  else

 #endif

    for (i=0; i

      normal[i] = f_normal[i];

  if(p->type==DPM_TYPE_INERT)

    {

      if (fabs(P_USER_REAL(p, 2)) < fabs(P_USER_REAL(p, 3)))

    // Particle bounce: Aerodynamic force > Pressure force on particle. So particle must be resuspended to the flow.

      {

          for(i=0; i

            P_VEL(p)[i] = 0;

          P_USER_REAL(p, 7) += 1.;

            posx = P_POS(p)[0];

            posy = P_POS(p)[1];

            posz = P_POS(p)[2];

            P_USER_REAL(p, 4) = posx;

            P_USER_REAL(p, 5) = posy;

            P_USER_REAL(p, 6) = posz;

            C_UDMI(c, tcell, 0) += 1.0;

            C_UDMI(c, tcell, 1) += P_FLOW_RATE(p);

            //P_USER_REAL(p, 8) += P_FLOW_RATE(p);

            tem_MFR = P_FLOW_RATE(p);

            F_AREA(A,f,t);

            Area=NV_MAG(A);

            // why do you have tem_MFR instead of

            Flux_density = tem_MFR/ Area;

            C_UDMI(c, tcell, 2) = Flux_density;

          fp2=fopen ("Impact-Stick-mod.txt" , "a");

          fprintf(fp2 , "tem_MFR: %e, Area: %e, Particle flux density: %f, Force due to aero (Fa): %f, Friction force (ff): %f, How many stick? = %f, x coordinate: %f, y coorfdinate: %f, z coordinate: %f\n", tem_MFR, Area, Flux_density, P_USER_REAL(p, 2), P_USER_REAL(p, 3), P_USER_REAL(p, 7), P_USER_REAL(p, 4), P_USER_REAL(p, 5), P_USER_REAL(p, 6));

          fflush(fp2);

          fclose(fp2);

          return PATH_ABORT;

        }

      else

      {

        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

          vn += P_VEL(p)[i]*normal[i];

        /* Subtract off normal velocity. */

        for(i=0; i

          P_VEL(p)[i] -= vn*normal[i];

        /* Apply tangential coefficient of restitution. */

        for(i=0; i

          P_VEL(p)[i] *= tan_coeff;

        /* Add reflected normal velocity. */

        for(i=0; i

          P_VEL(p)[i] -= nor_coeff*vn*normal[i];  

        /* Store new velocity in P_VEL0 of particle */

        for(i=0; i

          P_VEL0(p)[i] = P_VEL(p)[i];

         

        C_UDMI(c, tcell, 3) += 1.0;

        C_UDMI(c, tcell, 4) += P_FLOW_RATE(p);

        P_USER_REAL(p, 8) += P_FLOW_RATE(p);

        tem_MFR_2 = P_FLOW_RATE(p);

        F_AREA(A,f,t);

        Area_2=NV_MAG(A);

        Flux_density_2 = tem_MFR_2/ Area_2;

        C_UDMI(c, tcell, 5) = Flux_density_2;

        P_USER_REAL(p, 9) += 1.;

        fp3=fopen ("Impact-Bounce-mod.txt" , "a") ;

        fprintf(fp3 ,"Force due to aero (Fa): %f, Friction force (ff): %f, Number of particles bounced: %f\n", P_USER_REAL(p, 2), P_USER_REAL(p, 3), P_USER_REAL(p, 8));

        fflush(fp3);

        fclose(fp3);

        return PATH_ACTIVE;

      }

  }

  return PATH_ABORT;

}

 

 

I tried debugging this but I do not seem to understand why this might happen, your suggestions would be greatly appreciated.

 

Thanks,
Kiran

[Rob]

You're tracking parcels, not particles, but that's a definition issue. Is the model transient? Have you run the particle tracking more than once? Note, I'm aware that particle and parcel are used interchangeably (and incorrectly) within Fluent, but as that's historic won't be changing anytime soon. 

[kiran.purushothamakeshavan]

If its a parcel, how would I know how many particles are being released at the inlet?

Yes this is a transient model. I do not think I ran the particle tracking more than once. 
This is my DPM setting:

Thanks,
Kiran

[Rob]

You're tracking parcels, each parcel follows the trajectory of the particle but has mass etc of the parcel. Check the DPM part of the manuals and possibly some background reading. Note, a parcel usually has significantly more mass than a single particle but that's not always the case. In Fluent particle is used incorrectly in place of parcel in many instances. 

I don't know what your injection is set as, but you're releasing some parcels every flow time step. If it's a surface release from the inlet you're releasing a parcel for every facet on the inlet every fluid time step. That may explain the parcel count you're seeing.