求各位大神指点迷津一个添加源项的UDF，小弟感激不尽……

hu521pan

#include "udf.h"
#include "sg_mphase.h"
#define A 16.0114
#define B 3456.80
#define C -8.67
#define L 0.00000125.0
#define POP 101325.0
#define K 0.1
#define C10H22_NUM 0
#define LIQ_PHASE_ID 0
#define GAS_PHASE_ID 1

DEFINE_SOURCE(liq_zhengkuiwan,c,liq_thread,dS,eqn) //添加c10h22液滴蒸发燃烧的源项，计算的是二维
{
int i;
cell_t c;
Thread*gas_thread,*mixture_thread;
Domain*mixture_domain,*liq_domain,*gas_domain;
Material *gas_material,*species;
real c10h22_mole_fract,P_c10h22,total_mole,Mw[MAX_SPE_EQNS];
mixture_domain=Get_Domain(1);
liq_domain=DOMAIN_SUB_DOMAIN(mixture_domain,LIQ_PHASE_ID);
gas_domain=DOMAIN_SUB_DOMAIN(mixture_domain,GAS_PHASE_ID);  //感觉没什么用，第一次编写udf就都写上了
mixture_thread=THREAD_SUPER_THREAD(liq_thread);
gas_thread=THREAD_SUB_THREAD(mixture_thread,GAS_PHASE_ID);
gas_material=THREAD_MATERIAL(gas_thread);
real m_dot_liq,c_s,c_g,m_kuosan;
real x[ND_ND];
real ps,v;  //ps为液滴表面饱和压力 v为单个网格的体积

thread_loop_c(mixture_thread,mixture_domain) //在整个二维控制区查询混合物的单元线
{
begin_c_loop(c,mixture_thread)        //在单元线上查询单元
{if(0.01<=C_VOF(c,liq_thread)&&C_VOF(c,liq_thread)<1.0)   //查询气液界面处的网格，我选了0.01到1之间，因为当液滴的体积分数小于0.01时会产生沸腾蒸发
{
ps=133*exp(A-B/(C_T(c,liq_th)+C);           //Antoine蒸气压方程求得液滴表面饱和压力ps
c_s=ps/(C_RGAS(c,gas_thread)*C_T(c,liq_thread));   //求得液滴表面浓度  C_RGAS(c,gas_thread)气体常数这样可以吗？？也可以直接输入8.3134

total_mole=0;
mixture_species_loop(gas_material,species,i)
{
Mw=MATERIAL_PROP(species,PROP_mwi);
total_mole+=C_YI(c,gas_thread,i)/Mw;           //求得气液界面网格内气体的总摩尔数
}
c10h22_mole_fract=(C_YI(c,gas_thread,C10H22_NUM)/Mw[C10H22_NUM])/total_mole;   //求得气液界面内c10h22的摩尔分数
P_c10h22=POP*c10h22_mole_fract;                         //求得气液界面内c10h22的分压
c_g=P_c10h22/(C_RGAS(c,gas_thread)*C_T(c,liq_thread));         //求得气液界面网格内气体中c10h22的浓度
m_kuosan=c_s-c_g;            //气液界面网格内液滴表面c10h22浓度和气体中c10h22浓度的差值，从而产生扩散
v=L*L*L;
C_CENTROID(x,c,mixture_thread);
y=x[1];
m_dot_liq=142*K*y*m_kuosan*C_DIFF_EFF(c,gas_thread,i)/v*C_VOF(c,liq_th)*pow(1-C_VOF(c,liq_thread),2.0)*(1-(1-C_VOF(c,liq_thread))*c10h22_mole_fract);
dS[eqn]=          //最终源项      dS还没有求，比较麻烦
}
}
end_c_loop(c,mixture_thread)
return m_dot_liq;
}