You can use DEFINE_DPM_BC to specify your own boundary conditions for particles. The function is executed every time a particle touches a boundary of the domain, except for symmetric or periodic boundaries. You can define a separate UDF (using DEFINE_DPM_BC) for each boundary.

DEFINE_DPM_BC (name, tp, t, f, f_normal, dim)

Function returns

int

There are six arguments to DEFINE_DPM_BC: name, tp, t, f, f_normal, and dim. You supply name, the name of the UDF. tp, t, f, f_normal, and dim are variables that are passed by the Ansys Fluent solver to your UDF. Your UDF will need to compute the new velocity of a particle after hitting the wall, and then return the status of the particle track (as an int), after it has hit the wall.

For the return status PATH_ACTIVE, the particle continues to track. For the return status PATH_ABORT, the particle will be stopped and considered to be aborted. For the return status PATH_END, the particle will be stopped as well, but considered to have escaped from the domain.

This example shows the usage of DEFINE_DPM_BC for a simple reflection at walls. It is similar to the reflection method executed by Ansys Fluent except that Ansys Fluent accommodates moving walls. The function must be executed as a compiled UDF.

The function assumes an ideal reflection for the normal velocity component (nor_coeff 1) while the tangential component is damped (tan_coeff 0.3). First, the angle of incidence is computed. Next, the normal particle velocity, with respect to the wall, is computed and subtracted from the particles velocity. The reflection is complete after the reflected normal velocity is added. The new particle velocity has to be stored in TP_VEL0 to account for the change of particle velocity in the momentum balance for coupled flows. The function returns PATH_ACTIVE for inert particles while it stops particles of all other types.

