4.2.4　厌氧消化过程的理论模型

4.2.4.1　污泥厌氧消化产甲烷动力学原理

（1）概述

污泥的厌氧消化过程是污泥所参与的、以微生物（包括水解菌、产酸菌及产甲烷菌等这些微小的生物体）为主体所催化的生化反应过程。厌氧消化过程有以下主要特征。

1）微生物是反应过程的主体　首先，微生物是反应过程的催化剂，它摄取了原料中的养分，通过微生物内特定的酶系统进行复杂的生化反应，把污泥转化为有用的产品（生物能）；同时，它又如同微小的反应器，原料中的反应物通过渗透作用经由微生物细胞壁和细胞膜进入到微生物体内，在酶的作用下进行催化反应，把反应物转化为产物，接着产物又被释放出来。

2）厌氧消化过程的本质是复杂的酶催化反应体系　这一复杂的反应体系就是通过微生物的代谢作用来维持微生物内物质和能量的自身平衡的。一方面，微生物通过同化作用将外界摄取的营养物质转化为自身的组成物质；另一方面，又将微生物内的组成物质不断地分解排出，这称为异化作用。

3）厌氧消化过程是动态的　即在反应进行的同时，微生物也得到生长，其形态、组成、活性都处在一个动态变化的过程。从微生物的组成分析，它包含蛋白质、脂肪、碳水化合物等，这些成分的含量也随着环境的变化而变化。微生物能通过代谢机制进行定量调节以适应外界环境的变化^［2］。

以上这些因素，再加上实际操作中有机类物质的来源不同，造成了厌氧消化过程描述和控制的复杂性。

针对不同的厌氧消化设备和消化工艺，在做出必要的简化和合理假设的条件下，依据经验、理论和对试验结果的分析，运用适当的数学工具能定量描述整个厌氧消化过程各参数的特性变化。

模型建立后，模型参数估值成为解决问题的关键。参数估值的合理性和准确性关系到模型的可靠性和实用性。本章将简单介绍厌氧消化的一些基本理论，然后重点介绍厌氧消化过程建模和模型参数求解的常用方法。

（2）厌氧消化动力学原理

微生物的生长、繁殖代谢是一个复杂的生物化学过程，既包括微生物体内的生化反应，也包括微生物体内和体外的物质交换，同时还包括微生物体外的物质传递和反应。每个微生物都经历着生长、成熟和衰老的过程，同时伴随有变异和退化。该体系具有多相、多组分、非线性的特点，要对这样一个反应过程进行准确的描述几乎是不可能的。因此，建立过程的数学模型，首先要进行合理的简化。模型假设有以下几点：a.微生物反应动力学是对微生物群体动力学行为的描述，而不是对单一微生物进行描述；b.不考虑微生物之间的差别，认为其内部结构是一致的；c.微生物为单组分，有确定的化学方程式描述；d.发酵罐内基质（污泥）的物性认为是均匀一致的。

1）厌氧消化过程计量学方程及产率系数　反应计量学是对反应物的组成和反应转化程度的数量化研究。对厌氧消化过程而言，由于众多组分参与反应和代谢途径的错综复杂，同时在微生物生长的过程中还伴随着代谢产物生成的反应，因此，需要通过一些特殊的方式加以简化处理。

为了表示消化过程中各物质和各组分之间的数量关系，最常用的方法是建立各元素之间的原子衡算方程式。建立原子衡算方程式首先要确定微生物的元素组成及其分子式。不同的微生物，其组成是不一样的，即便是同一微生物，所处生长阶段不同，其组成也不一样。在实际中，常采用平均微生物组成，并给出了经验分子式。一般将微生物的分子式定义为C_αH_βO_δN_γ，忽略像P、S和灰分等微量元素或成分。为此，周少奇^［10］根据生化反应电子流守恒原理，建立了有机类垃圾厌氧消化过程的生化计量学方程。根据计量方程式，我们就可以对初始有机物、中间产物、终产物以及微生物之间的关系进行定量描述。该计量方程式将蛋白质、死菌体和其他动植物残渣归结为含氮有机物（C_aH_bO_cN_d），将淀粉、纤维素和木质素归结为烃类化合物（C_αH_βO_γ），中间产物为乙酸（CH₃COOH）、丙酸（C₂H₅COOH）和丁酸（C₃H₇COOH），终产物为H₂、CH₄和CO₂。在应用计量方程定量描述之前，需要阐明一下产率系数的概念。

