24小时热门版块排行榜

返回列表

jimoxingzhe

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 17.8
帖子: 4
在线: 6.7小时
虫号: 2109486
注册: 2012-11-06

[求助] 求组大侠给我看看鄙人编写的UDF，感激不尽……小弟没有金币啊呜呜

#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)
{
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],y=x[1];
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;
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;
}

[ Last edited by jimoxingzhe on 2013-4-9 at 16:14 ]

回复此楼