当前，一些学者对厌氧反应器流场和水力学性能进行了研究，例如，王卫京等(2010；2011)应用湍流模型模拟了侧伸搅拌式厌氧反应器流场，季军远等(2012)采取多釜串联及扩散模型模拟了分段组合式厌氧反应器行为，邹尧等(2012)通过测定停留时间分布研究了厌氧反应器中不同反应区的流动特征.这些研究取得了一定成果，但由于模型相对简单，对于构造复杂的内循环厌氧反应器不能完全适用.同时，将动力学方程和流动模型结合以模拟厌氧反应器中反应过程的研究则较少有文献(Diamantis et al., 2007；张良佺等，2011；张倩倩等，2013)涉及.基于此，本文在建立厌氧反应动力学方程后，针对自主设计的内循环厌氧反应器采取组合模型进行模拟研究.厌氧反应过程中由于有气体生成，气体从产生到形成气泡，以及气泡在反应器中的向上运动必然会对整个反应器的流动型式产生重大的影响，使得反应器的流型偏离于理想状况.为此，本文通过测定实际反应器中流体的停留时间分布，根据停留时间分布的特征提出流动的组合模型，并据此研究实际的反应过程，以期为内循环厌氧反应器的进一步工业化放大提供有益的开发思路和方法.

接种污泥为某啤酒厂UASB反应器中活性污泥；营养液参照产甲烷活性测定营养液配制方法进行配制；牲畜废水为牲畜粪便经过淘洗、过滤去除大颗粒泥沙及部分悬浮物后所得的液体，粪便取自某养猪场；示踪剂为硝酸钠(分析纯)，由杭州化学试剂有限公司提供.

主要仪器：FA1104N型电子天平、DFG-9053A型电热恒温鼓风干燥箱、TC-15型恒温电热套、SC-15型数控超级恒温槽、BT01-YZ1515型蠕动泵、GC-1690B型气相色谱仪

为建立厌氧反应器的流动模型，实验研究了反应器停留时间分布，并考察了不同气体流量的影响，获得了停留时间分布函数.内循环厌氧反应器由主、次两反应区组成，由于反应主要在主反应区4进行，因此，实验主要针对反应区4进行停留时间分布测定.实验流程如图 1所示.

图1(Fig.1)

图1 工艺流程示意图(1.进水罐，2.蠕动泵，3.气泵，4.主反应区，5.内反应区，6.出水口，7.气体出口，8.反应液出口，9.环隙反应区) Fig.1 Schematic diagram of process flow

以水为工作介质，将水加入进水罐至一定高度，开启蠕动泵，保持水流量为15 L · h^-1，流量计安装在泵出口处；同时开启气泵，气泵流量由安装在泵出口处的气体流量计检测，分别为0、90、180 mL · h^-1.水从下而上流过反应器，反应器分为上下两部分，下部为主反应区，内径为110 mm，高550 mm，几何体积为5 L；上部为次反应区，又分为两部分，一是内反应区5，另一为环隙反应区9，环隙反应区几何体积V_R4为2.44 L，内反应区几何体积V_R5为2.86 L.在环隙反应区的上部安装有气体出口管，由于测定的是主反应区停留时间分布，故液体在主反应区上部出水口处排出，出水口处同时安装一电导仪，通过检测出口液体的电导率而推算出示踪物浓度的高低.

当气液两相流动达到稳定后，在进水罐上方快速加入一定量的示踪剂，同时计时，间隔一定时间，于反应器出口处记录下水的电导率随时间的变化.

验证实验工艺流程示意图如图 2所示.新鲜料液加入储槽中，料液由蠕动泵输送，进入换热器中，换热器由10根有机玻璃管组成，管外径为10 mm，长65 mm，料液走管程，与循环水换热后进入反应器，在次反应区顶部进行气液分离，液体进入废液槽，气体经过气体流量计后由缓冲瓶进入液封槽，液封槽容积为1500 mL，最后进入皂沫流量计计量放空.

图2(Fig.2)

图2 验证实验工艺流程示意图(1.原料液储槽，2.蠕动泵，3.流量计，4.换热器，5.反应器，6.反应废液槽，7.气体流量计，8.液封槽，9.缓冲瓶，10.皂沫流量计) Fig.2 Schematic diagram of technological process in confirmatory experiment

