FLUENT沉降池模拟

#include "udf.h"#include "turb.h"
#define SC_T 0.7
#define RHOP 1450. #define RHOW 1000.
#define GX 9.81
#define US0 0.005#define RH 0.7#define RP 5#define CMIN 0.01#define USMAX 0.002
static real settling_velocity(real C){  real US;
if(CUS=0.;elseUS=US0*(exp(-RH*(C-CMIN))-exp(-RP*(C-CMIN)));
if(US>0.002)US=0.002;
return US;}
DEFINE_UDS_FLUX(settling_flux,f,tf,i){  cell_t c0=F_C0(f,tf),c1=F_C1(f,tf);Thread *t0=THREAD_T0(tf),*t1=THREAD_T1(tf);  real face_concentration;  real gravity_vector[]={1,0,0},area[ND_ND];
if ((t1==NULL)||(THREAD_TYPE(tf)==THREAD_F_WALL))      face_concentration=F_UDSI(f,tf,0);else    face_concentration=.5*(C_UDSI(c0,t0,0)+C_UDSI(c1,t1,0)); F_AREA(area,f,tf);
if (THREAD_TYPE(tf)==THREAD_F_WALL)return 0.;elsereturn F_FLUX(f,tf)+C_R(c0,t0)*settling_velocity(face_concentration)*NV_DOT(gravity_vector,area);  }
DEFINE_DIFFUSIVITY(turbulent_diff,c,tc,i){ return 1e-10+C_MU_T(c,tc)/SC_T;}
DEFINE_SOURCE(X_mom_src,c,tc,ds,eqn){ C_UDMI(c,tc,0)=settling_velocity(C_UDSI(c,tc,0));
C_UDMI(c,tc,1)=C_MU_T(c,tc)/SC_T;  
C_UDMI(c,tc,2) = GX*C_UDSI(c,tc,0)*(RHOP-C_R(c,tc))/RHOP;
return C_UDMI(c,tc,2) ;}
DEFINE_SOURCE(turb_k_source,c,tc,ds,eqn){  real beta=(RHOP-RHOW)/(RHOP*RHOW); C_UDMI(c,tc,3)=-GX*beta*C_MU_T(c,tc)/SC_T*C_UDSI_G(c,tc,0)[0];  ds[eqn]=-2*GX*beta*0.09*C_K(c,tc)/(C_D(c,tc)*SC_T)*C_UDSI_G(c,tc,0)[0];
return C_UDMI(c,tc,3);}
DEFINE_SOURCE(turb_e_source,c,tc,ds,eqn){  real C3_eps=0.;  C3_eps=tanh(fabs(C_U(c,tc)/C_V(c,tc))+fabs(SMALL_VP));  C_UDMI(c,tc,4)=C3_eps*1.44*C_D(c,tc)/C_K(c,tc)*C_UDMI(c,tc,3);
return C_UDMI(c,tc,4);}
DEFINE_PROPERTY(viscosity,c,tc){if (C_UDSI(c,tc,0)>0.7)return 0.00327*pow(10.,0.132*C_UDSI(c,tc,0));elseif (C_UDSI(c,tc,0)>0.01)return 0.00404527222+0.00304527222*(C_UDSI(c,tc,0)-0.7)/0.69;elsereturn 0.001;}