/***********************************************************************
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))) /* NULLP returns TRUE if storage is not allocated foruser-defined storage variable */
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, p, 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] = P_POS(p)[1]; /* Radial location vector. */
radi_pos[1] = P_POS(p)[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] = P_VEL(p)[0]; /* Axial particle velocity component. */
imp_vel[1] = NVD_DOT(radi_pos, P_VEL(p)[1], P_VEL(p)[2], 0.);
}
else
#endif
NV_V(imp_vel, =, P_VEL(p));
/* 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_PARTICLE(p, 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) = (P_DIAM(p)+ 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();
}
评论0