产率系数分宏观产率系数和理论产率系数。假定发酵过程中所消耗的基质总量为ΔS_T，用于微生物生长的基质数量为ΔS_G，用于生成代谢产物的量为ΔS_R。若定义Y_x/S=Δx/-ΔS_T=Δx/-（ΔS_G+ΔS_R），此时，求得的对基质的微生物产率系数称为宏观产率系数；若定义Y_x/S=Δx/-ΔS_R，此时，求得的对基质的微生物产率系数为理论产率系数。在实际中，用于微生物生长的量很难准确测量，而这个值和用于生成代谢产物的量相比又很小。因此，通常采用理论产率系数的简化式（忽略ΔS_G）来求解产率系数。我们仅以水解阶段的计量方程式为例，说明各参数之间的定量关系。

C₆H₁₀O₅·nNH₃y_eC₆H₁₀O₅+（1-y_e）C₆H₁₀O₅·mNH₃+［n-m（1-y_e）］NH₃

此计量方程式是污泥（C₆H₁₀O₅·nNH₃）水解为可溶有机物（C₆H₁₀O₅）和不可溶有机物（C₆H₁₀O₅·mNH₃）以及氨氮的表达式。通过此方程式，可以计算出可溶性有机物的产率系数为1/y_e；也可以计算不可溶有机物的产率系数为1/（1-y_e）；氨氮的产率系数为1/［n-m（1-y_e）］，其中y_e为污泥中的可降解系数。

2）污泥厌氧消化微生物生长动力学模型　厌氧消化过程，包括微生物的生长、基质的消耗和代谢产物（氢气、甲烷）的生成，要定量描述厌氧消化过程的速率，显然微生物生长动力学是其核心。

① 微生物比生长速率和浓度表达式。不同的污泥厌氧消化工艺有一个共同的特点，就是都必须经过水解、酸化过程才能将污泥降解并产生生物气体。对于水解过程，用Monod模型来表述；对于接下来的几个消化过程，常采用形式较复杂的Haldane模型来阐述厌氧菌的生长，即对于仅受基质限制的微生物，其生长动力学模型通常采用Monod模型来描述；对于同时受基质限制和基质抑制的微生物，采用反竞争性抑制模型来表达。图4-13所示为Monod模型和反竞争性抑制模型的比生长速率对比。

从图4-13可以看出，在反竞争抑制模型中，当基质浓度低时，微生物比生长速率（μ）随基质浓度的提高而增大，并达到最大值；当基质浓度继续提高时，μ反而会降低。底物的抑制作用影响因子可以采用以下方程形式来表述：

　　（4-31）

式中，K_I、C分别表示底物（或产物）对厌氧菌的抑制系数和底物（产物）浓度。

水解菌的比生长速率表达式如下：

　　（4-32）

式中，μ_h为水解菌比生长速率，d^-1；μ_hmax为水解菌最大比生长速率，d^-1；C_S为挥发性固体浓度，g/L；K_SS为挥发性固体的半速率系数，g/L。

产酸菌的比生长速率表达式如下：

　　（4-33）

式中，μ₂为产酸菌的比生长速率，d^-1；μ_2，max为产酸菌的最大比生长速率，d^-1；S₂为挥发性脂肪酸浓度，g/L；K_S，2为挥发性脂肪酸的半速率系数，g/L；为产酸菌生长抑制系数，g/L。

产甲烷菌的比生长速率表达式如下：

　　（4-34）

式中，μ₃为产甲烷菌或产氢菌的比生长速率，d^-1；μ_3，max为产甲烷菌或产氢菌的最大比生长速率，d^-1；S₃为乙酸浓度，g/L；K_S，3为乙酸的半速率系数，g/L；为产甲烷菌或产氢菌的生长抑制系数，为氨氮的抑制系数，g/L。