/* reflect boundary condition for inert particles */
 #include "udf.h"
 DEFINE_DPM_BC(bc_reflect,tp,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(TP_POS(tp)[1]*TP_POS(tp)[1] +
              TP_POS(tp)[2]*TP_POS(tp)[2]);
         if (R > 1.e-20)
           {
              idim = 3;
              normal[0] = f_normal[0];
              normal[1] = (f_normal[1]*TP_POS(tp)[1])/R;
              normal[2] = (f_normal[1]*TP_POS(tp)[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(TP_TYPE(tp) == DPM_TYPE_INERT)
      {
         alpha = M_PI/2. - acos(MAX(-1.,MIN(1.,NV_DOT(normal,TP_VEL(tp))/
              MAX(NV_MAG(TP_VEL(tp)),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 += TP_VEL(tp)[i]*normal[i];

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

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

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

         /* Store new velocity in TP_VEL0 of particle */
         for(i=0; i<idim; i++)
           TP_VEL0(tp)[i] = TP_VEL(tp)[i];

         return PATH_ACTIVE;
      }
  return PATH_ABORT;
 } 

This example shows how to use DEFINE_DPM_BC for a wall impingement model. The function must be executed as a compiled UDF.

 #include "udf.h"
 #include "dpm.h"
 #include "surf.h"
 #include "random.h"
 
 /* define a user-defined dpm boundary condition routine
 * bc_reflect: name
 * tp:  the tracked particle
 * t:   the touched face thread
 * f:   the touched face
 * f_normal:  normal vector of touched face
 * dim:  dimension of the problem (2 in 2d and 2d-axi-swirl, 3 in 3d)
 *
 * return is the status of the particle, see enumeration of Path_Status
 * in dpm.h
 */
 
 #define V_CROSS(a,b,r)\
     ((r)[0] = (a)[1]*(b)[2] - (b)[1]*(a)[2],\
     (r)[1] = (a)[2]*(b)[0] - (b)[2]*(a)[0],\
     (r)[2] = (a)[0]*(b)[1] - (b)[0]*(a)[1])
  
 DEFINE_DPM_BC(bc_wall_jet, tp, thread, f, f_normal, dim)
 {
    /*
      Routine implementing the Naber and Reitz Wall
      impingement model (SAE 880107)
    */
    real normal[3];
    real tan_1[3];
    real tan_2[3];
    real rel_vel[3];
    real face_vel[3];
    real alpha, beta, phi, cp, sp;
    real rel_dot_n, vmag, vnew, dum;
    real weber_in, weber_out;
    int i, idim = dim;
    cxboolean moving = (SV_ALLOCATED_P (thread,SV_BND_GRID_V) &&
    SV_ALLOCATED_P (thread,SV_WALL_V));
 
 #if RP_2D
    if (rp_axi_swirl)
      {
         real R = sqrt(TP_POS(tp)[1]*TP_POS(tp)[1] +
            TP_POS(tp)[2]*TP_POS(tp)[2]);
         if (R > 1.e-20)
           {
              idim = 3;
              normal[0] = f_normal[0];
              normal[1] = (f_normal[1]*TP_POS(tp)[1])/R;
              normal[2] = (f_normal[1]*TP_POS(tp)[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];

      /*
        Set up velocity vectors and calculate the Weber number
        to determine the regime.
      */

    for (i=0; i<idim; i++)
      {
         if (moving)
           face_vel[i] = WALL_F_VV(f,thread)[i] + BOUNDARY_F_GRID_VV(f,thread)[i];
         else
           face_vel[i] = 0.0;
           rel_vel[i] = TP_VEL(tp)[i] - face_vel[i];
      }

    vmag = MAX(NV_MAG(rel_vel),DPM_SMALL);
    rel_dot_n = MAX(NV_DOT(rel_vel,normal),DPM_SMALL);
    weber_in = TP_RHO(tp) * SQR(rel_dot_n) * TP_DIAM(tp) /
      MAX(DPM_SURFTEN(tp), DPM_SMALL);

    /*
    Regime where bouncing occurs (We_in < 80).
    (Data from Mundo, Sommerfeld and Tropea
     Int. J. of Multiphase Flow, v21, #2, pp151-173, 1995)  */

    if (weber_in <= 80.)
      {
         weber_out = 0.6785*weber_in*exp(-0.04415*weber_in);
         vnew = rel_dot_n * (1.0 + sqrt(weber_out /
         MAX(weber_in, DPM_SMALL)));

     /*
     The normal component of the velocity is changed based
     on the experimental paper above (that is the Weber number
     is based on the relative velocity).
       */

         for (i=0; i < idim; i++)
           TP_VEL(tp)[i] = rel_vel[i] - vnew*normal[i] + face_vel[i];
      }

    if (weber_in > 80.)
      {
         alpha = acos(-rel_dot_n/vmag);

     /*
     Get one tangent vector by subtracting off the normal
     component from the impingement vector, then cross the
     normal with the tangent to get an out-of-plane vector.
     */

     for (i=0; i < idim; i++)
       tan_1[i] = rel_vel[i] - rel_dot_n*normal[i];

     UNIT_VECT(tan_1,tan_1);
     V_CROSS(tan_1,normal,tan_2);

     /*
     beta is calculated by neglecting the coth(alpha)   
     term in the paper (it is approximately right).
     */

     beta = MAX(M_PI*sqrt(sin(alpha)/(1.0-sin(alpha))),DPM_SMALL);
     phi= -M_PI/beta*log(1.0-cheap_uniform_random()*(1.0-exp(-beta)));
     if (cheap_uniform_random() > 0.5)
       phi = -phi;

     vnew = vmag;

     cp = cos(phi);
     sp = sin(phi);

     for (i=0; i < idim; i++)
       TP_VEL(tp)[i] = vnew*(tan_1[i]*cp + tan_2[i]*sp) + face_vel[i];

      }

  /*
    Subtract off from the original state.
  */
    for (i=0; i < idim; i++)
        TP_VEL0(tp)[i] = TP_VEL(tp)[i];

    if (DPM_STOCHASTIC_P(TP_INJECTION(tp)))
      {
         /* Reflect turbulent fluctuations also */
         /* Compute normal velocity. */

       dum = 0;
       for(i=0; i<idim; i++)
         dum += tp->V_prime[i]*normal[i];

           /* Subtract off normal velocity. */

        for(i=0; i<idim; i++)
           tp->V_prime[i] -= 2.*dum*normal[i];
      }
    return PATH_ACTIVE;
 } 

This example shows how to use DEFINE_DPM_BC for a simple filter model. It illustrates the proper way to handle particles that cross the boundary zone and enter a new cell. If the particle leaves the current cell, this information is stored in TP_CCELL(tp).

#include "udf.h"
#include "cxsurf.h"


/**
 *  filter particles that are above a specified cut-off size; allow others to pass
 **/

#define CUT_OFF_SIZE 10.0e-6   /* define the cut-off size in meters */

DEFINE_DPM_BC(my_filter, tp, t, f, f_normal, dim)
{
  Thread *t0 = THREAD_T0(t);
  Thread *t1 = THREAD_T1(t);
  cxboolean in_t0;
  CX_Cell_Id neigh;
  
  /* determine what cell/thread the particle is currently in */
  if ( (P_CELL_THREAD(tp)->id == t0->id) && (P_CELL(tp) == F_C0(f, t)) )
    in_t0 = TRUE;
  else
    in_t0 = FALSE;
  
  /* filter particles above cut-off size */
  if (P_DIAM(tp) > CUT_OFF_SIZE)
    return PATH_END;

  /* particle diameter is below cut-off size - set the new cell and thread */
  if (in_t0)
    STORE_CX_CELL(&neigh, F_C1(f, t), t1);
  else
    STORE_CX_CELL(&neigh, F_C0(f, t), t0);

  COPY_CELL_ID(TP_CCELL(tp), &neigh);
  calcFaceEquations(tp);

  return PATH_ACTIVE;
}

This example shows how to use DEFINE_DPM_BC to reflect a Lagrangian wall film particle off of a boundary for cases that use the subtet particle tracking algorithm. If your UDF modifies one of the particle macros (mainly, TP_POS(tp), TP_CCELL(tp), TP_FILM_FACE(tp), or TP_FILM_THREAD(tp)), you must ensure that the values of these particle variables are consistent with each other. For example, the particle position, TP_POS(tp), must be within the cell specified by TP_CCELL(tp), and the cell face specified by TP_FILM_FACE(tp) must be one of the faces of that cell. If the UDF does not modify any of these macros, then the current values will still be valid after the UDF is executed.

#include "udf.h"


/**
 *  reflect wall film particles off of a boundary when using subtet particle tracking
 **/


DEFINE_DPM_BC(my_reflect_bc, tp, tf, f, wall_normal, dim)
{
  int i;

  double dot_product;
  double normal[3];          /* wall film face normal */
  real edge_normal[3];       /* unit edge normal */

#if RP_2D
  normal[2] = 0.0;
#endif

  /* this UDF is only valid when using subtet tracking */
  if ( !dpm_par.use_subtet_tracking )
    return PATH_ABORT;

  /* get the wall film face normal */
  memcpy(normal, tp->gvtp.subtet->normal, 3 * sizeof(double));

  /* get the dot product of the film face normal and the reflected wall normal */
  dot_product = NV_DOT(normal, wall_normal);

  /* get the projection of the wall normal onto the wall film face normal*/
  for (i = 0; i < 3; i++)
    normal[i] *= dot_product;

  /* subtract the normal component from the original wall normal */
  for (i = 0; i < 3; i++)
    edge_normal[i] = wall_normal[i] - (real) normal[i];  /* in-plane (edge) normal */
  
  /* reflect the particle in the plane of the wall film face */
  UNIT_VECT(edge_normal, edge_normal);
  Reflect_Particle(tp, edge_normal, NULL, f, tf, NULL, -1);

  return PATH_ACTIVE;
}

The following example shows the usage of the DEFINE_DPM_BC user-defined function (UDF) to:

#include "udf.h"
DEFINE_DPM_BC(reinj, tp, t, f, normal, dim)
{
  /* prepare some variables used below:
   */
  cxboolean reinj_unsteady_tracking = dpm_par.unsteady_tracking;
  cxboolean reinj_update = TRUE;
  cxboolean reinj_update_summary = (((0 == c_par.skip) &&  ! reinj_unsteady_tracking) ||
                                    (reinj_unsteady_tracking && reinj_update)) &&
                                   (dpm_par.is_pathline == 0);


  /* Determine which injection to use for the reinjection --
   * in this simple example, particles whose Y coordinate (TP_POS(tp)[1])
   * exceeds 0.02_m will use injection "single-inj-no-mfr", all others "inlet-inj-no-mfr":
   */
  char *this_reinj_inj_name = TP_POS(tp)[1] > 0.02 ? "single-inj-no-mfr" : "inlet-inj-no-mfr";
  Injection *this_reinj_inj = Pick_Injection(this_reinj_inj_name);


  /* Now "this_reinj_inj" should be a pointer to the injection struct.
   * Let us check whether we've actually managed to choose one:
   */
  if (NULLP(this_reinj_inj))
  {
    DELAY_ERROR("reinj UDF: cannot find injection: %s", this_reinj_inj_name);
    return PATH_END;     /* Fluent will call "Escape_Tracked_Particle(..) */
  }


  /* OK, some injection was found -- assign it to the Tracked Particle:
   */
  tp->reinj_inj = this_reinj_inj;


  /* Add the particle to the dpm summary report for this boundary zone:
   */
  if (reinj_update_summary)
    add_to_dpm_summary(tp, FATE_BC_UDF, t);


  /* Have Fluent do some things that ought to be done
   * for every particle that escapes or is reinjected:
   */
  Escape_OR_Reinject_Tracked_Particle(tp, reinj_update_summary, f, t);


  /* EXAMPLE: Delay the reinjection by 0.5_s:
   */
  TP_TIME(tp) += 0.5;
  PP_TIME(tp->pp) += 0.5;

  return PATH_REINJECT;
}

This function performs the following:

The UDF shows how you can delay particle reinjection, thus accounting for a significant residence time of particles inside a recirculation line outside the meshed CFD domain. To achieve this, the UDF adds the time offset to both TP_TIME(tp) and PP_TIME(tp->pp). The particles will be moved to their new position and assigned their new velocity and a start time, and until the simulation reaches that time, the particles will remain still at their positions and not participate in the simulated processes. Note that you can also specify Reinjection Time Delay during particle reinjection in the Set Injection Properties dialog box (Point Properties tab). See Injection Types in the Fluent User's Guide for more information.

After the UDF that you have defined using DEFINE_DPM_BC is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the appropriate boundary condition dialog box (for example, the Velocity Inlet dialog box) in Ansys Fluent. See Hooking DEFINE_DPM_BC UDFs for details on how to hook your DEFINE_DPM_BC UDF to Ansys Fluent.

You can use DEFINE_DPM_BODY_FORCE to specify a body force other than a gravitational or drag force on the particles.

DEFINE_DPM_BODY_FORCE (name, tp, i)

Function returns

real

There are three arguments to DEFINE_DPM_BODY_FORCE: name, tp, and i. You supply name, the name of the UDF. tp and i are variables that are passed by the Ansys Fluent solver to your UDF. Your UDF will need to return the real value of the acceleration due to the body force (in m/s²) to the Ansys Fluent solver.

The following UDF computes the electrostatic force on a charged particle. It consists of two functions: one that enables the continuous interpolation of the electrostatic field to the particle position, and another that calculates the acceleration of the particle due to the electric field. This example assumes that the electric field components at the cell centroid are stored in user-defined memory.

DEFINE_DPM_BODY_FORCE is called at every particle time step in Ansys Fluent and requires a significant amount of CPU time to execute. For this reason, the UDF should be executed as a compiled UDF.

#include "udf.h"
#include "dpm_interpolation.h"

/* setup interpolation of the electrostatic field */
DEFINE_DPM_INTERPOLATION(dpm_setup_interp, tp)
{
  /* interpolate the electrostatic field to the particle position */
  if ( HIGH_RESOLUTION_TRACKING_ENABLED )
  {
    
    if (N_UDM < 7)
      return;
    
      /* specify now many and which UDM should be interpolated */
    int n_udm = 3;
    int udm_components[] = { 4, 5, 6 };

    /* interpolate UDM 4, 5, and 6 to the particle position */
    DPM_Interpolate_UDM_from_Nodes(tp, n_udm, udm_components);
  }
}

/* calculate the electrostatic force (acceleration) on charged particles */
DEFINE_DPM_BODY_FORCE(particle_E_force, tp, i)
{
  /* From Bell and Hockberg correlation at 30k RPM and 70kV */
  const real chargedensity = 3.884e-3;   /* C/kg */
  real electric_field;
  
  if (N_UDM < 7)
    return 0.0;

  if ( HIGH_RESOLUTION_TRACKING_ENABLED )
    electric_field = TP_UDMI(tp, 4 + i);  /* E-field at the particle position */
  else
  {
    Thread *t = TP_CELL_THREAD(tp);
    cell_t c  = TP_CELL(tp);

    electric_field = C_UDMI(c, t, 4 + i); /* E-field at the cell centroid */
  }
  
  /* return the acceleration due to the electric force */
  return electric_field * chargedensity;
}

After the UDF that you have defined using DEFINE_DPM_BODY_FORCE is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_DPM_BODY_FORCE UDFs for details on how to hook your DEFINE_DPM_BODY_FORCE UDF to Ansys Fluent.

You can use DEFINE_DPM_DRAG to specify the drag between particles and fluid as a dimensionless group () as it appears in the drag force per unit particle mass:

(2–19)

DEFINE_DPM_DRAG (name, Re, tp)

Function returns

real

There are three arguments to DEFINE_DPM_DRAG: name, Re, and tp. You supply name, the name of the UDF. Re and tp are variables that are passed by the Ansys Fluent solver to your UDF. Your UDF will need to return the real value of the drag force on a particle. The value returned to the solver must be dimensionless and represent 18 * Cd * Re / 24.

The following UDF, named particle_drag_force, computes the drag force on a particle and is a variation of the body force UDF presented in DEFINE_DPM_BODY_FORCE . The flow is the same, but a different curve is used to describe the particle drag. DEFINE_DPM_DRAG is called at every particle time step in Ansys Fluent, and requires a significant amount of CPU time to execute. For this reason, the UDF should be executed as a compiled UDF.

/***********************************************************************
   UDF for computing particle drag coefficient (18 Cd Re/24)
   curve as suggested by R. Clift, J. R. Grace and M. E. Weber
   "Bubbles, Drops, and Particles" (1978)
 ************************************************************************/
 
 #include "udf.h"
 
 DEFINE_DPM_DRAG(particle_drag_force,Re,tp)
 {
    real w, drag_force;
    if (Re < 0.01)
    {
      drag_force=18.0;
      return (drag_force);
    }
    else if (Re < 20.0)
    {
      w = log10(Re);
      drag_force = 18.0 + 2.367*pow(Re,0.82-0.05*w) ;
      return (drag_force);
    }
    else
    /* Note: suggested valid range 20 < Re < 260 */
    {
      drag_force = 18.0 + 3.483*pow(Re,0.6305) ;
      return (drag_force);
    }
 } 

After the UDF that you have defined using DEFINE_DPM_DRAG is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Set Injection Properties dialog box in Ansys Fluent. See Hooking DEFINE_DPM_DRAG UDFs for details on how to hook your DEFINE_DPM_DRAG UDF to Ansys Fluent.

You can use DEFINE_DPM_EROSION to specify the erosion and accretion rates calculated as the particle stream strikes a wall surface or hits a porous jump surface. The function is called when the particle encounters a reflecting surface.

DEFINE_DPM_EROSION (name, tp, t, f, normal, alpha, Vmag, mdot)

Function returns

void

There are eight arguments to DEFINE_DPM_EROSION: name, tp, t, f, normal, alpha, Vmag, and mdot. You supply name, the name of the UDF. tp, t, f, normal, alpha, Vmag, and mdot are variables that are passed by the Ansys Fluent solver to your UDF. Your UDF will need to compute the values for the erosion rate and/or accretion rate and store the values at the faces in F_STORAGE_R_XV(f, t, SV_DPMS_EROSION, EROSION_UDF) and F_STORAGE_R (f,t,SV_DPMS_ACCRETION), respectively.

The following is an example of a compiled UDF that uses DEFINE_DPM_EROSION to extend postprocessing of wall impacts in a 2D axisymmetric flow. It provides additional information on how the local particle deposition rate depends on the diameter and normal velocity of the particles. It is based on the assumption that every wall impact leads to more accretion, and, therefore, every trajectory is removed at its first wall impact. (This is done by setting the flag P_FL_REMOVED within DEFINE_DPM_EROSION.) User-defined memory locations (UDMLs) are used to store and visualize the following:

/***********************************************************************
   UDF for extending postprocessing of wall impacts
 ************************************************************************/
 #include "udf.h"
 
#define MIN_IMPACT_VELO -1000.
/* Minimum particle velocity normal to wall (m/s) to allow Accretion.*/
 
Domain *domain; /* Get the domain pointer and assign it later to domain*/
 
enum   /* Enumeration of used User-Defined Memory Locations. */
  {
    NUM_OF_HITS, /* Number of particle hits into wall face considered.*/
    AVG_DIAMETER, /* Average diameter of particles that hit the wall. */
    AVG_RADI_VELO, /* Average radial velocity of "" "" ------------ */
    NUM_OF_USED_UDM
  };
 
int UDM_checked = 0;  /* Availability of UDMLs checked? */
 
void reset_UDM_s(void); /* Function to follow below. */
 
int check_for_UDM(void)  /* Check for UDMLs availability... */
{
  Thread *t;
  if (UDM_checked)
    return UDM_checked;
 
  thread_loop_c(t,domain) /* We require all cell threads to */
  {       /* provide space in memory for UDML */
     if (FLUID_THREAD_P(t))   if (NULLP(THREAD_STORAGE(t,SV_UDM_I)))
        return 0;
  }
  UDM_checked = 1;  /* To make the following work properly... */
  reset_UDM_s();   /* This line will be executed only once, */
  return UDM_checked; /* because check_for_UDM checks for */
}      /* UDM_checked first. */
 
void reset_UDM_s(void)
{
  Thread *t;
  cell_t c;
  face_t f;
  int  i;
  if (!check_for_UDM()) /* Dont do it, if memory is not available. */
    return;
  Message("Resetting User Defined Memory...\n");
  thread_loop_f(t, domain)
    {
      if (NNULLP(THREAD_STORAGE(t,SV_UDM_I)))
        {
          begin_f_loop(f,t)
	         {
	           for (i = 0; i < NUM_OF_USED_UDM; i++)
		           F_UDMI(f,t,i) = 0.;
	         }
          end_f_loop(f, t)
        }
     else
       {
	         Message("Skipping FACE thread no. %d..\n", THREAD_ID(t));
       }
    }
  thread_loop_c(t,domain)
    {
      if (NNULLP(THREAD_STORAGE(t,SV_UDM_I)))
	       {
          begin_c_loop(c,t)
	         {
	           for (i = 0; i < NUM_OF_USED_UDM; i++)
		           C_UDMI(c,t,i) = 0.;
	         }
          end_c_loop(c,t)
	       }
      else
	       {
          Message(" Skipping CELL thread no. %d..\n", THREAD_ID(t));
	       }
    }    /* Skipping Cell Threads can happen if the user */
  /* uses reset_UDM prior to initializing. */
  Message(" --- Done.\n");
}
 
DEFINE_DPM_EROSION(dpm_accr, tp, t, f, normal, alpha, Vmag, Mdot)
{
  real A[ND_ND], area;
  int num_in_data;
  Thread *t0;
  cell_t c0;
  real imp_vel[3], vel_ortho;

#if RP_2D  
  if (rp_axi) 
    {
      real radi_pos[3], radius;
      /* The following is ONLY valid for 2d-axisymmetric calculations!!! */
      /* Additional effort is necessary because DPM tracking is done in */
      /* THREE dimensions for TWO-dimensional axisymmetric calculations. */
      
      radi_pos[0] = TP_POS(tp)[1];  /* Radial location vector. */
      radi_pos[1] = TP_POS(tp)[2];  /* (Y and Z in 0 and 1...) */
      radius = NV_MAG(radi_pos);
      NV_VS(radi_pos, =, radi_pos, /, radius);
      /* Normalized radius direction vector.*/
      imp_vel[0] = TP_VEL(tp)[0]; /* Axial particle velocity component. */
      imp_vel[1] = NVD_DOT(radi_pos, TP_VEL(tp)[1], TP_VEL(tp)[2], 0.);
    }
  else
#endif
    NV_V(imp_vel, =, TP_VEL(tp));
  
  /* Dot product of normalized radius vector and y & z components */
  /* of particle velocity vector gives _radial_ particle velocity */
  /* component */
  vel_ortho = NV_DOT(imp_vel, normal); /*velocity orthogonal to wall */

  if (vel_ortho < MIN_IMPACT_VELO) /* See above, MIN_IMPACT_VELO */
    return;
  
  if (!UDM_checked)   /* We will need some UDMs, */
    if (!check_for_UDM()) /* so check for their availability.. */
      return;    /* (Using int variable for speed, could */
  /* even just call check_for UDFM().) */
  
  c0 = F_C0(f,t);
  t0 = THREAD_T0(t);
  
  F_AREA(A,f,t);
  area = NV_MAG(A);
  F_STORAGE_R(f,t,SV_DPMS_ACCRETION) += Mdot / area;
  
  MARK_TP(tp, P_FL_REMOVED); /* Remove particle after accretion */

  /* F_UDMI not allocated for porous jumps */
  if (THREAD_TYPE(t) == THREAD_F_JUMP)
    return;

  num_in_data  = F_UDMI(f,t,NUM_OF_HITS);
  
  /* Average diameter of particles that hit the particular wall face:*/
  F_UDMI(f,t,AVG_DIAMETER) = (TP_DIAM(tp)
			      + num_in_data * F_UDMI(f,t,AVG_DIAMETER))
                              / (num_in_data + 1);
  C_UDMI(c0,t0,AVG_DIAMETER) = F_UDMI(f,t,AVG_DIAMETER);
 
  /* Average velocity normal to wall of particles hitting the wall:*/
  F_UDMI(f,t,AVG_RADI_VELO) = (vel_ortho
			       + num_in_data * F_UDMI(f,t,AVG_RADI_VELO))
                               / (num_in_data + 1);
  C_UDMI(c0,t0,AVG_RADI_VELO) = F_UDMI(f,t,AVG_RADI_VELO);
  
  F_UDMI(f, t, NUM_OF_HITS) = num_in_data + 1;
  C_UDMI(c0,t0,NUM_OF_HITS) = num_in_data + 1;

  
}
 
DEFINE_ON_DEMAND(reset_UDM)
{
  /* assign domain pointer with global domain */
  domain = Get_Domain(1);
  reset_UDM_s();
} 

After the UDF that you have defined using DEFINE_DPM_EROSION is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_DPM_EROSION UDFs for details on how to hook your DEFINE_DPM_EROSION UDF to Ansys Fluent.

You can use DEFINE_DPM_HEAT_MASS to specify the heat and mass transfer of multicomponent particles to the gas phase.

When a DEFINE_DPM_HEAT_MASS UDF is activated, then the number of species that can be referenced and interact with the particles via the UDF is limited to those with a species index less than the maximum UDF species number, defined using the TUI command define/models/dpm/options/maximum-udf-species. The default number for maximum-udf-species is 50.

DEFINE_DPM_HEAT_MASS (name, tp, C_p, hgas, hvap, cvap_surf, Z, dydt, dzdt)

Function returns

void

There are eight arguments to DEFINE_DPM_HEAT_MASS: name, tp, C_p, hgas, hvap, cvap_surf, Z, dydt, and dzdt. You supply name, the name of the UDF. tp, C_p, hgas, hvap, cvap_surf, and Z are variables that are passed by the Ansys Fluent solver to your UDF. Your UDF will need to compute the particle and gas phase source terms and store the values in dydt and dzdt, respectively.

The following is an example of a compiled UDF that uses DEFINE_DPM_HEAT_MASS. It implements the source terms for the following:

/***********************************************************************
   UDF for defining the heat and mass transport for
   multicomponent particle vaporization
 ***********************************************************************/
 #include "udf.h"
 #include "dpm_mem.h"
 
 DEFINE_DPM_HEAT_MASS(multivap,tp,Cp,hgas,hvap,cvap_surf,Z,dydt,dzdt)
 {
   int ns;
   Material *sp;
   real dens_total = 0.0;     /* total vapor density*/
   real P_total = 0.0;      /* vapor pressure */
   int nc = TP_N_COMPONENTS(tp);   /* number of particle components */
   Thread *t0 = TP_CELL_THREAD(tp);   /* thread where the particle is in*/
   Material *gas_mix = THREAD_MATERIAL(DPM_THREAD(t0, tp)); /* gas mixture
   material */
   Material *cond_mix = TP_MATERIAL(tp); /* particle mixture material*/
   cphase_state_t *c = &(tp->cphase[0]); /* cell information of particle location*/
   real molwt[MAX_SPE_EQNS]; /* molecular weight of gas species */
   real Tp = TP_T(tp);   /* particle temperature */
   real mp = TP_MASS(tp);   /* particle mass */
   real molwt_bulk = 0.;  /* average molecular weight in bulk gas */
   real Dp = DPM_DIAM_FROM_VOL(mp / TP_RHO(tp)); /* particle diameter */
   real Ap = DPM_AREA(Dp);      /* particle surface */
   real Pr = c->sHeat * c->mu / c->tCond;   /* Prandtl number */
   real Nu = 2.0 + 0.6 * sqrt(tp->Re) * pow(Pr, 1./3.); /* Nusselt number */
   real h = Nu * c->tCond / Dp;     /* Heat transfer coefficient*/
   real dh_dt = h * (c->temp - Tp) * Ap;  /* heat source term*/
   dydt[0] += dh_dt / (mp * Cp);
   dzdt->energy -= dh_dt;
   mixture_species_loop(gas_mix,sp,ns)
   {
      molwt[ns] = MATERIAL_PROP(sp,PROP_mwi); /* molecular weight of gas
         species */
    molwt_bulk += c->yi[ns] / molwt[ns]; /* average molecular weight */
   }

 /* prevent division by zero */
 molwt_bulk = MAX(molwt_bulk,DPM_SMALL);

 for (ns = 0; ns < nc; ns++)
   {
      int gas_index = TP_COMPONENT_INDEX_I(tp,ns);  /* gas species index of
       vaporization */
  if(gas_index >= 0)
    {
       /* condensed material */
       Material * cond_c = MIXTURE_COMPONENT(cond_mix, ns);
       /* vaporization temperature */
       real vap_temp = MATERIAL_PROP(cond_c,PROP_vap_temp);
       /* diffusion coefficient */
       real D = DPM_BINARY_DIFFUSIVITY(tp,cond_c,TP_T(tp));
       /* Schmidt number */
       real Sc = c->mu / (c->rho * D);
       /* mass transfer coefficient */
       real k = (2. + 0.6 * sqrt(tp->Re) * pow(Sc, 1./3.)) * D / Dp;
       /* bulk gas concentration (ideal gas) */
       real cvap_bulk = c->pressure / UNIVERSAL_GAS_CONSTANT / c->temp
       * c->yi[gas_index] / molwt_bulk / solver_par.molWeight[gas_index];
       /* vaporization rate */
       real vap_rate = k * molwt[gas_index] * Ap
       * (cvap_surf[ns] - cvap_bulk);
       /* no vaporization below vaporization temperature, no condensation */
       if (Tp < vap_temp || vap_rate < 0.0)
         vap_rate = 0.;

       dydt[1+ns] -= vap_rate;   
       dzdt->species[gas_index] += vap_rate;
       /* dT/dt = dh/dt / (m Cp)*/
       dydt[0] -= hvap[gas_index] * vap_rate / (mp * Cp);
       /* gas enthalpy source term */
       dzdt->energy += hgas[gas_index] * vap_rate;

       P_total += cvap_surf[ns];
       dens_total += cvap_surf[ns] * molwt[gas_index];
      }
   }
   /* multicomponent boiling */
   P_total *= Z * UNIVERSAL_GAS_CONSTANT * Tp;
   if (P_total > c->pressure && dydt[0] > 0.)
     {
        real h_boil = dydt[0] * mp * Cp;
        /* keep particle temperature constant */
        dydt[0] = 0.;
        for (ns = 0; ns < nc; ns++)
          {
             int gas_index = TP_COMPONENT_INDEX_I(tp,ns);
             if (gas_index >= 0)
               {
                  real boil_rate = h_boil / hvap[gas_index] * cvap_surf[ns] *
                     molwt[gas_index] / dens_total;
                  /* particle component mass source term */
                  dydt[1+ns] -= boil_rate;
                  /* fluid species source */
                  dzdt->species[gas_index] += boil_rate;
                  /* fluid energy source */
                  dzdt->energy += hgas[gas_index] * boil_rate;
               }
          }
     }
 }

After the UDF that you have defined using DEFINE_DPM_HEAT_MASS is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument (for example, multivap) will become visible in the Set Injection Properties dialog box in Ansys Fluent. See Hooking DEFINE_DPM_HEAT_MASS UDFs for details on how to hook your DEFINE_DPM_HEAT_MASS UDF to Ansys Fluent.

You can use DEFINE_DPM_INJECTION_INIT to initialize a particle’s injection properties such as location, diameter, and velocity.

DEFINE_DPM_INJECTION_INIT (name, I)

Function returns

void

There are two arguments to DEFINE_DPM_INJECTION_INIT: name and I. You supply name, the name of the UDF. I is a variable that is passed by the Ansys Fluent solver to your UDF.

The following UDF, named init_bubbles, initializes particles on a surface injection due to a surface reaction. This function must be executed as a compiled UDF. Note that if you are going to use this UDF in a transient simulation to compute transient particles, you will need to replace loop(p, I->p) with loop(p, I->p_init). Transient particle initialization cannot be performed with a loop over I->p.

/**********************************************************************
   UDF that initializes particles on a surface injection due
   to a surface reaction
 ***********************************************************************/
 
 #include "udf.h"
 #include "surf.h" /* RP_CELL and RP_THREAD are defined in surf.h */
 
 #define REACTING_SURFACE_ID 2
 #define MW_H2 2
 #define STOIC_H2 1
 
 /* ARRHENIUS CONSTANTS */
 #define PRE_EXP 1e+15
 #define ACTIVE 1e+08
 #define BETA 0.0
 
 real arrhenius_rate(real temp)
 {
   return PRE_EXP*pow(temp,BETA)*exp(-ACTIVE/(UNIVERSAL_GAS_CONSTANT*temp));
 }
 
 /* Species numbers. Must match order in Ansys Fluent dialog box */
 #define HF 0
 
 /* Reaction Exponents */
 #define HF_EXP 2.0
 
 
 /* Reaction Rate Routine used in UDF */
 
 real reaction_rate(cell_t c, Thread *cthread,real mw[],real yi[])
 
 /* Note that all arguments in the reaction_rate function call in your .c source file 
   MUST be on the same line or a compilation error will occur */ 
 {
   real concenHF = C_R(c,cthread)*yi[HF]/mw[HF];
   return arrhenius_rate(C_T(c,cthread))*pow(concenHF,HF_EXP);
 }
 
 real contact_area(cell_t c,Thread *t,int s_id,int *n);
 
 
 DEFINE_DPM_INJECTION_INIT(init_bubbles,I)
 {
   int count,i;
   real area, mw[MAX_SPE_EQNS], yi[MAX_SPE_EQNS];
   /* MAX_SPE_EQNS is an Ansys Fluent constant in materials.h */
   Particle *p;
   cell_t cell;
   Thread *cthread;
   Material *mix, *sp;
   Message("Initializing Injection: %s\n",I->name);
   loop(p,I->p) /* Standard Ansys Fluent Looping Macro to get particle
          streams in an Injection */
     {
        cell = PP_CELL(p);  /* Get the cell and thread that
                             * the particle is currently in */
        cthread = PP_CELL_THREAD(p);
        /* Set up molecular weight & mass fraction arrays */
        mix = THREAD_MATERIAL(cthread);
        mixture_species_loop(mix,sp,i)
            {
             mw[i] = MATERIAL_PROP(sp,PROP_mwi);
             yi[i] = C_YI(cell,cthread,i);
            }
       area = contact_area(cell, cthread, REACTING_SURFACE_ID,&count);
       /* Function that gets total area of REACTING_SURFACE faces in
         contact with cell */
       /* count is the number of contacting faces, and is needed
         to share the total bubble emission between the faces  */
        if (count > 0) /* if cell is in contact with REACTING_SURFACE */
          {
           PP_FLOW_RATE(p) = (area *MW_H2* STOIC_H2 *
              reaction_rate(cell, cthread, mw, yi))/
              (real)count; /* to get correct total flow
               rate when multiple faces contact the same cell */
           PP_DIAM(p) = 1e-3;
           PP_RHO(p) = 1.0;
           PP_MASS(p) = PP_RHO(p)*M_PI*pow(PP_DIAM(p),3.0)/6.0;
          }
        else
           PP_FLOW_RATE(p) = 0.0;
         }
     }
 
 real contact_area(cell_t c, Thread *t, int s_id, int *n)
 {
   int i = 0;
   real area = 0.0, A[ND_ND];
   *n = 0;
   c_face_loop(c,t,i)
     {
      if(THREAD_ID(C_FACE_THREAD(c,t,i)) == s_id)
          {
          (*n)++;
         F_AREA(A,C_FACE(c,t,i), C_FACE_THREAD(c,t,i));
         area += NV_MAG(A);
          }
     } return area;
 }

When you are using a DPM particle sampling file or a VOF-to-DPM lump conversion transcript file as an unsteady injection file in a separate simulation, you can add more injection property data to such a file in extra columns. To do this, you can, for example, use a DEFINE_DPM_OUTPUT UDF during either sampling or VOF-to-DPM lump conversion (see Sampling of Trajectories in the Fluent User's Guide and DEFINE_DPM_OUTPUT for details). You can then use a DEFINE_DPM_INJECTION_INIT UDF to read that additional information from the file’s extra columns and initialize the particles that Ansys Fluent generates from the same unsteady injection file.

In the context of parallel processing, different particles generated from the unsteady injection file may be located in different partitions. Therefore, the DEFINE_DPM_INJECTION_INIT UDF must first obtain the contents of the file on all compute-node processes and then, for each new particle found in any partition, identify the corresponding line in the unsteady injection file. Only then can the additional information from the file be used, for example, to initialize user-defined particle properties.

The easiest way to obtain the contents of the file on every compute-node process is to read the file independently on every process, as shown in the example below. Note that Ansys Fluent itself reads the file only on compute-node 0 and then uses inter-process communication to automatically distribute the standard column data across all compute-node processes. You could implement a similar approach in your UDF for reading and sending information from the extra columns to other compute-node processes, but it involves extra programming that will not be further discussed here.

The following example can be used as a template for many purposes. The injudf UDF first reads the particle id (after the colon near the end of each line in the injection file) and then uses this id to initialize the first DPM particle user real in the new particle. For a typical purpose, there would be additional columns in the file, and you would use some extra fscanf(...) statement in this UDF to read them from the file.

#define USE_MY_OWN_FILE _NT

#if USE_MY_OWN_FILE
# define USE_FLUENT_IO_API  0
# define USE_UDF_PAR_IO  0
#endif

#include "udf.h"
#include "dpm_injections.h"

#define  REL_RESOL_IN_INJ_FILE 0.0001
/* If, for example, the file format is...:
 *  (( 4.9821e+00  -1.0061e+00  4.2000e+00  -1.6938e-05 ...
 * then there are five significant digits -->
 * the "relative resolution" is 0.0001 = 1.e-4
 */

/* Following definition divides by EPSILON to cancel
 * that from the definition of REAL_EQUAL_SCALED:
 */
#define MY_REAL_EQUAL(a, b)  \
   REAL_EQUAL_SCALED((a), (b), (REL_RESOL_IN_INJ_FILE / EPSILON))

DEFINE_DPM_INJECTION_INIT(injudf, I)
{
  Particle *p = NULL;
  char rest[1024] = "";
  int n_times_rewound = 0;
  FILE *my_inj_file_ptr = I->inj_fil_ptr;

  if ( ! I_DO_DPM)
    return;

  if (dpm_par.n_user_reals < 1)
  {
    DELAY_ERROR("node-%d, UDF '%s', Injection '%s', file '%s':\n"
                "Need at least one DPM user-real.",
                myid, __FUNCTION__, I->name, I->injection_file);
    return;
  }

#if USE_MY_OWN_FILE
  my_inj_file_ptr = fopen(I->injection_file, "r");
  /* Cannot use the open file pointer provided
   * by the Fluent executable.
   */
#endif

  if (NULLP(my_inj_file_ptr))
  {
    DELAY_ERROR("node-%d, UDF '%s', Injection '%s', file '%s':\n"
                "Cannot read from the injection file.",
                myid, __FUNCTION__, I->name, I->injection_file);
    return;
  }

#if USE_MY_OWN_FILE
  /* "read away" the file header: */
  fgets(rest, 1023, my_inj_file_ptr);
  fgets(rest, 1023, my_inj_file_ptr);
#elif 00
  rewind(I->inj_fil_ptr);
#endif

  int rewind_counter = 0;
  DPM_FPOS inj_file_start_pos = DPM_FTELL(my_inj_file_ptr);

  loop (p, I->p_init)
  {
    while (TRUE)
    {
      real file_x, file_y, file_z,
           file_u, file_v, file_w,
           file_d, file_T, file_M,
           file_m, file_n, file_t, file_f = 0.;

      int one = fgetc(my_inj_file_ptr);
      int two = fgetc(my_inj_file_ptr);

      if (feof(my_inj_file_ptr))
      {
        if (2 < rewind_counter)
          break;
        else
        {
          DPM_FSEEK(my_inj_file_ptr, inj_file_start_pos, SEEK_SET);
          ++rewind_counter;
          continue;   /* start over for this particle... */
        }
      }

      if ((one != '(') ||
          (two != '('))
      {
        DELAY_ERROR("node-%d, UDF '%s', Injection '%s', file '%s':\n"
                    "Cannot read from the file: '((' not found.",
                    myid, __FUNCTION__, I->name, I->injection_file);
        rewind_counter = 3;
        break;
      }

      int nmatch = fscanf(my_inj_file_ptr,
                      "%lg%lg%lg%lg%lg%lg%lg%lg%lg%lg%lg%lg%lg",
                      &file_x, &file_y, &file_z,
                      &file_u, &file_v, &file_w,
                      &file_d, &file_T, &file_M,
                      &file_m, &file_n, &file_t, &file_f);

      if (13 != nmatch)
      {
        DELAY_ERROR("node-%d, UDF '%s', Injection '%s', file '%s':\n"
                    "Cannot read at least 13 values from a line.",
                    myid, __FUNCTION__, I->name, I->injection_file);
        rewind_counter = 3;
        break;
      }

      /* Read additional numerical columns here... */

      /* Read the rest of the line into a string buffer: */
      fgets(rest, 1023, my_inj_file_ptr);

      /* if ( ! Check_Inj_File_Type(I, FALSE)) --- may change the file read pointer..! */
      /* The following block is relevant for *unsteady* injection files only...(!!)..: */
      {
        file_f -= I->start_flow_time_in_unst_file;

        if (file_f < 0.)   /* Before the start time, */
        {
          char dummyline[1024] = "";
          /* advance to end of line */
          do ;  /* nothing */
          while (( ! NULLP(fgets(dummyline, 1024, I->inj_fil_ptr))) &&
                 *dummyline &&   /* At least one character read? */
                 NULLP(strchr(dummyline, '\n'))); /* no LF read? */

          continue;   /* go on searching the file for the right entry */
        }

        if (0. < I->repeat_intrvl_from_unst_file)
        {
          if (file_f > I->repeat_intrvl_from_unst_file)
          {
            if (2 < ++n_times_rewound)
            {
              DELAY_ERROR("node-%d, UDF '%s', Injection '%s', file '%s':\n"
                          "Cannot find all expected entries in the file.",
                          myid, __FUNCTION__, I->name, I->injection_file);
              rewind_counter = 3;
              break;
            }

#if USE_MY_OWN_FILE
            rewind(my_inj_file_ptr);
            /* "read away" the file header: */
            fgets(rest, 1023, my_inj_file_ptr);
            fgets(rest, 1023, my_inj_file_ptr);
#else
            DPM_FSEEK(I->inj_fil_ptr, I->inj_fil_loc_bfr_rpt_intvl, SEEK_SET);
#endif
            file_f = 0.;
            continue;   /* go on searching the file for the right entry */
          }
        }

        file_f += I->unsteady_start;

        if (0. < I->repeat_intrvl_from_unst_file)   /* periodically repeating? */
        {
          /* file_f += I->repeat_intrvl_from_unst_file * (real) I->inj_fil_repeat_count; */
          while (file_f < dpm_par.time)
            file_f += I->repeat_intrvl_from_unst_file;
        }

        if (PP_INJ_FLOWTIME(p) > (dpm_par.time + dpm_par.time_step))
        {
          ASSERT(NULLP(p->next));   /* last particle -- the one that has been "read ahead" */
          break;  /* done for this particle from I->p_init, break out of loop through file */
        }
      }

      if (MY_REAL_EQUAL(PP_POS(p)[0],         file_x) &&
          MY_REAL_EQUAL(PP_POS(p)[1],         file_y) &&
          MY_REAL_EQUAL(PP_POS(p)[2],         file_z) &&
          MY_REAL_EQUAL(PP_VEL(p)[0],         file_u) &&
          MY_REAL_EQUAL(PP_VEL(p)[1],         file_v) &&
          MY_REAL_EQUAL(PP_VEL(p)[2],         file_w) &&
          MY_REAL_EQUAL(PP_DIAM(p),           file_d) &&
          MY_REAL_EQUAL(PP_T(p),              file_T) &&
          MY_REAL_EQUAL(PP_FLOW_RATE(p),      file_M) &&
       /* The following PP_...(p) values aren't yet known in
        * parallel and will be calculated from P_FLOW_RATE(p)
        * later:
        * MY_REAL_EQUAL(PP_MASS(p) * PP_N(p), file_M) &&
        * MY_REAL_EQUAL(PP_MASS(p),           file_m) &&
        * MY_REAL_EQUAL(PP_N(p),              file_n) &&
        */
          MY_REAL_EQUAL(PP_INJ_DURATIME(p),   file_t) &&
          MY_REAL_EQUAL(PP_INJ_FLOWTIME(p),   file_f))
      {
        /* From the rest of the line, determine some old part_id,
         * add 7000 and assign to the new particle's user real:
         */
        sscanf(strchr(rest, ':') + 1, "%lg", &(PP_USER_REAL(p, 0)));
        PP_USER_REAL(p, 0) += 7000.;
	/* Message("Value of scalar 0 \t %lg\n",PP_USER_REAL(p,0)); */
        break;  /* done for this particle from I->p_init, break out of loop through file */
      }
    }   /* while(TRUE) */

    /* End of file, but yet another particle to do? */
    if (2 < rewind_counter &&  ! NULLP(p))
    {
      DELAY_ERROR("node-%d, UDF '%s', Injection '%s', file '%s':\n"
                  "Could not process all new particles.",
                  myid, __FUNCTION__, I->name, I->injection_file);
      break;
    }
  }

#if USE_MY_OWN_FILE
  fclose(my_inj_file_ptr);
#endif
}

As mentioned above, in a parallel environment, the UDF must identify for every particle on each compute-node process which line from the injection file describes the particle. Therefore, the UDF reads all information from the file into every compute-node process and then compares the values from the standard columns against the properties of each particle in order to find the matching file entry for every particle. Once the correct line has been found, the UDF processes additional particle property information and assigns it to the particle.

The preprocessor macro _NT, which is used in the definition of the macro USE_MY_OWN_FILE, is defined to TRUE (1) by udf.h on Windows only and not on Linux. This is needed because on Windows, the UDF must read the file independently of the file handle already established by the Fluent executable. The reason for this is that, as long as you compile your UDF with the recommended compiler, it cannot use the file handle data structure provided by the Fluent executable.

After the UDF that you have defined using DEFINE_DPM_INJECTION_INIT is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Set Injection Properties dialog box in Ansys Fluent.

See Hooking DEFINE_DPM_INJECTION_INIT UDFs for details on how to hook your DEFINE_DPM_INJECTION_INIT UDF to Ansys Fluent.

You can use DEFINE_DPM_LAW to customize laws for particles. For example your UDF can specify custom laws for heat and mass transfer rates for droplets and combusting particles. Additionally, you can specify custom laws for mass, diameter, and temperature properties as the droplet or particle exchanges mass and energy with its surroundings.

When a DEFINE_DPM_LAW UDF is activated, then the number of species that can be referenced and interact with the particles via the UDF is limited to those with a species index less than the maximum UDF species number, defined using the TUI command define/models/dpm/options/maximum-udf-species. The default number for maximum-udf-species is 50.

DEFINE_DPM_LAW (name, tp, ci)

Function returns

void

There are three arguments to DEFINE_DPM_LAW: name, tp, and ci. You supply name, the name of the UDF. tp and ci are variables that are passed by the Ansys Fluent solver to your UDF.

The following UDF, named Evapor_Swelling_Law, models a custom law for the evaporation swelling of particles. The source code can be interpreted or compiled in Ansys Fluent. See Example DEFINE_DPM_LAW usage.

/**********************************************************************
   UDF that models a custom law for evaporation swelling of particles
 ***********************************************************************/
 
 #include "udf.h"
 
 DEFINE_DPM_LAW(Evapor_Swelling_Law,tp,ci)
 {
    real swelling_coeff = 1.1;

    /* first, call standard evaporation routine to calculate
    the mass and heat transfer       */
    VaporizationLaw(tp);

    /* compute new particle diameter and density */
    TP_DIAM(tp) = TP_INIT_DIAM(tp)*(1. + (swelling_coeff - 1.)*
   (TP_INIT_MASS(tp)-TP_MASS(tp))/(TP_VOLATILE_FRACTION(tp)*TP_INIT_MASS(tp)));
    TP_RHO(tp) = TP_MASS(tp) / (3.14159*TP_DIAM(tp)*TP_DIAM(tp)*TP_DIAM(tp)/6);
    TP_RHO(tp) = MAX(0.1, MIN(1e5, TP_RHO(tp)));
 } 

After the UDF that you have defined using DEFINE_DPM_LAW is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Custom Laws dialog box in Ansys Fluent. See Hooking DEFINE_DPM_LAW UDFs for details on how to hook your DEFINE_DPM_LAW UDF to Ansys Fluent.

You can use DEFINE_DPM_OUTPUT to modify what is written to the sampling device output. This function allows access to the variables that are written to a file as a particle passes through a sampler (see Sampling of Trajectories in the Fluent User's Guide for details).

DEFINE_DPM_OUTPUT (name, header, fp, tp, t, plane)

Function returns

void

There are six arguments to DEFINE_DPM_OUTPUT: name, header, fp, tp, t, and plane. You supply name, the name of the UDF. header, fp, tp, t, and plane are variables that are passed by the Ansys Fluent solver to your UDF. The output of your UDF will be written to the file indicated by fp.

When using DEFINE_DPM_OUTPUT to write sample files, certain special file operations must be performed by Ansys Fluent. Therefore, the usual C output function fprintf cannot be used. Instead, you must use the functions par_fprintf and par_fprintf_head. For details on the use of these functions, refer to The par_fprintf_head and par_fprintf Functions and the following example. In particular, note that the first two arguments passed to par_fprintf must always be given as indicated, because they are used to sort the contents of the file but they will not appear in the final output.

The following UDF named discrete_phase_sample samples the size and velocity of discrete phase particles at selected planes downstream of an injection. For 2D axisymmetric simulations, it is assumed that droplets/particles are being sampled at planes (lines) corresponding to constant . For 3D simulations, the sampling planes correspond to constant .

To remove particles from the domain after they have been sampled, change the value of REMOVE_PARTICLES to TRUE. In this case, particles will be deleted following the time step in which they cross the plane. This is useful when you want to sample a spray immediately in front of an injector and you do not want to track the particles further downstream.

#include "udf.h"

/*******************************************************************/
/* UDF that samples discrete phase size and velocity distributions */
/* within the domain.                                              */
/*******************************************************************/

#define REMOVE_PARTICLES FALSE

DEFINE_DPM_OUTPUT(discrete_phase_sample,header,fp,tp,t,plane)
{
#if RP_2D
  real y;

  if(header)
  {
    par_fprintf_head(fp," #Time[s]  R [m]  X-velocity[m/s]");
    par_fprintf_head(fp," W-velocity[m/s] R-velocity[m/s] ");
    par_fprintf_head(fp,"Drop Diameter[m] Number of Drops  ");
    par_fprintf_head(fp,"Temperature [K] Initial Diam [m] ");
    par_fprintf_head(fp,"Injection Time [s] \n");
  }

  if(NULLP(tp))
    return;

  if (rp_axi && (sg_swirl || GVAR_TURB(rp, ke)))
    y = MAX(sqrt(SQR(TP_POS(tp)[1]) + SQR(TP_POS(tp)[2])),DPM_SMALL);
  else
    y = TP_POS(tp)[1];

  /* Note: The first two arguments to par_fprintf are used internally and */
  /* must not be changed, even though they do not appear in the final output.*/
  par_fprintf(fp, "%d %" int64_fmt " %e %f %f %f %f %e %e %f %e %f \n",
              P_INJ_ID(TP_INJECTION(tp)),
              TP_ID(tp),
              TP_TIME(tp),
              y,
              TP_VEL(tp)[0],
              TP_VEL(tp)[1],
              TP_VEL(tp)[2],
              TP_DIAM(tp),
              TP_N(tp),
              TP_T(tp),
              TP_INIT_DIAM(tp),
              TP_TIME_OF_BIRTH(tp));
#else
  real r, x, y;

  if(header)
  {
    par_fprintf_head(fp," #Time[s] R [m] x-velocity[m/s] ");
    par_fprintf_head(fp,"y-velocity[m/s] z-velocity[m/s]  ");
    par_fprintf_head(fp,"Drop Diameter[m]  Number of Drops ");
    par_fprintf_head(fp,"Temperature [K]  Initial Diam [m] ");
    par_fprintf_head(fp,"Injection Time [s] \n");
  }

  if(NULLP(tp))
    return;

  x = TP_POS(tp)[0];
  y = TP_POS(tp)[1];
  r = sqrt(SQR(x) + SQR(y));

  /* Note: The first two arguments to par_fprintf are used internally and */
  /* must not be changed, even though they do not appear in the final output.*/
  par_fprintf(fp,"%d %" int64_fmt " %e %f %f %f %f %e %e %f %e %f \n",
              P_INJ_ID(TP_INJECTION(tp)),
              TP_ID(tp),
              TP_TIME(tp),
              r,
              TP_VEL(tp)[0],
              TP_VEL(tp)[1],
              TP_VEL(tp)[2],
              TP_DIAM(tp),
              TP_N(tp),
              TP_T(tp),
              TP_INIT_DIAM(tp),
              TP_TIME_OF_BIRTH(tp));
#endif

#if REMOVE_PARTICLES
  MARK_TP(tp, P_FL_REMOVED);
#endif
}

The following example provides the source code used in Ansys Fluent simulations when you do not use a DEFINE_DPM_OUTPUT UDF. You can modify this template to adapt it to your needs.

#include "udf.h"

/*****************************************************************/
/* DPM sampling output UDF that does what Fluent does by default */
/*****************************************************************/

#define REMOVE_PARTICLES FALSE

DEFINE_DPM_OUTPUT(my_dpm_out, header, fp, tp, thread, plane)
{
  if (header)
  {
    char *sort_name;
    char sort_fn[4096];

    if (NNULLP(thread))
      sort_name = THREAD_HEAD(thread)->dpm_summary.sort_file_name;
    else if ( ! NULLP(plane))
      sort_name = plane->sort_file_name;
    else   /* This is not expected to happen for regular particle sampling.. */
    {
      if (dpm_par.unsteady_tracking)
        sort_name = "parcels";
      else
        sort_name = "tracks";
    }

    /* sort_name may contain "/" (Linux)
     * or ":" and "\" (Windows) --
     * replace them all by "_":
     */
    strcpy(sort_fn, sort_name);
    replace_path_chars_in_string(sort_fn);
    if (dpm_par.unsteady_tracking)
      par_fprintf_head(fp, "(%s %d)\n", sort_fn, 13);
    else
      par_fprintf_head(fp, "(%s %d)\n", sort_fn, 12);

#if RP_2D
    if (rp_axi_swirl)
      par_fprintf_head(fp, "(x    r    theta    u    v    w");
    else
#endif
      par_fprintf_head(fp, "(x    y    z    u    v    w");

    if (dpm_par.unsteady_tracking)
      par_fprintf_head(fp, "  diameter    t    parcel-mass  "
                           "  mass    n-in-parcel   time   flow-time   name)\n");
    else
      par_fprintf_head(fp, "  diameter    t    mass-flow  "
                           "  mass    frequency    time    name)\n");
  }
  else if ( ! NULLP(tp))
  {
    /* Do some preparatory calculations for later use:
     */
    real flow_rate = 0.;
    real V_vel = TP_VEL(tp)[1];
    real W_vel = TP_VEL(tp)[2];
    real Y = TP_POS(tp)[1];
    real Z = TP_POS(tp)[2];
    real strength = 0.;
    real mass = 0.;

    if (TP_INJECTION(tp)->type != DPM_TYPE_MASSLESS)
    {
      mass = TP_MASS(tp);

      if (dpm_par.unsteady_tracking)
        strength = TP_N(tp);
      else
      {
        strength = TP_FLOW_RATE(tp) / TP_INIT_MASS(tp);
        if (TP_STOCHASTIC(tp))
          strength /= (real)TP_STOCHASTIC_NTRIES(tp);
      }

      flow_rate = strength * mass;
    }

#if RP_2D
    if (rp_axi_swirl)
    {
      Y = MAX(sqrt(TP_POS(tp)[1] * TP_POS(tp)[1] + TP_POS(tp)[2] * TP_POS(tp)[2]), DPM_SMALL);
      V_vel = (TP_VEL(tp)[1] * TP_POS(tp)[1] + TP_VEL(tp)[2] * TP_POS(tp)[2]) / Y;
      W_vel = (TP_VEL(tp)[2] * TP_POS(tp)[1] - TP_VEL(tp)[1] * TP_POS(tp)[2]) / Y;
      if (Y > 1.e-20) Z = LIMIT_ACOS(TP_POS(tp)[1] / Y);
    }
#endif

    if ( ! dpm_par.unsteady_tracking)
      par_fprintf(fp, /* Note: The first two arguments to par_fprintf are */
                      /* used internally and must not be changed, even */
                      /* though they do not appear in the final output.*/
                  "%d %" int64_fmt " ((%e  %e  %e  %e  %e  %e "
                  " %e  %e  %e  %e  %e  %e) %s:%" int64_fmt ")\n",
                  P_INJ_ID(TP_INJECTION(tp)), tp->part_id,
                  TP_POS(tp)[0],
                  Y,
                  Z,
                  TP_VEL(tp)[0],
                  V_vel,
                  W_vel,
                  TP_DIAM(tp),
                  TP_T(tp),
                  flow_rate,
                  mass,
                  strength,
                  TP_TIME(tp) - tp->time_of_birth,
                  TP_INJECTION(tp)->name,
                  tp->part_id);
    else
      par_fprintf(fp, /* Note: The first two arguments to par_fprintf are */
                      /* used internally and must not be changed, even */
                      /* though they do not appear in the final output.*/
                  "%d %" int64_fmt " ((%e  %e  %e  %e  %e  %e "
                  " %e  %e  %e  %e  %e  %e  %e) %s:%" int64_fmt ")\n",
                  (int) TP_TIME(tp),
                  (int64_t) ((TP_TIME(tp) -
                              (real) ((int) TP_TIME(tp))) *
                             (real) 1000000000000.),
                  TP_POS(tp)[0],
                  Y,
                  Z,
                  TP_VEL(tp)[0],
                  V_vel,
                  W_vel,
                  TP_DIAM(tp),
                  TP_T(tp),
                  flow_rate,
                  mass,
                  strength,
                  TP_TIME(tp) - tp->time_of_birth,
                  TP_TIME(tp),
                  TP_INJECTION(tp)->name,
                  tp->part_id);

#if REMOVE_PARTICLES
    MARK_TP(tp, P_FL_REMOVED);
#endif
  }
}

For a VOF-to-DPM model transition simulation, if you use the option to write information about every liquid lump converted into DPM particle parcels to a file (as described in Setting up the VOF-to-DPM Model Transition in the Fluent User's Guide), you can use a DEFINE_DPM_OUTPUT function to define the file content. After you hook the DEFINE_DPM_OUTPUT function as described in Hooking a DPM Output UDF to Ansys Fluent, that function will be also called multiple times every time the VOF-to-DPM model transition mechanism is triggered.

You can also use the DEFINE_DPM_OUTPUT function to control whether a liquid droplet is represented by a single particle parcel or split into multiple particle parcels as follows:

#include "udf.h"
#include "vof_to_dpm.h"

/*****************************************************************/
/* DPM sampling output UDF that does what Fluent does by default */
/*****************************************************************/

DEFINE_DPM_OUTPUT(my_dpm_out, header, fp, tp, thread, plane)
{
  if (header)
  {
    char *sort_name;
    char sort_fn[4096];

    if (NNULLP(thread))
      sort_name = THREAD_HEAD(thread)->dpm_summary.sort_file_name;
    else if ( ! NULLP(plane))
      sort_name = plane->sort_file_name;
    else   /* This is not expected to happen for regular particle sampling.. */
    {
      if ( ! NULLP(convert_lump_args.myldps.injection))
       sort_name = convert_lump_args.myldps.injection->name;
      else
      if (dpm_par.unsteady_tracking)
        sort_name = "parcels";
      else
        sort_name = "tracks";
    }

    /* sort_name may contain "/" (Linux)
     * or ":" and "\" (Windows) --
     * replace them all by "_":
     */
    strcpy(sort_fn, sort_name);
    replace_path_chars_in_string(sort_fn);
    if (dpm_par.unsteady_tracking)
      par_fprintf_head(fp, "(%s %d)\n", sort_fn, 13);
    else
      par_fprintf_head(fp, "(%s %d)\n", sort_fn, 12);

#if RP_2D
    if (rp_axi_swirl)
      par_fprintf_head(fp, "(x    r    theta    u    v    w");
    else
#endif
      par_fprintf_head(fp, "(x    y    z    u    v    w");

    if (dpm_par.unsteady_tracking)
      par_fprintf_head(fp, "  diameter    t    parcel-mass  "
                           "  mass    n-in-parcel   time   flow-time  name)\n");
    else
      par_fprintf_head(fp, "  diameter    t    mass-flow  "
                           "  mass    frequency    time    name)\n");
  }
  else if ( ! NULLP(tp))
  {
    /* Do some preparatory calculations for later use:
     */
    real flow_rate = 0.;
    real V_vel = TP_VEL(tp)[1];
    real W_vel = TP_VEL(tp)[2];
    real Y = TP_POS(tp)[1];
    real Z = TP_POS(tp)[2];
    real strength = 0.;
    real mass = 0.;

    if (TP_INJECTION(tp)->type != DPM_TYPE_MASSLESS)
    {
      mass = TP_MASS(tp);

      if (dpm_par.unsteady_tracking)
        strength = TP_N(tp);
      else
      {
        strength = TP_FLOW_RATE(tp) / TP_INIT_MASS(tp);
        if (TP_STOCHASTIC(tp))
          strength /= (real)TP_STOCHASTIC_NTRIES(tp);
      }

      flow_rate = strength * mass;
    }

#if RP_2D
    if (rp_axi_swirl)
    {
      Y = MAX(sqrt(TP_POS(tp)[1] * TP_POS(tp)[1] + TP_POS(tp)[2] * TP_POS(tp)[2]), DPM_SMALL);
      V_vel = (TP_VEL(tp)[1] * TP_POS(tp)[1] + TP_VEL(tp)[2] * TP_POS(tp)[2]) / Y;
      W_vel = (TP_VEL(tp)[2] * TP_POS(tp)[1] - TP_VEL(tp)[1] * TP_POS(tp)[2]) / Y;
      if (Y > 1.e-20) Z = LIMIT_ACOS(TP_POS(tp)[1] / Y);
    }
#endif

    if ( ! NULLP(convert_lump_args.domain))
    {
      /* This is definitely VOF-to-DPM calling,
       * so demonstrate some things we can do:
       */
      if (TP_POS(tp)[0] > 0.07)   /* Silly condition for demonstration purposes. */
      {
        TP_N(tp) = -1.;           /* Tell the code not even to convert the lump. */
        return;                   /* Do not go on to write anything to the file. */
      }
      else
      if (TP_POS(tp)[0] > 0.04)   /* Silly condition for demonstration purposes. */
      {
        MARK_TP(tp, P_FL_REMOVED);    /* Drop the particle, never track it. */
        return;              /* Do not go on to write anything to the file. */
      }
      else
      {
        /* Just for demonstration that it works,
         * mirror the particle at the XY plane,
         * or set its temperature to 777:
         */
        /* TP_POS(tp)[2] *= -1.; */
        TP_T(tp) = 777.;
      }
    }

    if (NULLP(convert_lump_args.domain) ||  /* For future use: Not VOF-to-DPM, but...: */
        ( ! NULLP(thread)) ||               /* ...thread AND plane may still be NULL... */
        ( ! NULLP(plane))  ||  /* When VOF-to-DPM calls this, "thread" and "plane" are NULL */
#if RP_NODE
        I_AM_NODE_ZERO_P   ||  /* and only ONE compute-node process must write to the file: */
#endif
        I_AM_NODE_HOST_P   ||  /* (This line for shared-memory parallel DPM tracking only) */
        I_AM_NODE_SERIAL_P)    /* (This line for SERIAL ["-t0"] execution [unsupported]) */
      if ((1. == strength) ||  /* VOF-to-DPM: want to record one single entry per lump! */
          (NULLP(convert_lump_args.domain)))   /* Otherwise, write every particle.. */
      {
        if ( ! dpm_par.unsteady_tracking)
          par_fprintf(fp, /* The first two arguments to par_fprintf are used internally and */
                          /* must not be changed, although they do not appear in the output.*/
                      "%d %" int64_fmt " ((%.6e  %.6e  %.6e "
                      " %.6e  %.6e  %.6e  %.6e  %.6e  %.6e "
                      " %.6e  %.6e  %.6e) %s:%" int64_fmt ")\n",
                      P_INJ_ID(TP_INJECTION(tp)), tp->part_id,
                      TP_POS(tp)[0],
                      Y,
                      Z,
                      TP_VEL(tp)[0],
                      V_vel,
                      W_vel,
                      TP_DIAM(tp),
                      TP_T(tp),
                      flow_rate,
                      mass,
                      strength,
                      TP_TIME(tp) - tp->time_of_birth,
                      TP_INJECTION(tp)->name,
                      tp->part_id);
        else
          par_fprintf(fp, /* The first two arguments to par_fprintf are used internally and */
                          /* must not be changed, although they do not appear in the output.*/
                      "%d %" int64_fmt " ((%.6e  %.6e  %.6e "
                      " %.6e  %.6e  %.6e  %.6e  %.6e  %.6e "
                      " %.6e  %.6e  %.6e  %.6e) %s:%" int64_fmt ")\n",
                      (int) TP_TIME(tp),
                      (int64_t) ((TP_TIME(tp) -
                                  (real) ((int) TP_TIME(tp))) *
                                 (real) 1000000000000.),
                      TP_POS(tp)[0],
                      Y,
                      Z,
                      TP_VEL(tp)[0],
                      V_vel,
                      W_vel,
                      TP_DIAM(tp),
                      TP_T(tp),
                      flow_rate,
                      mass,
                      strength,
                      TP_TIME(tp) - tp->time_of_birth,
                      TP_TIME(tp),
                      TP_INJECTION(tp)->name,
                      tp->part_id);
      }
  }
  else
  {
    /* tp is NULL: VOF-to-DPM calls us TWICE synchronously on all compute-node proc.
     * so that we have an opportunity e.g. to do our own lump characterisation (i.e.
     * calculate some lump properties that have not been calculated by Fluent yet),
     * or can interfere with the lump properties Fluent has determined so far:
     */
    Message0("\nmy_dpm_out: last_lump_id:  %d  =="
  "==  file yet opened? -- %s  ====\n",
          convert_lump_args.last_lump_id, NULLP(fp) ? "NO" : "yes!");
  }
}

After the UDF that you have defined using DEFINE_DPM_OUTPUT is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Sample Trajectories dialog box in Ansys Fluent. See Hooking DEFINE_DPM_OUTPUT UDFs for details on how to hook your DEFINE_DPM_OUTPUT UDF to Ansys Fluent.

You can use DEFINE_DPM_PROPERTY to specify properties of discrete phase materials. You can model the following dispersed phase properties with this type of UDF:

DEFINE_DPM_PROPERTY (name, c, t, tp, T)

Function returns

real

There are five arguments to DEFINE_DPM_PROPERTY: name, c, t, tp, and T. You supply name, the name of the UDF. c, t, tp, and T are variables that are passed by the Ansys Fluent solver to your UDF. Your UDF will need to compute the real value of the discrete phase property and return it to the solver.

If you are using DEFINE_DPM_PROPERTY to specify the specific heat for particle materials, your UDF will also need to set the value of the particle enthalpy in the Tracked_Particle *tp, tp->enthalpy, to the particle sensible enthalpy, which should be calculated as the temperature integral of the specific heat function from the reference temperature, T_REF, to the temperature, T.

In the following example, three discrete phase material property UDFs (named coal_emissivity , coal_scattering, and coal_cp, respectively) are concatenated into a single C source file. These UDFs must be executed as compiled UDFs in Ansys Fluent.

/*********************************************************************
 *    UDF that specifies discrete phase material properties
 *********************************************************************
 */

#include "udf.h"

DEFINE_DPM_PROPERTY(coal_emissivity, c, t, tp, T)
{
  real mp0;
  real mp;
  real vf, cf;

  if (NULLP(tp->pp) ||
      NULLP(TP_CELL_THREAD(tp)))
    return 1.0;                    /* initial value */

  mp0 = TP_INIT_MASS(tp);
  mp  = TP_MASS(tp);

  /* get the material char and volatile fractions and
   * store them in vf and cf:
   */
  vf = TP_VOLATILE_FRACTION(tp);
  cf = TP_CHAR_FRACTION(tp);

  if (!(((mp / mp0) >= 1) ||
        ((mp / mp0) <= 0)))
    {
      if ((mp / mp0) < (1 - vf - cf))
        {
          /* only ash left  */
          /* vf = cf = 0    */

          return .001;
        }
      else if ((mp / mp0) < (1 - vf))
        {
          /* only ash and char left               */
          /* cf = 1 - (1 - vf - cf) / (mp / mp0)  */
          /* vf = 0                               */

          return 1.0;
        }
      else
        {
          /* volatiles, char, and ash left   */
          /* cf =          cf  / (mp / mp0)  */
          /* vf = 1 - (1 - vf) / (mp / mp0)  */

          return 1.0;
        }
    }

  return 1.0;
}


DEFINE_DPM_PROPERTY(coal_scattering, c, t, tp, T)
{
  real mp0;
  real mp;
  real cf, vf;

  if (NULLP(tp->pp) ||
      NULLP(TP_CELL_THREAD(tp)))
    return 1.0;                    /* initial value */


  mp0 = TP_INIT_MASS(tp);
  mp  = TP_MASS(tp);

  /* get the original char and volatile fractions
  /* and store them in vf and cf:
   */
  vf = TP_VOLATILE_FRACTION(tp);
  cf = TP_CHAR_FRACTION(tp);

  if (!(((mp / mp0) >= 1) ||
        ((mp / mp0) <= 0)))
    {
      if ((mp / mp0) < (1 - vf - cf))
        {
          /* only ash left  */
          /* vf = cf = 0    */

          return 1.1;
        }
      else if ((mp / mp0) < (1 - vf))
        {
          /* only ash and char left               */
          /* cf = 1 - (1 - vf - cf) / (mp / mp0)  */
          /* vf = 0                               */

          return 0.9;
        }
      else
        {
          /* volatiles, char, and ash left   */
          /* cf =          cf  / (mp / mp0)  */
          /* vf = 1 - (1 - vf) / (mp / mp0)  */

          return 1.0;
        }
    }

  return 1.0;
}


DEFINE_DPM_PROPERTY(coal_cp, c, t, tp, T)
{
  real mp0;
  real mp;
  real cf;
  real vf;
  real af;
  real Cp;

  if (NULLP(tp->pp) ||
      NULLP(TP_CELL_THREAD(tp)))
    {
      Cp = 1600.;                       /* initial value */
    }
  else
    {
      mp0 = TP_INIT_MASS(tp);
      mp  = TP_MASS(tp);

      cf  = TP_CF(tp);                  /* char fraction */
      vf  = TP_VF(tp);                  /* volatiles fraction */
      af  = 1. - TP_VF(tp) - TP_CF(tp); /* ash fraction */

      Cp  = 2000. * af +
            1100. * vf +
            1300. * cf;
    }

  tp->enthalpy = Cp * (T - T_REF);

  return Cp;
}

After the UDF that you have defined using DEFINE_DPM_PROPERTY is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Create/Edit Materials dialog box in Ansys Fluent. See Hooking DEFINE_DPM_PROPERTY UDFs for details on how to hook your DEFINE_DPM_PROPERTY UDF to Ansys Fluent.

You can use DEFINE_DPM_SCALAR_UPDATE to update scalar quantities every time a particle position is updated. The function allows particle-related variables to be updated or integrated over the life of the particle. Particle values can be stored in an array associated with the Tracked_Particle (accessed with the macro TP_USER_REAL(tp,i)). Values calculated and stored in the array can be used to color the particle trajectory.

During Ansys Fluent execution, the DEFINE_DPM_SCALAR_UPDATE function is called at the start of particle integration and then after each time step for the particle trajectory integration. A value of 1 for initialize will be passed to the UDF when the particle is first injected or each time the particle tracker is called for unsteady particle tracking.

DEFINE_DPM_SCALAR_UPDATE (name, c, t, initialize, tp)

Function returns

void

There are five arguments to DEFINE_DPM_SCALAR_UPDATE: name, c, t, initialize, and tp. You supply name, the name of the UDF. c, t, initialize, and tp are variables that are passed by the Ansys Fluent solver to your UDF.

The following compiled UDF computes the melting index along a particle trajectory. The DEFINE_DPM_SCALAR_UPDATE function is called at every particle time step in Ansys Fluent and requires a significant amount of CPU time to execute.

The melting index is computed from

(2–20)

Also included in this UDF is an initialization function DEFINE_INIT that is used to initialize the scalar variables. DPM_OUTPUT is used to write the melting index at sample planes and surfaces. The macro NULLP(p), which expands to ((p) == NULL), checks if its argument p is a null pointer.

/*********************************************************************
   UDF for computing the melting index along a particle trajectory
 **********************************************************************/
 #include "udf.h"
 
 DEFINE_INIT(melt_setup,domain)
 {
    /* if memory for the particle variable titles has not been
   * allocated yet, do it now */
    if (NULLP(user_particle_vars)) Init_User_Particle_Vars();
      /* now set the name and label */
      strcpy(user_particle_vars[0].name,"melting-index");
      strcpy(user_particle_vars[0].label,"Melting Index");
      strcpy(user_particle_vars[1].name,"melting-index-0");
      strcpy(user_particle_vars[1].label,"Melting Index 0");
 }
 
 /* update the user scalar variables */
 DEFINE_DPM_SCALAR_UPDATE(melting_index,cell,thread,initialize,tp)
 {
    cphase_state_t *c = &(tp->cphase[0]);
    if (initialize)
      {
     /* this is the initialization call, set:
     * TP_USER_REAL(tp,0) contains the melting index, initialize to 0
     * TP_USER_REAL(tp,1) contains the viscosity at the start of a time step*/
     TP_USER_REAL(tp,0) = 0.;
     TP_USER_REAL(tp,1) = c->mu;
      }
    else
      {
       /* use a trapezoidal rule to integrate the melting index */
       TP_USER_REAL(tp,0) += TP_DT(tp) * .5 * (1/TP_USER_REAL(tp,1) + 1/c->mu);
       /* save current fluid viscosity for start of next step */   
      TP_USER_REAL(tp,1) = c->mu;
      }
   }
 
 /* write melting index when sorting particles at surfaces */
 DEFINE_DPM_OUTPUT(melting_output,header,fp,tp,thread,plane)
 {
    char name[100];
    if (header)
    {
     if (NNULLP(thread))
       par_fprintf_head(fp,"(%s %d)\n",THREAD_HEAD(thread)->
         dpm_summary.sort_file_name,11);
     else
       par_fprintf_head(fp,"(%s %d)\n",plane->sort_file_name,11);
       par_fprintf_head(fp,"(%10s %10s %10s %10s %10s %10s %10s"
        " %10s %10s %10s %10s %s)\n",
        "X","Y","Z","U","V","W","diameter","T","mass-flow",
        "time","melt-index","name");
      }
    else
      {
       sprintf(name,"%s:%d",TP_INJECTION(tp)->name,TP_ID(tp));
       /* Note: The first two arguments to par_fprintf are used internally and     */
       /* must not be changed, even though they do not appear in the final output. */
       par_fprintf(fp,
          "%d %d ((%10.6g %10.6g %10.6g %10.6g %10.6g %10.6g "
          "%10.6g %10.6g %10.6g %10.6g %10.6g) %s)\n",
        P_INJ_ID(TP_INJECTION(tp)), TP_ID(tp),
        TP_POS(tp)[0], TP_POS(tp)[1], TP_POS(tp)[2],
        TP_VEL(tp)[0], TP_VEL(tp)[1], TP_VEL(tp)[2],
        TP_DIAM(tp), TP_T(tp), TP_FLOW_RATE(tp), TP_TIME(tp),
        TP_USER_REAL(tp,0), name);
      }
 } 

After the UDF that you have defined using DEFINE_DPM_SCALAR_UPDATE is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Discrete Phase Model dialog box in Ansys Fluent.

See Hooking DEFINE_DPM_SCALAR_UPDATE UDFs for details on how to hook your DEFINE_DPM_SCALAR_UPDATE UDF to Ansys Fluent.

You can use DEFINE_DPM_SOURCE to specify particle source terms. The function allows access to the accumulated source terms for a particle in a given cell before they are added to the mass, momentum, and energy exchange terms for coupled DPM calculations.

When a DEFINE_DPM_SOURCE UDF is activated, then the number of species that can be referenced and interact with the particles via the UDF is limited to those with a species index less than the maximum UDF species number, defined using the TUI command define/models/dpm/options/maximum-udf-species. The default number for maximum-udf-species is 50.

DEFINE_DPM_SOURCE (name, c, t, S, strength, tp)

Function returns

void

There are six arguments to DEFINE_DPM_SOURCE: name, c, t, S, strength, and tp. You supply name, the name of the UDF. c, t, S, strength, and tp are variables that are passed by the Ansys Fluent solver to your UDF. The modified source terms, after they have been computed by the function, will be stored in S.

For an example of DEFINE_DPM_SOURCE usage with the species model, see Example.

For using DEFINE_DPM_SOURCE with the non-premixed model, see the example below.

The DEFINE_DPM_LAW named nonpremixed_example and the DEFINE_DPM_SOURCE named dpm_source provide a simplistic model for droplet vaporization for education purposes only. The DEFINE_DPM_LAW UDF also demonstrates how to define the additional terms required to include source term linearization in the DPM UDFs.

#include "udf.h"
DEFINE_DPM_LAW(nonpremixed_example, tp, coupled)
{
  real vapor_conc_bulk, vapor_conc_surf;
  real mp_dot = 0., remaining_mass=0.;

  const cphase_state_t *c = &(tp->cphase[0]);
  Material *m = TP_MATERIAL(tp);
  real area_old = DPM_SURFACE_AREA(TP_DIAM(tp));
  real Tp_old   = TP_T(tp);
  real dt = TP_DT(tp);

  real diffusion_coeff = DPM_BINARY_DIFFUSIVITY(tp, m, TP_T(tp));
  real Pr = (c->sHeat)*(c->mu)/(c->tCond);
  real Sc = (c->mu)/((c->rho)*(diffusion_coeff));
  real Nu = 2. + 0.6*sqrt(tp->Re)*pow(Pr,1./3.);
  real Sh = 2. + 0.6*sqrt(tp->Re)*pow(Sc,1./3.);
  real mass_transfer_coeff = Sh * diffusion_coeff / TP_DIAM(tp);

  real h_lat = DPM_LATENT_HEAT(tp, m);
  real p_saturation = DPM_VAPOR_PRESSURE(tp, m, TP_T(tp));

  vapor_conc_bulk = c->pressure / UNIVERSAL_GAS_CONSTANT
    / c->temp * c->yi[TP_EVAP_SPECIES_INDEX(tp)] * c->mol_wt_bulk
    / solver_par.molWeight[TP_EVAP_SPECIES_INDEX(tp)];
 
  if (p_saturation > c->pressure)
    p_saturation = c->pressure;
  else if (p_saturation < 0.0)
    p_saturation = 0.;
 
  vapor_conc_surf = p_saturation / UNIVERSAL_GAS_CONSTANT / TP_T(tp);
 
  mp_dot = mass_transfer_coeff * solver_par.molWeight[TP_EVAP_SPECIES_INDEX(tp)]
    * area_old * (vapor_conc_surf - vapor_conc_bulk);
 
  if (mp_dot < 0.)
    mp_dot = 0.0;
 
/* additional terms required for using DPM source term linearization */

  if (dpm_par.linearized_source_terms)
    tp->htcA = Nu * c->tCond * TP_DIAM(tp)* M_PI;
 
  if (dpm_par.linearize_mixture_fraction_source_terms)
    tp->source.pdfc[0]  = c->rho * area_old * mass_transfer_coeff;
 
  UpdateTemperature(tp, -mp_dot * h_lat, TP_MASS(tp), 1.);
 
  remaining_mass = TP_MASS(tp) - mp_dot * dt;
 
  /* if the particle is very small, or the particle temperature
     exceeds the boiling temperature dump the remaining mass */
  if ((remaining_mass < TP_DPM_MINIMUM_MASS(tp))||( TP_DIAM(tp) <= DPM_LOWEST_DIAMETER)
      || (TP_T(tp)> DPM_BOILING_TEMPERATURE(tp, m)))
  {
    TP_DIAM(tp) = 0.0;
    TP_T(tp) = Tp_old;
    TP_MASS(tp) = 0.;
  }
  else
  {
    TP_MASS(tp) -= mp_dot * dt;
    TP_DIAM(tp) = DPM_DIAM_FROM_VOL(TP_MASS(tp) / TP_RHO(tp));
  }
}

DEFINE_DPM_SOURCE(dpm_source, c, t, S, strength, tp)
{
   /* delta_m is the mass source to the continuous phase
    * (Difference in mass between entry and exit from cell)
    * multiplied by strength (Number of particles/s in stream)
    */
  real delta_m = (TP_MASS0(tp) - TP_MASS(tp)) * strength;

  if (TP_CURRENT_LAW(tp) == DPM_LAW_USER_1)
  {
       /* Sources relevant to the user law 1:
        *
        * add the mixture fraction source and adjust the energy source by
        * adding the latent heat at reference temperature
        *
        */
    S->pdf_s[0] += delta_m;
    S->energy = -delta_m * TP_INJECTION(tp)->latent_heat_ref;
  }
}

After the UDF that you have defined using DEFINE_DPM_SOURCE is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_DPM_SOURCE UDFs for details on how to hook your DEFINE_DPM_SOURCE UDF to Ansys Fluent.

You can use DEFINE_DPM_SPRAY_COLLIDE to side-step the default Ansys Fluent spray collision algorithm. When droplets collide they may bounce (in which case their velocity changes) or they may coalesce (in which case their velocity is changed, as well as their diameter and number in the DPM parcel). A spray collide UDF is called during droplet tracking after every droplet time step and requires that Droplet Collision is enabled in the Discrete Phase Model dialog box.

DEFINE_DPM_SPRAY_COLLIDE (name, tp, p)

Function returns

void

There are three arguments to DEFINE_DPM_SPRAY_COLLIDE: name, tp, and p. You supply name, the name of the UDF. tp and p are variables that are passed by the Ansys Fluent solver to your UDF. When collision is enabled, this linked list is ordered by the cell that the particle is currently in. As particles from this linked list are tracked, they are copied from the particle list into a Tracked_Particle structure.

The following UDF, named mean_spray_collide, is a simple (and non-physical) example that demonstrates the usage of DEFINE_SPRAY_COLLIDE. The droplet diameters are assumed to relax to their initial diameter over a specified time t_relax. The droplet velocity is also assumed to relax to the mean velocity of all droplets in the cell over the same time scale.

/***********************************************************
   DPM Spray Collide Example UDF
 ************************************************************/
 #include "udf.h"
 #include "dpm.h"
 #include "surf.h"
 DEFINE_DPM_SPRAY_COLLIDE(mean_spray_collide,tp,p)
 {
    /* non-physical collision UDF that relaxes the particle */
    /* velocity and diameter in a cell to the mean over the */
    /* specified time scale t_relax */
    const real t_relax = 0.001; /* seconds */

    /* get the cell and Thread that the particle is currently in */
    cell_t c = TP_CELL(tp);
    Thread *t = TP_CELL_THREAD(tp);
    /* Particle index for looping over all particles in the cell */
    Particle *pi;

    /* loop over all particles in the cell to find their mass */
    /* weighted mean velocity and diameter */
    int i;
    real u_mean[3]={0.}, mass_mean=0.;
    real d_orig = TP_DIAM(tp);
    real decay = 1. - exp(-t_relax);
    begin_particle_cell_loop(pi,c,t)
      {
         mass_mean += PP_MASS(pi);
         for(i=0;i<3;i++)
           u_mean[i] += PP_VEL(pi)[i]*PP_MASS(pi);
      }
    end_particle_cell_loop(pi,c,t)  /* relax particle velocity to the mean and diameter to the */
    /* initial diameter over the relaxation time scale t_relax */
    if(mass_mean > 0.)
      {
       for(i=0;i<3;i++)
         u_mean[i] /= mass_mean;
       for(i=0;i<3;i++)
         TP_VEL(tp)[i] += decay*(u_mean[i] - TP_VEL(tp)[i]);
         TP_DIAM(tp) += decay*(TP_INIT_DIAM(tp) - TP_DIAM(tp));
         /* adjust the number in the droplet parcel to conserve mass */
         TP_N(tp) *= CUB(d_orig/TP_DIAM(tp));
      }
 } 

After the UDF that you have defined using DEFINE_DPM_SPRAY_COLLIDE is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Discrete Phase Model dialog box in Ansys Fluent.

See Hooking DEFINE_DPM_SPRAY_COLLIDE UDFs for details on how to hook your DEFINE_DPM_SPRAY_COLLIDE UDF to Ansys Fluent.

You can use DEFINE_DPM_SWITCH to modify the criteria for switching between laws. The function can be used to control the switching between the user-defined particle laws and the default particle laws, or between different user-defined or default particle laws.

DEFINE_DPM_SWITCH (name, tp, ci)

Function returns

void

There are three arguments to DEFINE_DPM_SWITCH: name, tp, and ci. You supply name, the name of the UDF. tp and ci are variables that are passed by the Ansys Fluent solver to your UDF.

The following is an example of a compiled UDF that uses DEFINE_DPM_SWITCH to switch between DPM laws using a criterion. The UDF switches to DPM_LAW_USER_1 which refers to condenshumidlaw since only one user law has been defined. The switching criterion is the local humidity which is computed in the domain using a DEFINE_ON_DEMAND function, which again calls the function myHumidity for every cell. In the case where the humidity is greater than , condensation is computed by applying a simple mass transfer calculation. Otherwise, Ansys Fluent’s standard law for Inert Heating is applied. The UDF requires one UDML and needs a species called h2o to compute the local humidity.

/**********************************************************************
   Concatenated UDFs for the Discrete Phase Model including
   an implementation of a condensation model
   an example for the use of DPM_SWITCH
 ***********************************************************************/
 
 #include "udf.h"
 #include "dpm.h"

 #define UDM_RH 0           /* no. of UDM holding relative humidity */
 #define N_REQ_UDM 1        /* 1 more than UDM_RH */
 #define CONDENS 1.0e-4     /* a condensation rate constant */
 
 int h2o_index=0;           /* index of water vapor species in mixture material */
 real mw_h2o=18.;           /* molecular weight of water */
 
 real H2O_Saturation_Pressure(real T)
 {
   real ratio, aTmTp;
   T = MAX(T, 273);
   T = MIN(T, 647.286);
   aTmTp = .01 * (T - 338.15);
   ratio = (647.286 / T - 1.) *
           (-7.419242 + aTmTp * (.29721 +
                        aTmTp * (-.1155286 +
                        aTmTp * (8.685635e-3 +
                        aTmTp * (1.094098e-3 +
                        aTmTp * (-4.39993e-3 +
                        aTmTp * (2.520658e-3 -
                        aTmTp * 5.218684e-4)))))));
  return (22.089e6 * exp(MIN(ratio, 35.)));
 }
 
 real myHumidity(cell_t c, Thread *t)
 {
   int i;
   Material *m = THREAD_MATERIAL(t), *sp;
   real yi_h2o = 0;      /* water mass fraction */
   real r_mix = 0.0;     /* sum of [mass fraction / mol. weight] over all species */
   real humidity;

   if ((MATERIAL_TYPE(m) == MATERIAL_MIXTURE) && (FLUID_THREAD_P(t)))
     {
       yi_h2o = C_YI(c, t, h2o_index);     /* water vapor mass fraction */
       mixture_species_loop(m, sp, i)
         {
           r_mix += C_YI(c,t,i) / MATERIAL_PROP(sp, PROP_mwi);
         }

       humidity = op_pres * yi_h2o / (mw_h2o * r_mix) /
                  H2O_Saturation_Pressure(C_T(c,t));
       return humidity;
     }
   else
     return 0.;
 }
 
 DEFINE_DPM_LAW(condenshumidlaw, tp, coupled)
 {
   real area;
   real mp_dot;

   /* Get Cell and Thread from Particle Structure */
   cell_t c = TP_CELL(tp);    
   Thread *t = TP_CELL_THREAD(tp); 

   area = M_PI * (TP_DIAM(tp) * TP_DIAM(tp));

   /* Note This law only used if Humidity > 1.0 so mp_dot always positive*/
   mp_dot = CONDENS * sqrt(area) * (myHumidity(c, t) - 1.0);
   if (mp_dot > 0.0)
     {
       TP_MASS(tp) += mp_dot * TP_DT(tp); 
       TP_DIAM(tp) = pow(6.0 * TP_MASS(tp) / (TP_RHO(tp) * M_PI), 1./3.);
     }
   /* Assume condensing particle is in thermal equilibrium with fluid in cell */
   TP_T(tp) = C_T(c,t);
 }
 
 DEFINE_DPM_SOURCE(dpm_source, c, t, S, strength, tp) 
 { 
   real mp_dot; 

   /* mp_dot is the mass source to the continuous phase
    * (Difference in mass between entry and exit from cell)
    * multiplied by strength (Number of particles/s in stream)
    */
   mp_dot = (TP_MASS0(tp) - TP_MASS(tp)) * strength;

   if (TP_CURRENT_LAW(tp) == DPM_LAW_USER_1)
     {
       /* Sources relevant to the user law 1:
        * add the source to the condensing species
        * equation and adjust the energy source by 
        * adding the latent heat at reference temperature
        */
       S->species[h2o_index] += mp_dot;
       S->energy -= mp_dot * TP_INJECTION(tp)->latent_heat_ref;
     }
 } 
  
 DEFINE_DPM_SWITCH(dpm_switch, tp, coupled)
 {
   cell_t c = TP_CELL(tp);
   Thread *t = TP_CELL_THREAD(tp);
   Material *m = TP_MATERIAL(tp); 

   /* If the relative humidity is higher than 1 
    * and the particle temperature below the boiling temperature 
    * switch to condensation law
    */
   if ((C_UDMI(c,t,UDM_RH) > 1.0) && (TP_T(tp) < DPM_BOILING_TEMPERATURE(tp, m)))
     TP_CURRENT_LAW(tp) = DPM_LAW_USER_1;
   else
     TP_CURRENT_LAW(tp) = DPM_LAW_INITIAL_INERT_HEATING;
 }

 DEFINE_ADJUST(adj_relhum, domain)
 {
   cell_t cell;
   Thread *thread;

   if(N_UDM < N_REQ_UDM)
      Message("\nNot enough user defined memory allocated. %d required.\n",
      N_REQ_UDM);
   else
     {
       real humidity, min, max;
       min = 1e10;
       max = 0.0;
       thread_loop_c(thread, domain)
         {
           /* Check if thread is a Fluid thread and has UDMs set up on it */
           if (FLUID_THREAD_P(thread) && NNULLP(THREAD_STORAGE(thread, SV_UDM_I)))
             {
               Material *m = THREAD_MATERIAL(thread), *sp;
               int i;
               /* Set the species index and molecular weight of water */ 
               if (MATERIAL_TYPE(m) == MATERIAL_MIXTURE)
                 mixture_species_loop (m,sp,i)
                   {
                     if (0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"h2o") ||
                        (0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"H2O")))
                       {
                         h2o_index = i;
                         mw_h2o = MATERIAL_PROP(sp,PROP_mwi);
                       }
                   }

               begin_c_loop(cell,thread)
                 {
                   humidity = myHumidity(cell, thread);
                   min = MIN(min, humidity);
                   max = MAX(max, humidity);
                   C_UDMI(cell, thread, UDM_RH) = humidity;
                 }
               end_c_loop(cell, thread)
             }
         }
       Message("\nRelative Humidity set in udm-%d", UDM_RH);
       Message(" range:(%f,%f)\n", min, max);
     }/* end if for enough UDSs and UDMs */
 }
 
 DEFINE_ON_DEMAND(set_relhum)
 {
   adj_relhum(Get_Domain(1));
 } 

After the UDF that you have defined using DEFINE_DPM_SWITCH is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible in the Custom Laws dialog box in Ansys Fluent. See Hooking DEFINE_DPM_SWITCH UDFs for details on how to hook your DEFINE_DPM_SWITCH UDF to Ansys Fluent.

You can use DEFINE_DPM_TIMESTEP to change the time step for DPM particle tracking based on user-specified inputs. The time step can be prescribed for special applications where a certain time step is needed. It can also be limited to values that are required to validate physical models.

DEFINE_DPM_TIMESTEP (name, tp, ts)

Function returns

real

There are three arguments to DEFINE_DPM_TIMESTEP: name, tp, and ts. You supply the name of your user-defined function. tp and ts are variables that are passed by the Ansys Fluent solver to your UDF. Your function will return the real value of the DPM particle timestep to the solver.

The following compiled UDF named limit_to_e_minus_four sets the time step to a maximum value of 1e-4. If the time step computed by Ansys Fluent (and passed as an argument) is smaller than 1e-4, then Ansys Fluent’s time step is returned.

/* Time step control UDF for DPM */
 
 #include "udf.h"
 #include "dpm.h"
 
 DEFINE_DPM_TIMESTEP(limit_to_e_minus_four,tp,dt)
 {
   if (dt > 1.e-4)
    {
      /* TP_NEXT_TIME_STEP(tp) = 1.e-4; */
      return 1.e-4;
    }
 return dt;
 } 

The following compiled UDF named limit_to_fifth_of_prt computes the particle relaxation time based on the formula:

(2–21)

where

(2–22)

The particle time step is limited to a fifth of the particle relaxation time. If the particle time step computed by Ansys Fluent (and passed as an argument) is smaller than this value, then Ansys Fluent’s time step is returned.

When executing the below UDF with high-resolution particle tracking, it is recommended that the continuous phase viscosity be interpolated to the particle position if it varies with position. This can be enabled with the following text command:

define/models/dpm/numerics/high-resolution-tracking/barycentric-interpolation/interpolate-flow-viscosity? yes

/* Particle time step control UDF for DPM */
 
 #include "udf.h"
 #include "dpm.h"
 
 DEFINE_DPM_TIMESTEP(limit_to_fifth_of_prt,tp,dt)
 {
   real drag_factor = 0.;
   real p_relax_time;
   cphase_state_t *c = &(tp->cphase[0]);
   /* compute particle relaxation time */
   if (TP_DIAM(tp) != 0.0)
      drag_factor = DragCoeff(tp) * c->mu / (TP_RHO(tp) * TP_DIAM(tp) * TP_DIAM(tp));
   else
      drag_factor = 1.;
   p_relax_time = 1./drag_factor;
   /* check the condition and return the time step */
   if (dt > p_relax_time/5.)
    {
      return p_relax_time/5.;
    }
   return dt;
 } 

After the UDF that you have defined using DEFINE_DPM_TIMESTEP is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible and selectable for DPM Timestep in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_DPM_TIMESTEP UDFs for details on how to hook your DEFINE_DPM_TIMESTEP UDF to Ansys Fluent.

You can use DEFINE_DPM_VP_EQUILIB to specify the equilibrium vapor pressure of vaporizing components of multicomponent particles.

DEFINE_DPM_VP_EQUILIB (name, tp, T, cvap_surf, Z)

Function returns

void

There are five arguments to DEFINE_DPM_VP_EQUILIB: name, tp, T, cvap_surf, and Z. You supply the name of your user-defined function. tp and T are passed by the Ansys Fluent solver to your UDF. Your UDF will need to compute the equilibrium vapor concentration and the compressibility factor and store their values in cvap_surf and Z, respectively.

The following UDF named raoult_vpe computes the equilibrium vapor concentration of a multicomponent particle using the Raoult law. The vapor pressure in the law is proportional to the molar fraction of the condenses material. DEFINE_VP_EQUILIB is called several times every particle time step in Ansys Fluent and requires a significant amount of CPU time to execute. For this reason, the UDF should be executed as a compiled UDF.

/***********************************************************************
   UDF for defining the vapor particle equilibrium
   for multicomponent particles
 ***********************************************************************/
 #include <udf.h>
 DEFINE_DPM_VP_EQUILIB(raoult_vpe,tp,Tp,cvap_surf,Z)
 {
    int is;
    real molwt[MAX_SPE_EQNS];
    Thread *t0 = TP_CELL_THREAD(tp);  /* cell thread of particle location */
    Material *gas_mix = THREAD_MATERIAL(t0); /* gas mixture material */
    Material *cond_mix = TP_MATERIAL(tp); /* particle mixture material */
    int nc = TP_N_COMPONENTS(tp);    /* number of particle components */
    real molwt_cond = 0.;  /* reciprocal molecular weight of the particle */
    for (is = 0; is < nc; is++)
      {
         int gas_index = TP_COMPONENT_INDEX_I(tp,is); /* index of vaporizing
            component in the gas phase */
         if (gas_index >= 0)
           {
            /* the molecular weight of particle material */
            molwt[gas_index] =
           MATERIAL_PROP(MIXTURE_COMPONENT(gas_mix,gas_index),PROP_mwi);
            molwt_cond += TP_COMPONENT_I(tp,is) / molwt[gas_index];
           }
        }
    /* prevent division by zero */
    molwt_cond = MAX(molwt_cond,DPM_SMALL);
    for (is = 0; is < nc; is++)  {
     /* gas species index of vaporization */
     int gas_index = TP_COMPONENT_INDEX_I(tp,is);
     if(gas_index >= 0)
       {
        /* condensed material */
        Material * cond_c = MIXTURE_COMPONENT(cond_mix, is);
        /* condensed component molefraction */
        real xi_cond = TP_COMPONENT_I(tp,is)/(molwt[gas_index]*molwt_cond);
        /* particle saturation pressure */
        real p_saturation = DPM_vapor_pressure(tp, cond_c, Tp);
        if (p_saturation < 0.0)
          p_saturation = 0.0;
        /* vapor pressure over the surface, this is the actual Raoult law */
        cvap_surf[is] = xi_cond * p_saturation / UNIVERSAL_GAS_CONSTANT / Tp;
       }
    }
  /* compressibility for ideal gas */
  *Z = 1.0;
 } 

After the UDF that you have defined using DEFINE_DPM_VP_EQUILIB is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name of the argument that you supplied as the first DEFINE macro argument will become visible and selectable in the Create/Edit Materials dialog box in Ansys Fluent. Note that before you hook the UDF, you’ll need to create particle injections in the Injections dialog box with the type Multicomponent chosen. See Hooking DEFINE_DPM_VP_EQUILIB UDFs for details on how to hook your DEFINE_DPM_VP_EQUILIB UDF to Ansys Fluent.

You can use DEFINE_IMPINGEMENT to customize the impingement regime selection criteria.

DEFINE_IMPINGEMENT (name, tp, rel_dot_n, f, t, y_s, E_imp)

Function returns

int

There are seven arguments to DEFINE_IMPINGEMENT: name, tp, real_dot_n, f, t, y_s, and E_imp. You supply the name of your user-defined function. tp, real_dot_n, f, and t are passed by the Ansys Fluent solver to your UDF. Your function will:

The impingement regime that you return can either be one of the four predefined impingement regimes in Ansys Fluent, described in Interaction During Impact with a Boundary in the Fluent Theory Guide:

or one of the following user-defined regimes defined in a DEFINE_FILM_REGIME UDF:

The following UDF, named dry_impingement, returns one of the pre-defined regimes, or the user-defined regime FILM_USER_0 for high wall temperature and high impact energy conditions. This function must be run as a compiled UDF.

#include "udf.h"
#define E_crit_0 16.
#define E_crit_1 3329.
#define E_crit_2 7500.
#define T_crit 1.5
DEFINE_IMPINGEMENT(dry_impingement, tp, rel_dot_n, f, t, y_s, E_imp)
{
  real one = 1.0, zero = 0.0;
  real fh = F_WALL_FILM_HEIGHT(f,t);
  real Re, denom;
  real abs_visc  = MAX(DPM_MU(tp),DPM_SMALL);
  real sigma = MAX(DPM_SURFTEN(tp),DPM_SMALL);
  real p_rho = MAX(TP_RHO(tp),DPM_SMALL);
  real d_0   = TP_DIAM(tp);
  real T_w = rf_energy ? F_T(f,t) : 300.0;
  real T_b = dpm_par.Tmax;
  int regime = FILM_STICK;
  
  *y_s = zero;
  Re = p_rho * d_0 * rel_dot_n / abs_visc;
  denom = sigma * (MIN(fh/d_0,1.0) + pow(Re,-0.5));
  *E_imp = SQR(rel_dot_n) * d_0 * p_rho / MAX(DPM_SMALL, denom);

  /* rf_energy is defined if Energy model is enabled. It is possible to use inert particle
  with isothermal conditions (Energy off). In that case the UDF returns FILM_STICK. */
  if (rf_energy && TP_LIQUID(tp))
    {
      Material *mb;
      if (TP_WET_COMBUSTION(tp))
        mb = TP_WET_COMB_MATERIAL(tp);
      else
        mb = TP_MATERIAL(tp);
      T_b = DPM_Boiling_Temperature(tp,mb);
    }
  else
    return FILM_STICK;
  
  if (T_w <= T_crit*T_b)
    {
      if (*E_imp <= E_crit_0)
        regime = FILM_STICK;
      else if (*E_imp > E_crit_0 && *E_imp <=  E_crit_1)
        regime = FILM_SPREAD;
      else if (*E_imp > E_crit_1)
        regime = FILM_SPLASH;
    }
  else
    {
      if (*E_imp < E_crit_1)
        regime = FILM_REBOUND;
      else
        regime = FILM_USER_0;
    }
  
  if ((regime == FILM_SPREAD) || (regime == FILM_STICK))
    *y_s = zero;
  else if ((regime == FILM_REBOUND)|| (regime == FILM_USER_0))
    *y_s = one;
  else if (regime == FILM_SPLASH)
    {
      if (*E_imp < E_crit_2)
        *y_s = MAX(0., 1.8e-4*(*E_imp - E_crit_1));
      else
        *y_s = 0.7;
    }
  return regime;
}

After the UDF that you have defined using DEFINE_IMPINGEMENT is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name that you supplied as the first DEFINE macro argument will become visible and selectable in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_IMPINGEMENT UDFs for details on how to hook your DEFINE_IMPINGEMENT UDF to Ansys Fluent.

You can use DEFINE_FILM_REGIME to set the particle variables for user-defined regimes of particle impingement. Fluent supports up to 10 user-defined regimes:

DEFINE_FILM_REGIME (name, regime, tp, pp, f, t, f_normal, update)

Function returns

void

There are eight arguments to DEFINE_FILM_REGIME: name, regime, tp, pp, f, t, f_normal, and update. You supply the name of your user-defined function. regime, tp, pp, f, t, f_normal, and update are passed by the Ansys Fluent solver to your UDF. Your function will need to compute the particle parameters and variables corresponding to the user-defined regime and modify the Tracked_Particle *tp and Particle *pp structures accordingly.

The following UDF, named dry_breakup, describes a user-defined impingement regime where the impinging droplet particle bounces off the wall surface and breaks up into smaller droplets without any film formation. The ratio of the initial droplet diameter to the diameter after reflection is defined by the parameter D_RATIO. The UDF uses the Ansys Fluent function Reflect_Particle to compute the particle velocities after impingement. The reflected droplet’s new diameter and mass, as well as the number in parcel and flowrate are updated in the Tracked_Particle *tp structure. Note that the mass and diameter of the current and previous time steps, as well as the initial values must be updated. The initial mass and diameter must also be adjusted in Particle structure *pp at the end of the tracking step.

#include "udf.h"
#define D_RATIO 5
#define DIAM_FROM_VOL(V) ((real)pow((6.0 * (V) / M_PI), 1./3.))

DEFINE_FILM_REGIME(dry_breakup, regime, tp, pp, f, t, f_normal, update)
{
  real d1 = TP_DIAM(tp);
  real d2 = MAX(TP_DIAM(tp)/D_RATIO,dpm_par.lowest_diam);
  real drat_cubed;

  if (regime == FILM_USER_0)
    { 
      UNMARK_TP(tp, P_FL_ON_WALL);
      tp->tracking_scheme = tp->free_stream_tracking_scheme;
      TP_FILM_FACE(tp) = NULL_FACE;
      TP_FILM_THREAD(tp) = NULL;
      tp->gvtp.n_rebound ++ ;
      Reflect_Particle(tp,f_normal,ND_ND,f,t, NULL, NULL_FACE);
      drat_cubed = CUB(d2/d1);      
      TP_DIAM(tp) = d2;
      TP_N(tp) /= drat_cubed;
      TP_MASS(tp) *= drat_cubed;
      TP_FLOW_RATE(tp) /= drat_cubed;
      /* after diameter change all relevant variables in previous and initial states must be consistent */
      /* initialize the previous state mass and diameter to the current */
      TP_MASS0(tp) = TP_MASS(tp);
      TP_DIAM0(tp) = TP_DIAM(tp);
      /* initialize the initial state mass and diameter to the current */
      TP_INIT_MASS(tp) *= drat_cubed;
      TP_INIT_DIAM(tp) = DIAM_FROM_VOL(  TP_INIT_MASS(tp) /  TP_INIT_RHO(tp) );
      /* at the end of the iteration step update the initial mass
         and diameter in the stored particle list */
      if (update)
        {
          PP_INIT_MASS(pp) *= drat_cubed;
          PP_INIT_DIAM(pp) = DIAM_FROM_VOL(  PP_INIT_MASS(pp) /  PP_INIT_RHO(pp) );      
        }
    }
  else
    {
      Error("User defined impingement regime not implemented: %d.\n",regime);  
    }
}

After the UDF that you have defined using DEFINE_FILM_REGIME is interpreted (Interpreting UDFs) or compiled (Compiling UDFs), the name that you supplied as the first DEFINE macro argument will become visible and selectable in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_FILM_REGIME UDFs for details on how to hook your DEFINE_FILM_REGIME UDF to Ansys Fluent.

You can use DEFINE_SPLASHING_DISTRIBUTION to customize the diameter, number in parcel, and velocity distribution for the splashed particles.

DEFINE_SPLASHING_DISTRIBUTION (name, tp, rel_dot_n, f_normal, n_splash, splashing_distribution_t *s)

Function returns

void

There are six arguments to DEFINE_SPLASHING_DISTRIBUTION: name, tp, real_dot_n, f_normal, n_splash, and s. You supply the name of your user-defined function. tp, real_dot_n, f_normal, and n_splash are passed by the Ansys Fluent solver to your UDF. Note that n_splash is the value entered for Number of Splashed Drops in the DPM tab of the boundary condition dialog box when Discreate Phase BC Function is set to wall-film. The output of your function is the array of pointers to the structure *s containing the distribution of the n_splash particle parcels.

The following UDF applies the O’Rourke & Amsden splashing model [10]. This function must be run as a compiled UDF.

#include "udf.h"
#define NUKIYAMA_TANASAWA_TABLE_SIZE 15

static double nukiyama_tanasawa_table[NUKIYAMA_TANASAWA_TABLE_SIZE][2] =
  {
    {0.000000000000, 0.000000000000},
    {0.081108621634, 0.500000000000},
    {0.193908811137, 0.700000000000},
    {0.266112345561, 0.800000000000},
    {0.345136495805, 0.900000000000},
    {0.427593287415, 1.000000000000},
    {0.510077497807, 1.100000000000},
    {0.589500744183, 1.200000000000},
    {0.663337543595, 1.300000000000},
    {0.729766711799, 1.400000000000},
    {0.787709698251, 1.500000000000},
    {0.877181682854, 1.700000000000},
    {0.934793421924, 1.900000000000},
    {0.953988295926, 2.000000000000},
    {1.000000000000, 4.000000000000},
  };

double get_interpolated_value(double table[][2], int table_size, double inval)
{
  /* table is table[tableSize][1] with table[i][0] being the lookup values
     and table[i][1] the return values */

  int lowerindex, upperindex, midindex, jump;
  double lowerval, upperval, midval, ratio;

  lowerindex = 0;
  upperindex = table_size-1;

  lowerval = table[lowerindex][0];
  upperval = table[upperindex][0];

  jump = (upperindex-lowerindex)/2;

  /* Find values for Interpolation */
        
  while (jump)
    {
      midindex = lowerindex + jump;
      midval = table[midindex][0];

      if (inval > midval)
	{
	  lowerval = midval;
	  lowerindex = midindex;
	}
      else
	{
	  upperval = midval;
	  upperindex = midindex;
	}
      jump /= 2;
    }

  ratio = (inval-lowerval)/(upperval - lowerval);
  lowerval = table[lowerindex][1];
  upperval = table[upperindex][1];
        
  return lowerval + ratio*(upperval - lowerval);
}

double nukiyama_tanasawa_random()
{
  return get_interpolated_value(nukiyama_tanasawa_table, NUKIYAMA_TANASAWA_TABLE_SIZE, uniform_random());
}

double nabor_reitz_random(double beta)
{
  double exp_beta;
  double rand;
  double angle_range;

  angle_range = M_PI;
  rand = 1.0 - 2.0*uniform_random(); /* -1.0 < rand <= 1.0 */

  if (beta == 0.0)
    return angle_range * rand;
  
  if (rand < 0.0)
    {
      angle_range = -M_PI;
      rand = -rand;    /* 0.0 <= rand <= 1.0 */
    }
  exp_beta = exp(beta);
  
  return angle_range*(beta - log(exp_beta + rand*(1.0 - exp_beta)))/beta;
}

/*  O'Rourke & Amsden splashing model BC definition */

#define E_2_CRIT 3329.29 /* O'Rourke & Amsden critical energy squared (57.7^2) */

DEFINE_SPLASHING_DISTRIBUTION(splash, tp, rel_dot_n, norm, n_splashes, s)
{
  int i;
  real p_diam, p_v_mag; /* Particle properties */
  real w0, v0, v_t;     /* Particle condition  */
  real E_2;
  real sin_alpha;
  
  real e_t[ND_ND], e_p[ND_ND];
  /* Splashed particle properties */
  real weber;
  real d_max, w_dash, v_dash;                 
  real beta, psi, component;
  real total_volume = 0.;
      
  p_diam = TP_DIAM(tp);

  /* norm into wall and TP_VEL should be so w0 >= 0.0 */
  w0 = NV_DOT(TP_VEL(tp),norm);
  
  /* Particle impingement energy */
  E_2 = TP_IMPINGEMENT_PARAMETER(tp);
  
  /* Calculate new splash stream droplet diameters */ 
  weber = TP_WEBER_IMP(tp);
  
  d_max = MAX(MAX((E_2_CRIT/E_2),(6.4/weber)),0.06)*p_diam;
  
  for (i=0;i<n_splashes;i++)
    {
      SPLASH_DIAM(s,i) = d_max * nukiyama_tanasawa_random();
    }
  
  /* Calculate droplet velocities */
  NV_CROSS(e_p,TP_VEL(tp),norm);
  p_v_mag = NV_MAG(TP_VEL(tp));
  
  NV_S(e_p,/=,p_v_mag);
  sin_alpha = NV_MAG(e_p);

  /* If particle hits exactly normal to surface tilt slightly */ 
  sin_alpha = MAX(1.E-6, sin_alpha);
  NV_S(e_p,/=,sin_alpha);
  
  beta = M_PI*sqrt(sin_alpha/MAX(DPM_SMALL,(1.0-sin_alpha)));
  
  NV_CROSS(e_t,norm,e_p);
  
  v0 = NV_DOT(TP_VEL(tp),e_t); /* Tangential to wall */
  v_t = 0.8*v0;
  
  for (i=0;i<n_splashes;i++)
    {
      w_dash = -0.2 * w0 * nukiyama_tanasawa_random();
      NV_VS(SPLASH_VEL(s,i),=,norm,*,w_dash);
      
      v_dash = 0.1 * w0 * fabs(gauss_random()); /* O'Rourke and Amsden: delta = 0.1 w0 */
      v_dash += 0.12*w0;
      psi = nabor_reitz_random(beta);
      component = v_dash * cos(psi);
      NV_VS(SPLASH_VEL(s,i),+=,e_t,*,component);
      component = v_dash * sin(psi);
      NV_VS(SPLASH_VEL(s,i),+=,e_p,*,component);
      
      NV_VS(SPLASH_VEL(s,i),+=,e_t,*,v_t);
    }
  
  /*  The flow rate for each diameter will be proportional to its mass (=volume as density equal) */
  for (i=0;i<n_splashes;i++)
    {
      total_volume += CUB( SPLASH_DIAM(s,i));
    }
  for (i=0;i<n_splashes;i++)
    {
      SPLASH_N(s,i) = TP_SPLASHED_FRACTION(tp)*CUB(SPLASH_DIAM(s,i))/total_volume*TP_N(tp);
    }
} 

After the UDF that you have defined using DEFINE_SPLASHING_DISTRIBUTION is compiled (Compiling UDFs) the name that you supplied as the first DEFINE macro argument will become visible and selectable in the Discrete Phase Model dialog box in Ansys Fluent. See Hooking DEFINE_SPLASHING_DISTRIBUTION UDFs for details on how to hook your DEFINE_SPLASHING_DISTRIBUTION UDF to Ansys Fluent.