Learning Forum

[kiran.purushothamakeshavan]

 

I think the issue was with the order of each force quantities, it was in the scale if E-09 or E-11 so I used a factor which multiplied these insignificant quantities to compare with the others.

anyway, now my question is simple, I get to see how many particles stick and how many bounce from the object. How do I visualise the contours for the same? I would want to display the particle adhesion on a surface as a heatmap.

here is my code:

 

#include “udf.h”

#include “dpm.h”

#define mu 0.5

#define PI 3.14

 

real FD(real diam, real velocidad)

{

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

}

FILE *fp1;

FILE *fp2;

FILE *fp3;

//real x[3];

//x[0] = 0.0;   // Assign value to the x component

//x[1] = 0.0;   // Assign value to the y component

//x[2] = 0.0;   // Assign value to the z component

DEFINE_DPM_DRAG(particle_drag_force, Re, p)

{

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

    if (Re < 1)

    {

        drag_force = 18;

    }

    else if (Re < 20.0)

    {

        w = log10(Re);

        drag_force = 18.0 + 2.367 * pow(Re, 0.82 – 0.05 * w);

    }

    else

    {

        drag_force = 18.0 + 3.483 * pow(Re, 0.6305);

    }

    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;

    // F = (1/2)*1.225*0.3*M_PI*(((P_DIAM(p))*(P_DIAM(p)))/4)*(dvtdn*dvtdn);

    // F = (1/2)*1.225*0.3*10000000*M_PI*(((P_DIAM(p))*(P_DIAM(p)))/4)*(dvtdn*dvtdn)*(10000000000);

    // P_USER_REAL(p, 2) = F;

    // nf = C_P(c, tcell) * M_PI (((P_DIAM(p))*(P_DIAM(p)))/4) *10000000000;

 

    fp1=fopen (“Impact-2-drag-data.txt” , “a”);

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

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

    return (drag_force);

}

DEFINE_DPM_BC(wall_model, p, t, f, f_normal, dim)

{

    real alpha;  /* angle of particle path with face normal */

    real vn=0.;

    real nor_coeff = 0.3;

    real tan_coeff = 0.3;

    real normal[3];

    int i;

    real posx, posy, posz;

    // real NV_VEC(x);

    real x[ND_ND];

   

 //#if !RP_HOST

 

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

        //Message (“Text \n”);

    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.

        {

            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;

            for(i=0; i

                P_VEL(p)[i] = 0;

            P_USER_REAL(p, 7) += 1.;

            fp2=fopen (“Impact-Stick.txt” , “a”);

            fprintf( fp2 , “Force due to aero (Fa): %f, Friction force (ff): %f, How many stick? = %f, x coordinate: %f, y coorfdinate: %f, z coordinate: %f\n”, 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_ACTIVE;

        }

        else

        {

            alpha = (PI/2.0) – 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);

            /* 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];

            P_USER_REAL(p, 8) += 1.;

            fp3=fopen (“Impact-Bounce.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;

        }

    }

    else

    {

        return PATH_END;

    }

 //#endif

}