由于新鲜料液配料较为困难，当系统稳定，即气体流量变化不大时，每隔4 h分别测量气液两相组成、流量，一般测量2~3次结束实验.

对反应器主反应区流动状况作出如下假设：示踪物A为瞬间脉冲进料；反应器中气体流量恒定不变，且均匀上升；进料清水流量亦恒定不变.现假设厌氧反应器主反应区流动型态可由一平推流区(PFR)与两个并联的全混区(CSTR)串联组合而成，模型示意图如图 3所示.

图3(Fig.3)

图3 中反应器流动模型组合示意图文名 Fig.3 Schematic diagram of the reactor combination of flow model

设反应器主反应区总体积为V_R，等于模型中各反应区体积V_Ri之和，设PFR的体积分率为f₁，CSTR2、CSTR3的体积分率分别为f₂、f₃ ，总流量为v₀，流经CSTR2的流量分率为I，某时刻流出CSTR2、CSTR3全混区的流体中示踪剂浓度分别为C_A2、C_A3，混合后出口总浓度为C_A.当t=0时，在PFR反应器入口处，脉冲法加入质量为M的饱和示踪剂A.对于PFR：

对CSTR2进行物料衡算，求取CSTR2出口浓度C_A2：①当t＜ _p时，由于示踪剂A仍在PFR中，对CSTR2，C_A2=0；②当t≥$\bar t$_p时，对CSTR2作物料衡算，根据输入量=输出量+累计量，则有：

求解该方程，得

同理，对CSTR3进行物料衡算，可求取CSTR3出口浓度C_A3：①当t＜ $\bar t$_p时，由于CSTR3中无示踪剂，故C_A3=0；②当t≥ $\bar t$_p时，对CSTR3作物料衡算，可得：

最后，对CSTR2、CSTR3出口处作物料衡算，由于该处进行的是混合过程，则：①当t＜ $\bar t$_p时，C_A=0；②当t≥$\bar t$_p时，C_A=IC_A2+(1-I)C_A3.

停留时间分布密度函数E(t)的求取公式如下：

当t ＜$\bar t$_p时，E(t)=0；当t≥$\bar t$_p时，E(t)= .

最后：

对模型求解转化为求解一个三参数I、f₁、f₂的E(t)方程，其中，$\bar t$_p不是独立的.

实验是在清水流量为15 L · h^-1，气体平均流量分别为0、90.0、182.1 mL · h^-1条件下进行测定的，示踪剂脉冲加入，部分实验数据如表 1所示.

表1(Table 1)

表1 不同气体流量时分布电导率随时间变化关系 Table 1 Different gas flow distribution of electrical conductivity with time changing relationship