水解菌浓度变化速率为：

　　（4-35）

式中，K_dh为水解菌的比死亡速率，d^-1；X_h为水解菌浓度，g/L。

产酸菌浓度变化速率为：

　　（4-36）

式中，K_da为产酸菌衰亡速率系数，d^-1；X_a为产酸菌浓度，g/L。

产甲烷菌和产氢菌浓率变化速率用相同表达式表述如下：

　　（4-37）

式中，K_dm为不同过程相应厌氧菌的衰减常数。

② 系统pH值和温度的影响因子。在污泥厌氧消化过程中，维持微生物生长需要适宜的环境条件，其中pH值和温度是最重要的两个因素。

根据温度，厌氧消化常分为常温消化、中温消化和高温消化3种。当温度偏低时，微生物的生长缓慢；随着温度升高，生长速率变大；当温度超过一定范围时，可能导致微生物热死亡。通常，温度对微生物生长速率的影响因子可以通过下面的关系式来表达：

　　（4-38）

　　（4-39）

式中，μ₀为某一温度下的比生长速率，d^-1；为细菌死亡速率，d^-1；E_a为细菌活化能，kJ/kmol。

随着厌氧消化过程的进行，生成有机酸和氨氮，系统pH值会产生变化。正常情况下，微生物在生长阶段，pH值有上升或下降；在产物生成阶段，pH值趋于稳定；在细胞自溶阶段，pH值会上升。微生物本身具有调节pH值的能力，但当外界条件变化剧烈时，微生物失去调节能力，pH值就会发生波动。

pH值对微生物生长速率的影响因子可用式（4-40）描述：

　　（4-40）

③ 氨氮的抑制模式及影响因子。在厌氧消化过程的研究中，氨氮抑制产气的过程普遍存在四阶段模式，该四阶段模式如下：

　　（4-41）

　　（4-42）

　　（4-43）

　　（4-44）

随着游离氨浓度的增加，f（抑制因子）呈现不同形式的下降。当游离氨的浓度为1.1mg/L时，f是稳定的；当浓度从1.1mg/L增加到1.16mg/L时，f从1.0降到0.67；当浓度持续增加到1.34mg/L时，f保持稳定下降。非离子化NH₃的浓度主要取决于3个因素，即总氨氮浓度、温度和pH值。pH值对游离氨所占的比例有很大影响，当pH值上升到8时，游离氨仅占总氨氮的0.1%左右。

3）污泥厌氧消化中底物消耗动力学　污泥的消耗速率可以通过微生物的比生长速率和产率系数联系起来，污泥的降解速率可以通过式（4-45）～式（4-47）表达：

　　（4-45）

式中，μ_h为水解菌的比生长速率，d^-1；X_h为水解菌的浓度，g/L；Y_h为S_h的产率系数（水解菌/基质）。

　　（4-46）

式中，μ_a为产酸菌的比生长速率，d^-1；S_a为产酸菌的基质浓度，g/L；X_a为产酸菌的浓度，g/L；Y_a为S_a的降解系数（产酸菌/基质）；Y_vh为S_a的产率（产酸菌/基质）。

　　（4-47）

式中，X_m为产甲烷菌质量浓度，g/L；μ_m为产甲烷菌的比生长系数，d^-1；A为系统乙酸总质量浓度，g/L。

4）产物生成动力学　污泥厌氧消化过程的甲烷生成速率表达式如下：

　　（4-48）

式中，V_mmax为标准状态下（0℃，1.1325×10⁵Pa）CH₄的体积。

氢气生成速率表达式如下：

　　（4-49）

式中，V_hymax为标准状态下（0℃，1.1325×10⁵Pa）H₂的体积。

根据Keshtkar的理论，在厌氧过程，只有水解阶段生成氨氮，但整个过程都会消耗氨氮，因此氨氮浓度随时间的变化可表示为：

（4-50）

当μ_i-K_di≤0，Y_N为0。

　　（4-51）

式中，为NH₃的摩尔质量，17g/mol；c（N）表示氨离子浓度，g/L。

