有偿求助DPM模型中heat mass transfer模块怎么应用。

yi_chao19791001

有偿求助DPM模型中heat mass transfer模块怎么应用。我写了段程序，并请大家修改，并告诉我怎么加载。
#include "udf.h"
#define B 1319
#define M 0.69
#define N 0.40
#define NP1 3.44
#define NP2 3.28
#define Sc 0.00296
#define R 8.314
DEFINE_DPM_HEAT_MASS(multivap,p,Cp,hgas,hvap,cvap_surf,dydt,dzdt)
{
real hm,pw, pva, tw, Re, xa, ta, ga, xw, Sm1, Sm2, Sm3, kh1, kh2, kh3, lef, vp, x[ND_ND];
cell_t c0 = RP_CELL(&(p->cCell));
Thread *t0 = RP_THREAD(&(p->cCell));
Material *cond_mix=p->injection->material;
cphase_state_t *c=&(p->cphase);
real rw = 2260.0;
real cpv = 2.014;
real Tp = P_T(p);
real mp = 0.011095187;
real Pr = 0.00061157;
real Dp = DPM_DIAM_FROM_VOL( mp / P_RHO(p) );
real Ap = DPM_AREA(Dp);
Material *cond_c = MIXTURE_COMPONENT( cond_mix, 1 );
real D=MATERIAL_PROP_POLYNOMIAL(cond_c,PROP_binary_diffusivity,c->temp);
xa = C_YI(c0,t0,1) / (1 - C_YI(c0,t0,1) );
ga = C_R(c0,t0)*fabs(C_V(c0,t0));
tw = P_T(p);
ta= C_T(c0,t0);
xw = 3.703 + 0.286*tw + 0.009164*pow(tw,2) + 0.0001446*pow(tw,3) + 0.000001741*pow(tw,4) + 0.00000005195*pow(tw,5);
Sm1 = B * (xw-xa) * pow(mp,M) * pow(ga,N);
C_CENTROID(x,c0,t0);
if (x[1]<=382.3143)
{
  vp = sqrt(19.614*(382.3143-x[1]));
}
if (x[1]<=386.8243&&x[1]>384.3243)
{
  vp = 39.3*(386.8243-x[1])/(sqrt(78.5*(386.8243-x[1])+100)-10)-5;
}
Re = 337.84 * sqrt( pow(C_U(c0,t0),2) + pow(C_V(c0,t0)+vp,2) );
hm = D/Dp * (2+0.6*pow(Re,0.5)*pow(Sc,1/3));
pw = pow(2.71828, 9.3876-3826.36/(tw-45.47));
pva = pow(2.71828, 9.3876-3826.36/(C_T(c0,t0)-45.47));
Sm2 = NP1*Ap*hm*(pw/(R*tw)-pva/(R*C_T(c0,t0)));
Sm3 = NP2*Ap*hm*(pw/(R*tw)-pva/(R*C_T(c0,t0)));
if (x[1]<=382.3143)
{
  dzdt->species[1]=Sm2;
  dydt[1]=0-Sm2;
}
if (x[1]<=384.324&&x[1]>382.3143)
{
  dzdt->species[1]=Sm1;
  dydt[1]=0-Sm1;
}
if (x[1]<=386.8243&&x[1]>384.3243)
{
  dzdt->species[1]=Sm3;
  dydt[1]=0-Sm3;
}
lef = pow(0.865,0.667)*((0.662+xw)/(xa+0.662)-1)/log((xw+0.662)/(xa+0.662));
kh1 = lef * cpv * B * pow(mp,M) * pow(ga,N);
kh2 = NP1 * Ap * 0.0242 * (2.0 + 0.6*pow(Re,0.5)*pow(Pr,1/3))/Dp;
kh3 = NP2 * Ap * 0.0242 * (2.0 + 0.6*pow(Re,0.5)*pow(Pr,1/3))/Dp;
if (x[1]<=382.3143)
{
  dzdt->energy = Sm2*cpv*(tw-ta) + kh2*(tw-ta);
  dydt[0]=(-Sm2*cpv*(tw-ta)-kh2*(tw-ta))/4.182;
}
if (x[1]<=384.324&&x[1]>382.3143)
{
  dzdt->energy = Sm1*cpv*(tw-ta) + kh1*(tw-ta);
  dydt[0]=(-Sm1*cpv*(tw-ta)-kh1*(tw-ta))/4.182;
}
if (x[1]<=386.8243&&x[1]>384.3243)
{  
  dzdt->energy = Sm3*cpv*(tw-ta) + kh3*(tw-ta);
  dydt[0]=(-Sm3*cpv*(tw-ta)-kh3*(tw-ta))/4.182;
}
}

yi_chao19791001

事情很急，如能解决此问题，必有重谢