表1 不同气体流量时分布电导率随时间变化关系 Table 1 Different gas flow distribution of electrical conductivity with time changing relationship 时间/min 电导率/ (mS · m-1) 气体流量/ (mL · h-1) 时间 /min 电导率/ (mS · m-1) 气体流量/ (mL · h-1) 时间/ min 电导率/ (mS · m-1) 气体流量/ (mL · h-1) 0 0 0 0 0 90.0 0 0 182.3 2 0 0 1.7 0 89.6 1.1 0 182.0 4 0 0 3.4 0 89.2 2.2 0 182.4 6 0 0 5.1 0 90.0 3.3 0 182.7 8 0 0 6.8 0.42 90.2 4.4 3.49 181.7 10 0.34 0 7.2 5.59 90.1 5.5 5.38 181.9 10.2 5.01 0 7.6 5.89 89.8 6.6 6.15 181.6 10.4 6.44 0 8.0 5.74 90.0 7.7 5.82 182.4 10.6 6.73 0 8.4 5.61 89.9 8.8 5.47 182.6 10.8 6.77 0 8.8 5.50 90.2 9.9 5.16 182.1 11.0 6.67 0 9.2 5.27 90.0 11.0 4.91 182.3 12.5 6.09 0 10.0 5.01 89.8 12.1 4.63 182.7 14.0 5.45 0 10.8 4.71 89.9 13.2 4.38 181.8 15.5 4.84 0 11.6 4.52 90.0 14.3 4.19 181.7 17.0 4.27 0 12.4 4.32 90.5 15.4 3.97 181.6 18.5 3.78 0 13.2 4.05 89.8 16.5 3.78 182.5 20.0 3.33 0 14.0 3.78 89.6 17.6 3.54 182.3 21.5 2.91 0 14.8 3.67 90.3 18.7 3.32 182.2 23.0 2.57 0 15.6 3.45 90.4 19.8 3.14 181.9 24.5 2.25 0 16.4 3.37 89.8 20.9 2.98 181.8 26.0 1.98 0 17.2 3.04 90.7 22.0 2.79 182.2 27.5 1.73 0 18.8 2.71 90.0 26.0 2.21 181.9 29.0 1.52 0 20.4 2.40 89.8 30.0 1.71 182.0 32.0 1.14 0 22.0 2.18 90.1 34.0 1.32 182.3 35.0 0.88 0 23.6 1.91 90.7 38.0 1.03 182.2 38.0 0.65 0 25.2 1.68 90.0 42.0 0.78 181.9 41.0 0.49 0 26.8 1.52 90.3 46.0 0.59 182.4 44.0 0.34 0 28.4 1.33 89.9 50.0 0.45 181.8 47.0 0.23 0 30.0 1.15 90.3 54.0 0.33 182.0 50.0 0.15 0 33.5 0.92 90.4 58.0 0.24 182.3 53.0 0.08 0 37.0 0.72 89.9 62.0 0.17 182.6 56.0 0.03 0 40.5 0.52 90.4 66.0 0.12 181.9 59.0 0.01 0 44.0 0.40 90.0 70.0 0.05 182.2 62.0 0 0 47.5 0.31 89.8 74.0 0 181.9 51.0 0.24 89.5 54.5 0.16 89.3 58.0 0.12 90.4 61.5 0.08 90.4 65.0 0.05 90.3 68.5 0.03 89.7 72.0 0 90.3

根据，各时间点的E(t)值如表 2所示.应用非线性最小二乘法求解，不同气体流量时的模型参数及方程如表 3所示.

表2(Table 2)

表2 停留时间分布密度函数E(t)随时间变化关系 Table 2 Residence time distribution density function E(t)changing with time

