发布时间 : 星期二 文章FLUENT UDF代码 - 颗粒（煤或生物质）床层多相流燃烧过程数值模拟-第2部分-完结更新完毕开始阅读
if(IS_ASH == 0)
SetSpeciesIndex(); if(MGAS_Gasif) {
double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;
double p_ch4 = RoRT * yi[IP_CH4][IS_CH4]/mw[IP_CH4][IS_CH4] / 101325.;
y_carbon = yi[IP_SOOT][IS_SOOT]; mol_weight = mw[IP_SOOT][IS_SOOT];
if(C_VOF(c, ts) >= eps_s_small) {
if (rp_ke)
rr_turb = Turbulent_rr(c, t, hr, yi);
*rr = rr_h2_gasif(c, t, ts, tp, p_h2, p_ch4, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */
if( direction > 0.0) /* positive value implies 1/2 C(s) + H2 ---> 1/2 CH4 */ *rr = 0.0; else /* negative value implies 1/2 CH4 ---> H2 + (0.5)*1/25 Soot */ { *rr = abs(*rr); *rr = MIN(*rr, rr_turb); } } } }
double rr_h2_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_h2, double p_ch4, double y_carbon, double mol_weight, double* direction) {
double rate = 0.0, prod;
double T_g = MIN((MAX(TMIN,C_T(c,tp))), TMAX);
double p_h2_star = pow ((p_ch4/(exp(-13.43 + 10999/T_g))), 0.5);
prod = y_carbon*C_R(c,ts)*1.e-3/mol_weight * C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */
if(MGAS_Gasif) {
*direction = p_h2-p_h2_star;
if(*direction < 0.0) /* this implies reverse H2 gasification */
9
prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */
rate = exp( -7.087 - 8078/T_g )* prod * *direction ; /* mol/cm^3.s */ }
if(PCCL_Gasif) {
*direction = p_h2;
rate = A_h2_gasification*exp(-E_h2_gasification/Rgas/T_g)*Annealing_h2_gasification * prod * *direction; /* mol/cm^3.s */ }
rate *= 1000.; /* kmol/(m^3 .s) */ return rate; }
DEFINE_HET_RXN_RATE(coal_combustion,c,t,hr,mw,yi,rr,rr_t) {
Thread **pt = THREAD_SUB_THREADS(t); Thread *tp = pt[0]; /* gas phase */
int index_phase = Get_Phase_Index(hr);
Thread *ts = pt[index_phase]; /* solid phase */ double mol_weight, y_carbon, y_ash; *rr = 0.0;
/* Set the phase and species indices. Ash species index is initialized to zero, with all other indices.
Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the first
reaction, defined in the heterogeneous reaction panel in FLUENT GUI. */
if(IS_ASH == 0)
SetSpeciesIndex();
if( C_YI(c,tp,IS_O2) >= spe_small) {
SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight); y_ash = yi[index_phase][IS_ASH];
*rr = rr_combustion(c, t, ts, tp, yi[IP_O2][IS_O2], y_ash, y_carbon); /* mol/(cm^3 .s) */ *rr *= 1000.; /* kmol/(m^3 .s) */ } }
double rr_combustion(cell_t c, Thread *t, Thread *ts, Thread *tp, double yi_O2, double y_ash,
10
double y_carbon) {
double rd, k_f, k_r, factor, k_a, rate = 0.0, vrel; double Pt = MAX(0.1, (op_pres+C_P(c,t))/101325); double gas_constant = 82.06; /* atm.cm^3/mol.K */
double T = C_T(c,tp), T_s = C_T(c,ts), D_p = C_PHASE_DIAMETER(c,ts)*100.;
double p_o2 = C_R(c,tp)*UNIVERSAL_GAS_CONSTANT* T *yi_O2/mw[IP_O2][IS_O2] / 101325.; /* atm */
if(fc_ar > 0.) {
if (y_ash > 0.) {
rd = pow( (y_carbon * ash_ar/100.)/(y_ash * fc_ar/100.), (1./3.) ); rd = MIN(1., rd); } else
rd = 1.; } else
rd = 0.;
double diff = MAX((4.26 * pow((T/1800.),1.75)/Pt), 1.e-10); /* cm^2/s */ double Sc1o3 = pow(C_MU_L(c,tp)/(C_R(c,tp) * diff * 1.e-4), 1./3.); #if RP_2D
vrel = pow(( (C_U(c,tp)-C_U(c,ts))*(C_U(c,tp)-C_U(c,ts)) + (C_V(c,tp)-C_V(c,ts))*(C_V(c,tp)-C_V(c,ts))), 0.5); #endif #if RP_3D
vrel = pow(( (C_U(c,tp)-C_U(c,ts))*(C_U(c,tp)-C_U(c,ts)) + (C_V(c,tp)-C_V(c,ts))*(C_V(c,tp)-C_V(c,ts)) +
(C_W(c,tp)-C_W(c,ts))*(C_W(c,tp)-C_W(c,ts)) ), 0.5); #endif
double Re = C_VOF(c,tp) * D_p/100. * vrel * C_R(c,tp)/(C_MU_L(c,tp)+SMALL_S); double N_sherwood = (7. - 10. * C_VOF(c,tp) + 5. * C_VOF(c,tp) * C_VOF(c,tp) )* (1. + 0.7 * pow(Re, 0.2) * Sc1o3) +
(1.33 - 2.4 * C_VOF(c,tp) + 1.2 * C_VOF(c,tp) * C_VOF(c,tp)) * pow(Re, 0.7) * Sc1o3;
if ( rd <= 0. || C_VOF(c, ts) <= 0. ) {
rate = 0.; }
11
else {
k_f = diff * N_sherwood / (D_p * gas_constant/mw[IP_O2][IS_O2] * T ); /* g/(atm.cm^2.s) */
k_r = A_c_combustion * exp( -E_c_combustion/Rgas/T_s ) * rd * rd; if ( rd >= 1.) {
rate = 1. / (1./k_f + 1./k_r); } else {
k_a = 2. * rd * diff * f_ep_a / (D_p * (1.-rd) * gas_constant/mw[IP_O2][IS_O2] * T_s );
rate = 1. / (1./k_f + 1./k_r + 1./k_a); }
factor = y_carbon / (y_carbon + 1.e-6);
rate *= p_o2 * 6. * C_VOF(c,ts) * factor / (D_p * 32.); /* mol/(cm^3 .s) */ } return rate; }
#if !RP_NODE || !PARALLEL
void volatile_mass_fractions() {
read_c3m_data();
/* pan2 : Oct 2012 ... added CX_Messages for debugging */ CX_Message(\ CX_Message(\ CX_Message(\
CX_Message(\ CX_Message(\
CX_Message(\ CX_Message(\
CX_Message(\ CX_Message(\ CX_Message(\ CX_Message(\
CX_Message(\ CX_Message(\
CX_Message(\
CX_Message(\
12