#include "udf.h"
#define WALL_ID 9 // the wall id to be computed collision angle
#define PARTICLE_PHASE_ID 1 // the phase id of particle phase
#define ANGLE_UDM_ID 0 // the UDM No. to store collision angle, the first UDM is 0
#define NUM_UDM 1
int count = 0;
real start_avg = 0.005;
DEFINE_EXECUTE_AT_END(execute_at_end)
{
if(CURRENT_TIME > start_avg)
{
face_t f;
cell_t c0;
real cos_angle,angle;
real NV_VEC(f_normal),NV_VEC(particle_velocity);
Domain*mixture_domain,*particle_domain;
Thread*wall_thread,*particle_t0;
mixture_domain=Get_Domain(1);
particle_domain=DOMAIN_SUB_DOMAIN(mixture_domain,PARTICLE_PHASE_ID);
wall_thread=Lookup_Thread(particle_domain,WALL_ID);
count = count+1;
if (FLUID_THREAD_P(wall_thread))
{
begin_c_loop(c0,wall_thread)
{
thread_loop_c(wall_thread,particle_domain)
{
begin_f_loop(f,wall_thread)
{
c0=F_C0(f,wall_thread);
particle_t0=THREAD_T0(wall_thread);
#if RP_3D
NV_D(particle_velocity,=,C_U(c0,particle_t0),C_V(c0,particle_t0),C_W(c0,particle_t0));
#else
NV_D(particle_velocity,=,C_U(c0,particle_t0),C_V(c0,particle_t0),0);
#endif
F_AREA(f_normal,f,wall_thread);
cos_angle=NV_DOT(particle_velocity,f_normal)/MAX(NV_MAG(particle_velocity)*NV_MAG(f_normal),1E-10);
angle=90.0-acos(cos_angle)/M_PI*180.0;
Message0("face id=%d,angle between particle velocity and wall face=%f degree\n",f,angle);
C_UDMI(c0,particle_t0,ANGLE_UDM_ID)=angle;
}
end_f_loop(f,wall_thread);
}
C_UDMI(c0,wall_thread,0) = C_UDMI(c0,particle_t0,ANGLE_UDM_ID)+C_UDMI(c0,wall_thread,0);
C_UDMI(c0,wall_thread,1) = C_UDMI(c0,wall_thread,0)/count ;
}
end_c_loop(c0,wall_thread)
}
}
}
评论4