污泥厌氧消化过程伴随着生物反应热现象，是致使系统温度改变的重要因素之一，它的存在造成保温过程中温控的复杂性。在发酵系统保温装置的设计中，为了简化设计，常常忽略生物反应热的影响，但是，由于温度是影响污泥厌氧消化过程的最重要的调控参数之一，所以在实际发酵中，通常要求温差控制在ΔT≤1℃。因此，对污泥厌氧消化过程中生物反应热的计算需引起一定的重视。

微生物参与的发酵过程的产热速率，可用与生成产物相似的方程来进行处理。若采用比产热速率概念，则有式（4-52）：

　　（4-52）

式中，q_Hv为比产热速率，W/（d·g）；Y_X/Hv为微生物热产率系数，W/g；Y_P/Hv为产物热产率系数，W/g；m_Hv为细菌维持能，W/（d·g）。

5）发酵系统内的相平衡　厌氧消化系统中的二氧化碳、甲烷等实际气体可视作理想溶体，描述其气相和溶解平衡的经验性公式可由亨利定律来表达。

亨利定律表述为稀薄气体的蒸气分压正比于其液相浓度。因此，我们可以用式（4-53）来描述二氧化碳、甲烷气体在液相的浓度和相应蒸气分压的关系：

　　（4-53）

式中，H_c、H_m分别为二氧化碳和甲烷的亨利常数；P_c、P_m分别为二氧化碳和甲烷的蒸气分压。

在厌氧消化过程中，伴随有甲烷的产生和逸出，还存在二氧化碳的消耗。它们的蒸气分压随时间的变化规律用式（4-54）来表示：

　　（4-54）

式中，F_t为在厌氧消化过程中产生的气体压力，忽略水蒸气分压力影响；R为气体常数；T为厌氧消化温度；V_g为消化罐内液体容积。

6）发酵系统内的液相电离平衡　厌氧消化系统液相内存在离子的电离平衡，例如二氧化碳和水可以形成弱碱根离子HC和氢根离子H⁺，氨根离子电离可以形成氨氮和氢离子。这些物质的电离，将会影响系统的pH值，进而影响整个发酵过程的进行。简化后的电离平衡总方程式如下：

　　（4-55）

式中，c（Ac^-）为乙酸根离子浓度。

相关组分的电离方程及平衡常数如式（4-56）～式（4-59）所示：

　　　　　　（4-56）

　　　　HCC+H⁺ 　　k₂=　　（4-57）

　　　　　HAcHAc^-+H⁺ 　　k₃=　　（4-58）

　　　　NNH₃+H⁺ 　　k₄=　　（4-59）

式中，k₁～k₄为对应组分的电离平衡常数。

污泥厌氧消化工艺可分为单相和两相工艺。在两相工艺中，实现相分离的途径有化学法、物理法和动力学控制方法。目前常用的是动力学控制法。该方法利用产酸菌和产甲烷菌生长速率上的差异，控制两个反应器的有机负荷率（OLR）、水力停留时间等参数来实现相的分离。这种相的分离不是绝对的，只是在发酵的不同阶段，系统中的优势菌群不同。

产酸相和产甲烷相的分离是通过控制两相的水力停留时间来实现的。定义一个参数R，该参数跟进料间隔时间t_R和水力停留时间t_HR相关。令t_HR=t_R （R+1），从而得到：R=t_HR/t_R-1，该参数将被应用于确定新的循环状态变量的求解中。

对半连续操作而言，假设物料在发酵罐内停留时间t_R，到达一定时间后，体积为V_s的发酵液排出系统，随后，同体积的新鲜物料加入反应罐，最后充分混合。新循环状态变量的确定可以通过以下方程来计算，用C_Af表示进料状态，用C_Ae表示出料状态，那么对于下一个t_R而言，新循环状态变量的确定方式如下：

　　（4-60）

4.2.4.2　污泥厌氧消化模型的参数求解

（1）模型参数求解的方法

模型参数的估算关系到数值模拟的准确性和可靠性，是模型能否推广应用的重要影响因素，因此，需要对模型参数进行准确求解和优化。对实际的发酵过程，其动力学方程参数必须辅以实验方法来确定。

