FLUENT-UDF代码-颗粒(煤或生物质)床层多相流燃烧过程数值模拟-【第2部分-完结】

1、Fluent-UDF_Coal-Combustion-Multiphase-Flow-Processes-Codes for various kinds of heterogeneous reactions【第二部分-Part2】DEFINE_HET_RXN_RATE(SteamGasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t) Thread *pt = THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ int index_phase = Get_Phase_Index(hr); Thread *ts = p

2、tindex_phase; /* solid phase */ *rr = 0; double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/* CO + H2 -> H2O + 1/25 Soot The reverse steam gasification reaction, which is CO + H2 -> 1/25 Soot + H2O is written as the reaction shown above. So a negative rate means CO and H2 is consum

3、ed and H2O and Soot are generated. Note that no C(s) is generated, the stoich coeff for C(s) in the above reation is zero. 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. Th

4、is is done by the first reaction, defined in the heterogeneous reaction panel in FLUENT GUI. */ if(IS_ASH = 0) SetSpeciesIndex(); if(MGAS_Gasif) double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); double p_h2o = RoRT * yiIP_H2OIS_H2O/mwIP_H2OIS_H2O/ 101325.; double p_co = RoRT * yiIP_COIS_

5、CO/mwIP_COIS_CO / 101325.; double p_h2 = RoRT * yiIP_H2IS_H2/mwIP_H2IS_H2 / 101325.; y_carbon = yiIP_SOOTIS_SOOT;mol_weight = mwIP_SOOTIS_SOOT; if (rp_ke) rr_turb = Turbulent_rr(c, t, hr, yi); if(C_VOF(c, ts) >= eps_s_small) *rr = rr_steam_gasif(c, t, ts, tp, p_h2o, p_co, p_h2, y_carbon, mol_weig

6、ht, &direction); /* mol/(cm3 .s) */ if( direction > 0.0) /* positive value implies C(s) + H2O -> CO + H2 */ *rr = 0.0; else /* negative value implies CO + H2 -> H2O + 1/25 Soot */ *rr = abs(*rr); *rr = MIN(*rr, rr_turb); double rr_steam_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp,

7、 double p_h2o, double p_co, double p_h2, double y_carbon, double mol_weight, double* direction) double rate, prod, T_g = MIN(MAX(TMIN,C_T(c,tp),TMAX); double p_h2o_star = p_h2 * p_co / ( exp(17.29 - 16326/T_g) ); if(MGAS_Gasif) *direction = p_h2o - p_h2o_star; if(PCCL_Gasif) *direction = pow(p_h2o,

8、N_steam_gasification)/(1.+K_steam_gasification*p_h2); prod = y_carbon*(C_R(c,ts)*1e-03)/mol_weight*C_VOF(c,ts); /*1e-3 is to convert density from kg/m3 to g/cm3 */ if(MGAS_Gasif && *direction < 0.0) /* this implies reverse H2O gasification */ prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C

9、_VOF(c,tp); /*1e-3 is to convert density from kg/m3 to g/cm3 */ rate = A_steam_gasification*exp(-E_steam_gasification/Rgas/T_g)* Annealing_steam_gasification * prod * *direction; /* mol/cm3.s */ rate *= 1000.; /* kmol/(m3 .s) */ return rate;DEFINE_HET_RXN_RATE(Co2Gasif,c,t,hr,mw,yi,rr,rr_t) Thread *

10、pt = THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ int index_phase = Get_Phase_Index(hr); Thread *ts = ptindex_phase; /* solid phase */ *rr = 0; double direction = 0.0, mol_weight, y_carbon; /* C(s) + CO2 -> 2CO Set the phase and species indices. Ash species index is initialized to zer

11、o, 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(); double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); doubl

12、e p_co = RoRT * yiIP_COIS_CO/mwIP_COIS_CO / 101325.; double p_co2 = RoRT * yiIP_CO2IS_CO2/mwIP_CO2IS_CO2 / 101325.; SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight); if(C_VOF(c, ts) >= eps_s_small) *rr = rr_co2_gasif(c, t, ts, tp, p_co, p_co2, y_carbon, mol_weight, &direction); /*

13、 mol/(cm3 .s) */ if( direction < 0.0) /* negative implies reverse steam gasification */ *rr = 0.0; DEFINE_HET_RXN_RATE(Soot_CO2_Gasif,c,t,r,mw,yi,rr,rr_t) Thread *pt = THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ *rr = 0; double rr_turb = 1e+20; double T_g = MIN(MAX(TMIN,C_T(c,tp),TMA

14、X); /* 1/25 Soot + CO2 -> 2CO 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 GU

15、I. */ if(IS_ASH = 0) SetSpeciesIndex(); double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); double p_co = RoRT * yiIP_COIS_CO/mwIP_COIS_CO/ 101325.; double p_co2 = RoRT * yiIP_CO2IS_CO2/mwIP_CO2IS_CO2 / 101325.; if (rp_ke) rr_turb = Turbulent_rr(c, t, r, yi); double prod = yiIP_SOOTIS_SOOT

16、*(C_R(c,tp)*1e-03)/mwIP_SOOTIS_SOOT*C_VOF(c,tp); /*1e-3 is to convert density from kg/m3 to g/cm3 */ *rr = A_soot_co2_gasification*exp(-E_soot_co2_gasification/Rgas/T_g)* Annealing_soot_co2_gasification * prod * pow(p_co2, N_soot_co2_gasification)/(1.+K_soot_co2_gasification*p_co); /* mol/cm3.s */ *

17、rr *= 1000.; /* kmol/(m3 .s) */ *rr = MIN(*rr, rr_turb); DEFINE_HET_RXN_RATE(Co2Gasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t) Thread *pt = THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ int index_phase = Get_Phase_Index(hr); Thread *ts = ptindex_phase; /* solid phase */ *rr = 0; double direction =

18、0.0, mol_weight, y_carbon, rr_turb = 1e+20; /* 2CO -> CO2 + 1/25 Soot The reverse CO2 gasification reaction, which is 2CO -> 1/25 Soot + CO2 is written as the reaction shown above. So a negative rate means CO is consumed and CO2 and Soot are generated. Note that no C(s) is generated, the stoic

19、h coeff for C(s) in the above reation is zero. 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 pane

20、l in FLUENT GUI. */ if(IS_ASH = 0) SetSpeciesIndex(); if(MGAS_Gasif) double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); double p_co = RoRT * yiIP_COIS_CO/mwIP_COIS_CO / 101325.; double p_co2 = RoRT * yiIP_CO2IS_CO2/mwIP_CO2IS_CO2 / 101325.; y_carbon = yiIP_SOOTIS_SOOT; mol_weight = mwIP_S