表2 停留时间分布密度函数E(t)随时间变化关系 Table 2 Residence time distribution density function E(t)changing with time 时间 /min E(t)/实验值 E(t)/计算值 时间 /min E(t)/实验值 E(t)/计算值 时间 /min E(t)/实验 E(t)/计算值 0 0 0 0 0 0 0 0 0 2 0 0 1.7 0 0 1.1 0 0 4 0 0 3.4 0 0 2.2 0 0 6 0 0 5.1 0 0 3.3 0 0 8 0 0 6.8 4.92×10-3 0 4.4 28.97×10-3 0 10 3.92×10-3 0 7.2 65.45×10-3 0 5.5 44.66×10-3 0 10.2 57.78×10-3 0 7.6 68.96×10-3 66.02×10-3 6.6 51.06×10-3 51.32×10-3 10.4 74.27×10-3 0 8.0 67.20×10-3 64.26×10-3 7.7 48.32×10-3 48.49×10-3 10.6 77.61×10-3 0 8.4 65.68×10-3 62.54×10-3 8.8 45.41×10-3 45.81×10-3 10.8 78.07×10-3 78.33×10-3 8.8 64.39×10-3 60.88×10-3 9.9 42.84×10-3 43.28×10-3 11.0 76.92×10-3 77.09×10-3 9.2 61.70×10-3 59.26×10-3 11.0 40.76×10-3 40.89×10-3 12.5 70.23×10-3 68.38×10-3 10.0 58.66×10-3 56.14×10-3 12.1 38.44×10-3 38.64×10-3 14.0 62.85×10-3 60.68×10-3 10.8 55.14×10-3 53.20×10-3 13.2 36.36×10-3 36.51×10-3 15.5 55.82×10-3 53.87×10-3 11.6 52.92×10-3 50.41×10-3 14.3 34.78×10-3 34.50×10-3 17.0 49.24×10-3 47.83×10-3 12.4 50.58×10-3 47.78×10-3 15.4 32.96×10-3 32.60×10-3 18.5 43.59×10-3 42.49×10-3 13.2 47.42×10-3 45.28×10-3 16.5 31.38×10-3 30.80×10-3 20.0 38.40×10-3 37.75×10-3 14.0 44.26×10-3 42.91×10-3 17.6 29.39×10-3 29.11×10-3 21.5 33.56×10-3 33.56×10-3 14.8 42.97×10-3 40.67×10-3 18.7 27.56×10-3 27.51×10-3 23.0 29.64×10-3 29.83×10-3 15.6 40.39×10-3 38.56×10-3 19.8 26.07×10-3 25.99×10-3 24.5 25.95×10-3 26.53×10-3 16.4 39.45×10-3 36.55×10-3 20.9 24.74×10-3 24.56×10-3 26.0 22.83×10-3 23.60×10-3 17.2 35.59×10-3 34.65×10-3 22.0 23.16×10-3 23.21×10-3 27.5 19.95×10-3 21.00×10-3 18.8 31.73×10-3 31.15×10-3 26.0 18.35×10-3 18.90×10-3 29.0 17.53×10-3 18.70×10-3 20.4 28.10×10-3 28.01×10-3 30.0 14.20×10-3 15.39×10-3 32.0 13.15×10-3 14.83×10-3 22.0 25.52×10-3 25.19×10-3 34.0 10.96×10-3 12.54×10-3 35.0 10.15×10-3 11.77×10-3 23.6 22.36×10-3 22.67×10-3 38.0 8.55×10-3 10.22×10-3 38.0 7.50×10-3 9.35×10-3 25.2 19.67×10-3 20.40×10-3 42.0 6.48×10-3 8.33×10-3 41.0 5.65×10-3 7.43×10-3 26.8 17.80×10-3 18.36×10-3 46.0 4.90×10-3 6.79×10-3 44.0 3.92×10-3 5.92×10-3 28.4 15.57×10-3 16.53×10-3 50.0 3.74×10-3 5.54×10-3 47.0 2.65×10-3 4.71×10-3 30.0 13.46×10-3 14.89×10-3 54.0 2.74×10-3 4.52×10-3 50.0 1.73×10-3 3.76×10-3 33.5 10.77×10-3 11.86×10-3 58.0 1.99×10-3 3.68×10-3 53.0 0.923×10-3 2.99×10-3 37.0 8.43×10-3 9.45×10-3 62.0 1.41×10-3 3.01×10-3 56.0 0.346×10-3 2.39×10-3 40.5 6.09×10-3 7.55×10-3 66.0 0.99×10-3 2.45×10-3 59.0 0.115×10-3 1.91×10-3 44.0 4.68×10-3 6.03×10-3 70.0 0.42×10-3 2.01×10-3 62.0 0 1.53×10-3 47.5 3.63×10-3 4.83×10-3 74.0 0 1.63×10-3 51.0 2.81×10-3 3.87×10-3 54.5 1.87×10-3 3.11×10-3 58.0 1.40×10-3 2.49×10-3 61.5 0.94×10-3 2.01×10-3 65.0 0.59×10-3 1.61×10-3 68.5 0.35×10-3 1.31×10-3 72.0 0 1.05×10-3

表3(Table 3)

表3 不同气体流量时流动模型参数计算值 Table 3 Different gas flow when the flow model parameters calculated

表3 不同气体流量时流动模型参数计算值 Table 3 Different gas flow when the flow model parameters calculated 气体流量vg/(mL · h-1) f1 f2 f3 I 0 0.54 0.11 0.35 0.22 90.0 0.38 0.40 0.22 0.43 182.1 0.33 0.66 0.01 0.63 注：f1=fp=1-f2-f3.

当v_g=0、90.0、182.1 mL · h^-1时，分别获得E(t)方程.E(t)方程获得后，将不同时刻t时的C_A值，代入方程，可获得E(t)的计算值，其值列于表 2中.由表 3中数据可知，随着气体流量的增大，反应器中存在较大的返混，PFR体积所占比率f₁不断减小，说明气体流量对反应器的流型有着较大的影响.

1)反应器反应过程模型建立

根据以上停留时间分布的特点，厌氧反应器由主、次反应区组成，因此，可用图 4中的组合模型模拟反应器反应过程行为.图 4所示组合模型中，根据停留时间分布规律，1、2、3拟合的是厌氧反应器主反应区；次反应区由中间回流区和环隙区构成，中间回流区可作为平推流处理，环隙区由于有大量产物气体流过，可作为全混区处理，所以反应器次反应区可由4、5并联进行拟合.