动力学参数的求解方法很多，常用的有微分法、积分法、回归法、神经网络法以及遗传算法等。积分法受方程积分难易程度影响，只应用于最简单的动力学形式。微分法是依据不同实验条件下的反应速率，直接由反应速率方程来估计参数值，结合线性化作图更方便。随着计算机技术的进步发展起来的回归法（常用的是最小二乘法），因为不受参数数目多少的限制，又不受动力学方程的形式和实验方法的约束，具有普遍的适用性，其缺点是往往只能得到局部最优解。遗传算法可以计算得到目标函数的全局最优解，且不受初始值的限制，其除了具有非线性回归方法的优点外，还具有更大的优势。本节中动力学参数的估值将采用微分法和遗传算法，对便于通过实验数据求取的动力学参数采用微分法，其余的参数通过遗传算法求解。因此，重点对以上两种方法的原理和步骤进行介绍^［2］。

（2）微分法求解动力学参数

对于幂函数型动力学，例如两边取对数，则有根据实验数据，以对作图，可确定参数和n的值。图4-14为微分法求解动力学参数的示意图，函数f（C）代表模型中假设的函数关系。

下面以Monod模型参数的求解为例，说明微分法的具体用法。微生物的比生长速率可在一定时间范围内直接用图解微分法求解：

　　（4-61）

式中，为测定时间区间Δt内微生物浓度的平均值。根据Monod模型的μ-S关系式可以得到式（4-62）：

　　（4-62）

以对作图，可确定动力学参数。从试验数据的获取到动力学参数的求解，整个过程表示在图4-15中。

当然，除了图4-15提到的L-B作图法外，根据不同的线性化方法，可以采用E-H作图法及Langmuir作图法等以外的求解方法，它们的计算过程类似，但计算的精度有所不同，具体求解过程这里不再详述，可参见文献。

（3）遗传算法在求解动力学参数中的应用

1）遗传算法简介　现存的生命是经过漫长的岁月，在适应了地球上各种各样的环境条件，逐渐进化和发展起来的。更确切地说，正是由于突然变异和种群交配，才使得更适合自然环境的种群和个体得以生存下来。近几十年来，随着计算机技术的不断发展，研究者们开始利用计算机来仿真生物的遗传和进化过程，后来即发展演变成非常著名的“遗传算法”。遗传算法（genetic algorithm，GA）就是通过模仿生物进化过程而开发出来的一种概率探索、自我适应、自我学习和最优化的方法。现在，遗传算法已经广泛应用于许多领域，如过程模型参数的确定、过程的最优化控制、系统工程中的优化求解等。

遗传算法最主要的应用就是对特定的过程和系统进行优化。该法具有计算精度高、收敛速度快等优点。使用遗传算法得到的全局最优解一般不受初始条件的影响和限制，特别适合具有复杂和高度非线性化特性的生物过程。遗传算法在生物过程中的应用主要包括过程的最优化控制——求解最优化轨道（最优温度轨道、最优pH值轨道等）、优化发酵培养基、确定生物反应模型的参数等。

2）遗传算法在求解动力学参数中的应用　遗传算法的计算方法多种多样，这里，以单纯的遗传算法（simplegeneticahhorithm，SGA）来介绍其求解动力学参数的具体方法。生物体的特征是由基本构成单位——基因（gene）所组成的染色体又称个体（chromosome）来加以显示的。而生物体本身又是由多个这样的染色体所形成的集团或种群（population）所构成的。每个染色体由多个基因构成，每一个基因所处的位置被称为基因座（10cus）。染色体的数量则称为种群数（populationsize）。所形成的生物体要通过其在外界环境下生存的适合度（fitness）来进行评价、筛选和生存竞争，即按照“优胜劣汰、适者生存”的原则来实现种群的选择。另外，染色体中的一部分基因还要经过突然变异（mutation）而发生变化。不断地重复上面的操作保证了生物个体进化的不断进行，而操作循环的次数被称为传代数（generation）。下面详细介绍遗传算法求解模型参数的具体步骤。

定义如下的目标函数：

　　（4-63）

式中，x_j和x'_j分别为模型参数的计算值和对应发酵时刻的实测值。