21、OOTIS_SOOT; if (rp_ke) rr_turb = Turbulent_rr(c, t, hr, yi); if(C_VOF(c, ts) >= eps_s_small) *rr = rr_co2_gasif(c, t, ts, tp, p_co, p_co2, y_carbon, mol_weight, &direction); /* mol/(cm3 .s) */ if( direction > 0.0) /* positive value implies C(s) + CO2 -> 2CO */ *rr = 0.0; else /* negativ

22、e value implies 2CO -> CO2 + 1/25 Soot */ *rr = abs(*rr); *rr = MIN(*rr, rr_turb); double rr_co2_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_co, double p_co2, double y_carbon, double mol_weight, double* direction) double T_g = MIN(MAX(TMIN,C_T(c,tp), TMAX), prod; double p_co2_star

23、 = p_co * p_co/(exp(20.92 - 20282/T_g); if(MGAS_Gasif) *direction = p_co2-p_co2_star; if(PCCL_Gasif) *direction = pow(p_co2, N_co2_gasification)/(1. + K_co2_gasification * p_co); prod = y_carbon*C_R(c,ts)*1.e-3/mol_weight* C_VOF(c,ts); /*1e-3 is to convert density from kg/m3 to g/cm3 */ if(MGAS_Gasi

24、f && *direction < 0.0) /* this implies reverse CO2 gasification */ prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m3 to g/cm3 */ double rate = A_co2_gasification*exp(-E_co2_gasification/Rgas/T_g)*Annealing_co2_gasification * prod * (*directi

25、on); /* mol/cm3.s */ rate *= 1000.; /* kmol/(m3 .s) */ return rate; DEFINE_HET_RXN_RATE(H2Gasif,c,t,hr,mw,yi,rr,rr_t) Thread *pt = THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ int index_phase = Get_Phase_Index(hr); Thread *ts = ptindex_phase; /* solid phase */ *rr = 0; double direction =

26、 0.0, mol_weight, y_carbon; /* 1/2 C(s) + H2 -> 1/2 CH4 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 r

27、eaction panel in FLUENT GUI. */ if(IS_ASH = 0) SetSpeciesIndex(); double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); double p_h2 = RoRT * yiIP_H2IS_H2/mwIP_H2IS_H2 / 101325.; double p_ch4 = RoRT * yiIP_CH4IS_CH4/mwIP_CH4IS_CH4 / 101325.; SolidFuel_Reactant(c, t, hr, &y_carbon, &mo

28、l_weight); if(C_VOF(c, ts) >= eps_s_small) *rr = rr_h2_gasif(c, t, ts, tp, p_h2, p_ch4, y_carbon, mol_weight, &direction); /* mol/(cm3 .s) */ if( direction < 0.0) /* negative implies reverse steam gasification */ *rr = 0.0; DEFINE_HET_RXN_RATE(Soot_H2_Gasif,c,t,r,mw,yi,rr,rr_t) Thread *pt

29、= THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ *rr = 0; double rr_turb = 1e+20; double T_g = MIN(MAX(TMIN,C_T(c,tp),TMAX); /* 1/25 Soot + 2H2 -> CH4 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

30、 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(); double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp); double p_h2 = RoRT * yiIP_H2IS_H2/mwIP_H2IS_H2 / 101325.; if (rp_ke

31、) rr_turb = Turbulent_rr(c, t, r, yi); double prod = yiIP_SOOTIS_SOOT*(C_R(c,tp)*1e-03)/mwIP_SOOTIS_SOOT*C_VOF(c,tp); /*1e-3 is to convert density from kg/m3 to g/cm3 */ *rr = A_soot_h2_gasification*exp(-E_soot_h2_gasification/Rgas/T_g)* Annealing_soot_h2_gasification * prod * pow(p_h2, N_soot_h2_ga

32、sification); /* mol/cm3.s */ *rr *= 1000.; /* kmol/(m3 .s) */ *rr = MIN(*rr, rr_turb); DEFINE_HET_RXN_RATE(H2Gasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t) Thread *pt = THREAD_SUB_THREADS(t); Thread *tp = pt0; /* gas phase */ int index_phase = Get_Phase_Index(hr); Thread *ts = ptindex_phase; /* solid phase */

33、 *rr = 0; double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20; /* 1/2 CH4 -> H2 + (0.5)*1/25 Soot The reverse H2 gasification reaction, which is 1/2 CH4 -> 1/25 Soot + H2 is written as the reaction shown above. So a negative rate means CH4 is consumed and H2 and Soot are generated. N

34、ote that no C(s) is generated, the stoich coeff for C(s) in the above reation is zero. 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, de

35、fined in the heterogeneous reaction panel in FLUENT GUI. */ 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 * yiIP_H2IS_H2/mwIP_H2IS_H2 / 101325.; double p_ch4 = RoRT * yiIP_CH4IS_CH4/mwIP_CH4IS_CH4 / 101325.; y_carbon

36、 = yiIP_SOOTIS_SOOT; mol_weight = mwIP_SOOTIS_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/(cm3 .s) */ if( direction > 0.0) /* positive value implies 1/2 C(s) +

37、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 =

38、 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/m3 to g/cm3 */ if(MGAS_Gasif) *direction = p_h2-p_h2_star; if(*direction < 0.0) /* t

39、his implies reverse H2 gasification */ prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m3 to g/cm3 */ rate = exp( -7.087 - 8078/T_g )* prod * *direction ; /* mol/cm3.s */ if(PCCL_Gasif) *direction = p_h2; rate = A_h2_gasification*exp(-E_h2_gasification/

40、Rgas/T_g)*Annealing_h2_gasification * prod * *direction; /* mol/cm3.s */ rate *= 1000.; /* kmol/(m3 .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 = pt0; /* gas phase */ int index_phase = Get_Phase_Index(hr); Thread *ts = p

41、tindex_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

42、 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 = yiindex_phaseIS_ASH; *rr = rr_combustion(c, t, ts, tp, yiIP_O2IS_O2, y_ash, y_carbon); /* mol/(cm3 .s)

43、*/ *rr *= 1000.; /* kmol/(m3 .s) */ double rr_combustion(cell_t c, Thread *t, Thread *ts, Thread *tp, double yi_O2, double y_ash, 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.cm3/mol.K */ doub

44、le 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/mwIP_O2IS_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.; el

45、se rd = 0.; double diff = MAX(4.26 * pow(T/1800.),1.75)/Pt), 1.e-10); /* cm2/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,t

46、p)-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.

47、+ 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.; else k_f = diff * N_sherwood / (D_p * gas_constant/mwIP_O2IS_O2 * T ); /* g/(atm.cm2.s) */ k_r = A_c_combustion * exp( -E_c_combus

48、tion/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/mwIP_O2IS_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.);

49、 /* mol/(cm3 .s) */ return rate;#if !RP_NODE | !PARALLELvoid volatile_mass_fractions() read_c3m_data(); /* pan2 : Oct 2012 . added CX_Messages for debugging */ CX_Message("PCCL_Devol = %dn",PCCL_Devol); CX_Message("MGAS_Devol = %dn",MGAS_Devol); CX_Message("CPD_Devol = %dn",CPD_Devol); CX_Message("FGDVC_Devol = %dn",FGDVC_