前期推送了多个关于燃烧和化学反应的案例，从应用范围来讲，通用有限速率模型适应更多的应用场景。比如在氢能行业，关于燃料电池电化学的模拟、关于重整制氢的模拟等等，就需要考虑详细的化学反应。因此，我们觉得有必要对详细化学反应的模拟做一些更深入的探讨，接下来我们以一个具体的例子，上一个案例研究了FLUENT自身模型求解详细化学反应的关键内容，今天用UDF来模拟化学反应并与软件自身计算结果进行对比。我们考虑以下的化学反应，这是甲烷蒸汽重整制氢的反应之一，反应焓等于摩尔加权的反应物和生成物的生成焓之差，本反应的反应焓为正直，说明是一个吸热反应。FLUENT材料物性的标准状态焓就是材料的生成焓，单位为J/kmol，因此本反应的反应焓转化为FLUENT的反应焓就是206.3×106J/kmol。这个反应的字面理解就是每消耗1kmol的CH4或者每消耗1kmol的H2O或者每生成1kmol的CO或者每生成3kmol的H2，需要吸收206.3×106J的能量。当然，也可以理解为每消耗1/3kmol的CH4或者每消耗1/3kmol的H2O或者每生成1/3kmol的CO或者每生成1kmol的H2，需要吸收1/3×206.3×106J的能量，也即是前文所述的反应焓是摩尔加权计算结果，注意数值上一定要对应。通用有限速率模型中，FLUENT通过在能量方程中增加由于化学反应产生的源项来计算反应热，式5.10中各项的单位为Sh,rxn:W/m3，hj:J/kmol、Mj:kg/kmol、Rj:kg/m3s。该式也等效于ΔH(J/kmol)×R(kmol/m3s)。我们假定反应速率如下，反应速率式的各个变量、各个指数的单位与FLUENT的单位吻合。另外，根据反应式，我们可以得出各个组分的摩尔速率关系为：rCH4=rH2O=rCO=rI，rH2=3rI。也就是各组分摩尔速率之比就等于化学计量数之比，这是化学反应基本常识。在化学反应实验中，通常是得出某个物质的反应速率，也就得到了整个反应以及其他物质的反应速率。FLUENT通过反应摩尔速率式乘以各反应组分的化学计量数，等到该组分的反应摩尔速率，注意反应物的速率为负（表示消耗）、生成物的速率为正（表示生成）。根据上面的反应速率式，在FLUENT中输入对应的参数。如果我们不通过UDF定义反应速率，那么FLUENT采用上面的参数进行计算。如果我们用DEFINE_VR_RATE这个宏来定义反应速率，那么上面的参数将不再决定反应速率，而由宏来决定。DEFINE_VR_RATE宏的组成如下，注意该宏指定的反应速率特指体积反应（volumetric reaction）。该宏有8个变量，分别如下，读者朋友自行阅读，懒得翻译了，能用到这个宏的一定看得懂如下英文了。除了第一个为UDF的名字以外，其他都是由FLUENT求解器传递给UDF的变量，当然不是每个变量都用得上，具体得看速率和这些变量有没有关系了。另外，这是一个无返回值的宏。当采用层流有限速率/涡耗散模型时，取两者的较小值作为反应速率，这两个速率都要计算，如果是层流有限速率模型，可以只计算层流速率。我们用如下的UDF定义该反应速率，在之前的案例中，我们讨论了本征动力学方程在软件中的输入问题，如果所有组分的速率指数都起到作用，那么该UDF其实就是完全定义上述的反应速率式（摩尔速率），然后由于我们采用的是层流有限速率模型，这里简单粗暴地将湍流速率和层流速率强行划等号，当然也可以把湍流速率这行代码删除。该UDF代码可读性很强，容易理解，我们就不过多解释了。需要说明的是，反应速率可以根据实际情况编写，也就是说反应速率不受输入面板上的阿累尼乌斯的参数以及各组分速率指数的限制，而由本UDF决定，这就可以应用到那些非标准的阿累尼乌斯形式的反应速率问题了。DEFINE_VR_RATE(vol_reac_rate,c,t,r,wk,yk,rate,rr_t){real ci, prod1,cj,prod2;int i,j;prod1 = 1.;prod2 = 1.;for(i = 0; i < r->n_reactants; i++){ci = C_R(c,t) * yk[r->reactant[i]] / wk[r->reactant[i]];prod1 *= pow(ci, r->exp_reactant[i]);}for(j = 0; j < r->n_products; j++){cj = C_R(c,t) * yk[r->product[j]] / wk[r->product[j]];prod2 *= pow(cj, r->exp_product[j]);}*rate = r->A * exp(- r->E / (UNIVERSAL_GAS_CONSTANT * C_T(c,t))) *pow(C_T(c,t), r->b) * prod1*prod2;
但是如果生成物的速率指数不起作用的话，UDF代码就需要如下了，也就是删除了生成物的影响。DEFINE_VR_RATE(vol_reac_rate1,c,t,r,wk,yk,rate,rr_t){real ci, prod;int i;prod = 1.;for(i = 0; i < r->n_reactants; i++){ci = C_R(c,t) * yk[r->reactant[i]] / wk[r->reactant[i]];prod *= pow(ci, r->exp_reactant[i]);}*rate = r->A * exp(- r->E / (UNIVERSAL_GAS_CONSTANT * C_T(c,t))) *pow(C_T(c,t), r->b) * prod;*rr_t = *rate;}
通过UDF外挂面板在体积反应速率加载该函数，注意用编译的方式加载。另外，特别说明一下，对于用户自定义的反应速率，特别是温度指数β不为零的情况，如果用stiff chemistry solver则经常容易出现dasac failure错误，有时候这种错误会导致发散，有时候迭代几步就消失，如果消失了则可以无视，如果一直存在，则是有问题的。可考虑将化学求解器改成None-Direct source，比如本案例，如果不做这个修改，则dasac failure一直存在。我们设定反应入口的组分为H2O和CH4，并且摩尔比重为0.6和0.4。其他的相关设置我们就不再细说了，读者朋友可以参考以前的案例。这里需要指出的是，为了能够顺利自启动反应，我们给反应入口设定一个高的温度773.15K。首先，检查一下质量守恒情况，进出口满足质量守恒，这是最基础的，没有这个平衡，整个模拟计算可执行一票否决。然后，我们看一下进出口的各组分的质量分数，和上一个案例FLUENT自身计算的几乎一样。最后，检查一下能量守恒情况。反应热-2.45W，和上一个案例FLUENT自身计算的几乎一样。进出口温度和壁面温度以及温度云图如下，和上一个案例FLUENT自身计算的几乎一样。 