适合度和选择概率的定义式分别如下：

　　（4-64）

　　（4-65）

遗传算法具体的计算步骤及计算规则如下。

① 基因生成。规定待优化模型参数的最大值，并将其4等分，由计算机随机产生二进制文字序列构成的初代染色体种群M个。每个染色体就代表一整套完整的模型参数向量。在本节中，染色体中基因序列长度均为8个字节，长度的选取原则上主要取决于求解变量的计算精度、总的计算量和收敛速度。

② 染色体各基因序列的解码化。依据图4-16所示的方式，将各染色体的各段基因序列由二进制文字序列转变为十进制的数字序列。

③ 计算目标函数值。将解码后的染色体代入建立的模型中迭代求解，然后代入目标函数公式，计算目标函数值。

④ 求解适合度。根据上一步计算得到的J，代入上面适合度的定义式，计算染色体对应的适合度。

⑤ 计算选择概率。计算出染色体适合度之后，根据选择概率的定义式，计算选择概率。

⑥ 选择染色体种群——Roulett圆盘选择法。在求出所有M个染色体的选择概率之后，利用Roulett圆盘选择法来确定进入下一代的染色体种群，见图4-17。

在圆盘上标记上各染色体的选择概率S以及对应的面积，将圆盘按逆时针方向旋转。将圆盘旋转M次，可以得到M个进入下一代的染色体种群。很显然，适合度高的染色体进入下一代的可能性较大，适合度较低的可能性较小，但也不是绝对没有可能。这种种群选择的方法充分体现了 “优胜劣汰，适者生存”的原则。在进行个体交叉和突然变异等“遗传操作”时，个体很容易受到破坏而消失。

为了防止适合度高的染色体i在下一步的遗传操作中被破坏，常采用E.litist保存选择方法对种群进行选择。该选择法是将种群中适合度最高的个体保留，让其无条件地进入下一代的种群中去。在一般情况下，遗传算法中常常结合这两种方法来对种群进行选择。

⑦ 交叉操作。从所选定的M个染色体种群中任选两个个体作为亲本交配个体，然后在其中一个亲本个体的基因序列中随机地选择一点或多点作为交叉操作的交叉点（cross-overpoint）。按照图4-18的方式，在交叉点处亲本个体的基因序列进行交换，当各随机数小于预先设定的交叉概率P时，交叉操作进行，否则交叉操作不进行。同样，剩下的M-2个个体也按照同样的步骤进行交叉操作。交叉方式主要有“一点交叉”（one-pointcrossover）或“多点交叉”（multi-pointcrossover）两种。

⑧ 突然变异操作。从优化计算的角度来看，交叉操作不可能寻找到复杂函数和系统的全局最优解，反而可能会使优化求解陷入局部最优解之中。突然变异操作可避免“交叉操作”导致的“近亲繁殖”的不良后果。突然变异操作按照图4-19进行，在染色体的某一基因序列上，用与其对立的基因序列与原有的基因序列进行置换，从而产生了完全不同于原有种群的新的个体。正是由于突然变异操作的存在，保证了利用遗传算法进行复杂函数或系统的优化计算时可以寻找到全局的最优解。从整个种群中任选一个个体作为亲本变异个体，然后在该亲本个体中随机地选取一段基因序列，用与其对立的基因序列进行置换，再由计算机产生［0，1］范围的随机数，当该随机数大于预先设定的突然变异概率的P_m时，该变异操作不进行；反之，进行上面的变异操作。剩下的M-1个染色体按照同样的步骤，由计算机来完成操作。为了不破坏由“交叉操作”产生的优良个体，一般情况下P_m应该小于P_c，如果P_m选择过大，优化计算就相当于随机搜索，计算很难收敛。

⑨ 重复②～⑧步骤，直到达到预先设定的最大传代数。计算中，种群数M=20，最大传代数1000，交叉方法为两点交叉，交叉概率P_c=1.0，种群的选择方式为Roulett圆盘法+Elitist保存选择法，突然变异概率P_m=0.07。图4-20给出了遗传算法求解模型参数的简单示意图，可以通过该图对遗传算法有一个总体的认识和了解。