图4(Fig.4)

图4 反应器反应过程模型组合示意图 Fig.4 Schematic diagram of reactor process model

模型中各参数的确定：针对本反应器，实际反应过程中生成的气量在90~180 mL · h^-1之间，因此，模型中主反应区各参数取其平均值，具体如下：f₁=0.36，f₂=0.53，f₃=0.11，I₁=0.53.

次反应区流量分率为：

各反应区有效体积V_Ri的确定：总污泥加入量为5 L，上层即次反应区污泥所占分率按0.35计算，则：V_R1+ V_R2+ V_R3=3.25L，V_R1=1.17L，V_R2=1.7225L，V_R3=0.3575L，V_R4+ V_R5=1.75L，V_R4=V₄/(V₄+ V₅)=0.8057L，V_R5=0.9443L.

2)模型验证

为验证上述模型是否正确，进行了3组不同初始浓度操作状态的实验，当厌氧反应器操作稳定后，即在一定时间能气体流量基本不变，此时，操作参数分别为：废液体积流量v_0L=15 L · h^-1，温度T=35 ℃，废液初始COD C_S0=3.3840、3.5326、3.8729 g · L^-1，由前期实验所得动力学方程为：

式中，C_s为反应器中废液浓度(g · L^-1)，-r_s为废液COD消耗速率(g · min^-1 · L^-1)，r_pCH4为甲烷生成速率(g · min^-1 · L^-1)，r_pCO2为二氧化碳生成速率(g · min^-1 · L^-1)，模拟计算了反应器单程转化过程，主、次反应区出口即B、C两点处计算值如表 4所示.以基质初始浓度C_S0=3.3840 g · L^-1为例，实验数据及处理如表 5所示.

表4(Table 4)

表4 主、次反应区出口处B、C两点状态数据 Table 4 The B，C two state data at main and secondary reaction zone outlet

表4 主、次反应区出口处B、C两点状态数据 Table 4 The B，C two state data at main and secondary reaction zone outlet CS0/(g · L-1) B点 C点 CSB/(g · L-1) CSC/(g · L-1) FCH4/(g · min-1) FCO2/(g · min-1) VgC/(mL · h-1) 3.384 3.3113 3.282 2.13×10-3 2.40×10-4 210.09 3.5326 3.4584 3.4285 2.18×10-3 2.45×10-4 214.72 3.8729 3.7953 3.764 2.27×10-3 2.56×10-4 224.39

表5(Table 5)

表5 验证实验各点状态数据一览 Table 5 A view on states of B and C point in confirmatory experiment

表5 验证实验各点状态数据一览 Table 5 A view on states of B and C point in confirmatory experiment 时间/h B点 C点 CSB/(g · L-1) CSC/(g · L-1) 气体体积流量/(mL · h-1) CH4体积分率yCH4 CO2体积分率yCO2 4 3.2603 3.1982 254.87 76.67% 3.65% 8 2.9793 2.8638 242.19 79.61% 3.41% 平均值 3.0698 3.0310 248.03 78.14% 3.53%

CH₄体积流量v_CH4实为193.81 mL · h^-1，CH₄质量流量F_CH4实为2.045×10^-3 g · min^-1.同理，可计算CO₂体积流量、质量流量分别为8.75 mL · h^-1、2.54×10^-4 g · min^-1.误差按(实验值-计算值)/实验值×100%计算，可知本组B点基质浓度C_SB误差最大，为7.86%.

由表中数据可见：不同时刻反应器出口C点的气体体积流量v大于模拟计算时的体积流量v_gC，其原因可能是该验证实验时间不够长，系统中存在一定的N₂，同时反应过程中生成其它气体所致.其余两组数据计算依此类推，最后得到不同浓度时，C点出口处CH₄体积流量最大误差为9.67%，CO₂体积流量最大误差为8.53%，可见所建立的组合模型能够反映反应器实际操作状况.

根据以上反应过程模型，模拟了初始COD、温度变化对厌氧反应器状态的影响.

1)初始COD变化对反应器操作状态影响的模拟研究

操作参数为：废液体积流量v_0L=15 L · h^-1，温度35 ℃，改变COD初始值，B、C两处模拟数据如表 6所示.由表 6中数据可知，当温度一定时，随着基质初始浓度的增高，B、C两处基质浓度也随之升高，产物CH₄、CO₂的质量流量也升高.

表6(Table 6)

表6 反应器出口处B、C两点状态随基质浓度变化关系 Table 6 The relationship between the reactor outlet B，C two state changes with the concentration of substrate

表6 反应器出口处B、C两点状态随基质浓度变化关系 Table 6 The relationship between the reactor outlet B，C two state changes with the concentration of substrate CS0/(g · L-1) B点 C点 CSB/(g · L-1) FCH4/(g · min-1) FCO2/(g · min-1) vgB/(mL · h-1) CSC/(g · L-1) FCH4/(g · min-1) FCO2/(g · min-1) vgC/(mL · h-1) 3.0 2.9317 1.31×10-3 1.48×10-4 129.36 2.9043 2.01×10-3 2.26×10-4 198.21 3.2 3.1294 1.35×10-3 1.52×10-4 133.70 3.1010 2.07×10-3 2.34×10-4 204.79 3.4 3.3272 1.39×10-3 1.56×10-4 137.08 3.2979 2.13×10-3 2.40×10-4 210.54 3.6 3.5250 1.43×10-3 1.61×10-4 141.12 3.4948 2.20×10-3 2.47×10-4 216.88 3.8 3.7230 1.47×10-3 1.65×10-4 144.84 3.6920 2.25×10-3 2.54×10-4 222.46 4.0 3.9212 1.50×10-3 1.69×10-4 148.40 3.8894 2.31×10-3 2.59×10-4 227.94

2)反应温度变化对反应器操作状态影响模拟研究

操作参数为：废液体积流量v_0L=15 L · h^-1，COD初始值C_S0=3.40 g · L^-1，改变温度，B、C两处模拟数据如表 7所示.由表 7中数据可知，在基质初始浓度一定时，随着温度的升高，B、C两处的基质浓度不变，这主要是由于在该温度范围内，基质消耗活化能很低，以致消耗速率几乎不随温度变化而形成的结果；但产物CH₄、CO₂的质量流量随温度的升高而升高，同时CO₂的增加速率要快于CH₄，这主要是由于CO₂生成活化能高于CH₄，故其生成速率对温度更为敏感.

表7(Table 7)

表7 反应器出口处B、C两点状态随基质温度变化关系 Table 7 The relationship between the reactor outlet B，C two state changes with the temperature of matrix

表7 反应器出口处B、C两点状态随基质温度变化关系 Table 7 The relationship between the reactor outlet B，C two state changes with the temperature of matrix T/K B点 C点 CSB/(g · L-1) FCH4/(g · min-1) FCO2/(g · min-1) vgB/(mL · h-1) CSC/(g · L-1) FCH4/(g · min-1) FCO2/(g · min-1) vgC/(mL · h-1) 304.15 3.3272 1.21×10-3 1.32×10-4 117.38 3.2979 1.85×10-3 2.03×10-4 180.27 306.15 3.3272 1.29×10-3 1.44×10-4 126.91 3.2979 1.99×10-3 2.21×10-4 194.92 308.15 3.3272 1.39×10-3 1.56×10-4 137.08 3.2979 2.13×10-3 2.40×10-4 210.54 310.15 3.3272 1.46×10-3 1.70×10-4 144.86 3.2979 2.23×10-3 2.61×10-4 221.68 312.15 3.3272 1.59×10-3 1.84×10-4 159.49 3.2979 2.45×10-3 2.83×10-4 244.96

1)通过对反应器主反应区停留时间分布的考察，主反应区流型可以用一平推流区与两个并联的全混区进行串联组成，各区所占体积分率随反应过程中所产气体流量不同而不同，随着气体流量的增大，平推流区所占体积分率减小.

2)次反应区由平推流区和全混区并联组合而成，整个厌氧反应器由主反应区和次反应区串联组合而成，根据实验数据可知，C点出口处CH₄体积流量最大误差为9.67%，在可允许的范围内，说明该组合模型较好地模拟了反应器的实际状态.

3)采取单因素分析方法，模拟了基质初始COD和温度变化对反应器出口状态的影响，发现产物CH₄、CO₂的质量流量随基质浓度、温度的升高而升高；在304.15~312.15 K范围内，温度的变化对出口处基质浓度的影响不大.