基于 ASM1 模型的污水处理工艺认知与代码实现

生化池（厌氧 + 缺氧 + 好氧）并没有把污染物真的从水里'拿出来'，而是把它们'转移'到了细菌体内。

第二部分：ASM1 理论篇

如果说 AAO 工艺是污水处理厂的'硬件架构'，那么 ASM1 模型就是运行在其中的'软件内核'。它用一组微分方程 (ODEs) 精确描述了：什么细菌 (Who)，在什么环境 (Where)，以多快的速度 (How Fast)，吃掉了多少污染物 (How Much)。

ASM1 模型的状态变量（13 个）

ASM1 模型的生化动力学反应（8 个）

化学计量常数（5 个）

'化学计量常数'理解成反应里各物质的'定量配方'——就像做饭时'1 份米配 2 份水'的比例，它规定了生物反应中'反应物消耗多少、产物生成多少'的定量关系。

动力学参数（15 个）

'动力学'研究的是'反应速度规律'，这些参数就是用来量化微生物反应'有多快'的——比如细菌繁殖多快、底物降解多快、环境变化对反应速度的影响有多大。就像给生物反应装了个'速度仪表盘'，所以叫'动力学参数'。

反应过程图解

反应速率

异养菌的好氧生长

异养菌在好氧环境下，根据溶解性底物和溶解氧的浓度，按照'底物越多、氧气越充足，生长速度越快'的规律繁殖，公式里的μ H 是繁殖速度上限，K S 和 K O,H 是半饱和浓度（到这个浓度时生长速度减半），最后乘以异养菌浓度（X B,H）得到生长速率。

异养菌的缺氧生长

异养菌在缺氧环境下，用**硝酸盐（S NO）**代替氧气，同样根据底物（S S）浓度繁殖，但此时氧气浓度（S O）越低、硝酸盐浓度越高，生长速度越快，η g 是缺氧生长的效率修正系数，最后乘以异养菌浓度得到生长速率。

自养菌的好氧生长

自养菌（硝化菌）在好氧环境下，根据**氨氮（S NH）**和溶解氧（S O）的浓度繁殖，μ H 是繁殖速度上限，K NH 和 K O,A 是半饱和浓度，乘以自养菌浓度（X B,A）得到生长速率。

异养菌的衰减

异养菌自然死亡的速率，由衰减系数（b H）和异养菌浓度（X B,H）决定，衰减系数越大、菌浓度越高，死亡速度越快。

自养菌的衰减

自养菌自然死亡的速率，由衰减系数和自养菌浓度决定，逻辑同异养菌衰减。

溶解性有机氮氨化

溶解性有机氮（S ND）被异养菌分解成氨氮的速率，由氨化速率常数（k a）、有机氮浓度和异养菌浓度决定，速率常数越大、有机氮和菌浓度越高，氨化越快。

慢速可生物降解有机物的水解

颗粒性有机物（X S）被分解成溶解性有机物的速率，受底物浓度（X S 与 X B,H 的比值）、溶解氧（S O）、硝酸盐（S NO）浓度影响，好氧时水解快，缺氧时需乘以修正系数（η h），最后乘以异养菌浓度得到水解速率。

颗粒性可生物降解有机氮的水解

颗粒性有机氮（X ND）的水解速率，是慢速有机物水解速率（ρ 7）乘以有机氮与有机物的比例（X ND/X S），逻辑是有机氮水解规律和有机物水解一致，只是比例不同。

各组分转化速率

• SI（溶解性惰性有机物）： 不变，因为细菌不吃它。
• SS（溶解性可生化底物）： 异养菌好氧、缺氧吃的量，减去水解生成的量，就是它的变化。
• XI（颗粒惰性有机物）： 不变，细菌也不吃它。
• XS（颗粒可生化底物）： 异养菌、自养菌死了产生的颗粒底物，减去水解成溶解性底物的量，就是它的变化。
• XBH（活性异养菌）： 异养菌好氧、缺氧繁殖的量，减去自己死亡的量，就是它的变化。
• XBA（活性自养菌）： 自养菌繁殖的量，减去自己死亡的量，就是它的变化。
• XP（微生物衰减污泥颗粒）： 异养菌、自养菌死了变成的惰性污泥，死得越多它越多。
• SO（溶解氧）： 异养菌好氧、自养菌硝化消耗的氧气，消耗越多它越少。
• SNO（硝酸盐氮）： 异养菌缺氧消耗的量，减去自养菌硝化生成的量，就是它的变化。
• SNH（氨氮）： 异养菌、自养菌消耗的量，加上有机氮氨化生成的量，就是它的变化。
• SND（溶解性可生物降解有机氮）： 氨化消耗的量，减去有机氮水解生成的量，就是它的变化。
• XND（颗粒性可生物降解有机氮）： 细菌死了产生的颗粒有机氮，减去水解成溶解性有机氮的量，就是它的变化。
• S ALK（碱度）： 异养菌、自养菌反应消耗或生成的碱度，加上氨化生成的碱度，就是它的变化。

生化反应池模型的建立

第三部分：单个池子单个组分的手动计算和代码实现

以氨氮为例：

假设我们有一个曝气池（CSTR 完全混合反应器），我们需要计算当前时刻氨氮的反应速率和浓度变化。

水质组分（状态变量）当前浓度：

动力学参数 (ASM1 典型值)：

第一步：计算动力学过程速率

第二步：计算氨氮的反应速率

第三步：计算质量平衡

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# --------------------------
# 解决中文显示问题
# --------------------------
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC", "Microsoft YaHei"]
plt.rcParams["axes.unicode_minus"] = False

# --------------------------
# 1. 参数设置（ASM1 标准，20℃）
# --------------------------
params = {
    # --- 动力学参数 (单位: /d) ---
    'mu_A': 0.5,          # 自养菌最大比生长速率 (1/d)
    'k_a': 0.05,          # 氨化速率 (m³/(g COD·d))
    'b_A': 0.05,          # 自养菌衰减系数 (1/d)
    # --- 半饱和常数 (单位: mg/L) ---
    'K_NH': 1.0,          # 氨氮半饱和常数
    'K_OA': 0.4,          # 自养菌氧半饱和常数
    # --- 化学计量参数 ---
    'Y_A': 0.24,          # 自养菌产率 (g VSS/g N)
    'i_XB': 0.086,        # 生物体含氮量 (g N/g VSS)
    # --- 水力与环境参数 ---
    'q_in': 100.0,        # 进水流量 (m³/h)
    'q_out': 100.0,       # 出水流量 (m³/h)
    'V': 1000.0,          # 反应池容积 (m³)
    # --- 进水条件 ---
    'S_NH_in': 30.0       # 进水氨氮浓度 (mg N/L)
}

# --------------------------
# 2. 辅助函数：计算反应速率
# --------------------------
def calc_rates_nh(S_NH, env_vars, params):
    """
    计算与氨氮相关的反应过程速率 (Process Rates)
    env_vars: 包含 [S_O, S_ND, X_BA, X_BH] 的环境常量
    """
    # 解包环境常量
    S_O = env_vars['S_O']      # 溶解氧
    S_ND = env_vars['S_ND']    # 溶解性有机氮
    X_BA = env_vars['X_BA']    # 自养菌 (硝化菌)
    X_BH = env_vars['X_BH']    # 异养菌 (用于氨化)

    # --- j=3: 自养菌好氧生长 (硝化作用) ---
    # 公式：rho3 = mu_A * (S_NH/(K_NH+S_NH)) * (S_O/(K_OA+S_O)) * X_BA
    # 注意：params['mu_A'] 是 /d，需除以 24 转换为 /h
    monod_NH = S_NH / (params['K_NH'] + S_NH)
    monod_O = S_O / (params['K_OA'] + S_O)
    rho3 = (params['mu_A'] / 24.0) * monod_NH * monod_O * X_BA

    # --- j=6: 溶解性有机氮氨化 ---
    # 公式：rho6 = k_a * S_ND * X_BH
    rho6 = (params['k_a'] / 24.0) * S_ND * X_BH

    # --- 计算氨氮的总反应速率 r_NH (mg N/(L·h)) ---
    # r_NH = -消耗 + 生成
    # 消耗 = (i_XB + 1/Y_A) * rho3 (细菌组成氮 + 氧化产能消耗的氮)
    # 生成 = rho6 (有机氮水解)
    consumption = (params['i_XB'] + 1.0 / params['Y_A']) * rho3
    generation = rho6
    r_NH = -consumption + generation
    return rho3, rho6, r_NH

# --------------------------
# 3. 氨氮动态模型
# --------------------------
def nh_dynamics(t, y, env_vars, params):
    """氨氮浓度随时间变化的微分方程"""
    S_NH = y[0]
    # 边界保护：浓度不能为负
    if S_NH < 0:
        S_NH = 0

    # 1. 计算生化反应项 (r_NH)
    _, _, r_NH = calc_rates_nh(S_NH, env_vars, params)

    # 2. 计算物理输运项 (Transport)
    # 公式：(Qin * Sin - Qout * S) / V
    transport = (params['q_in'] * params['S_NH_in'] - params['q_out'] * S_NH) / params['V']

    # 3. 总变化率 dS_NH/dt
    dSNH_dt = transport + r_NH
    return [dSNH_dt]

# --------------------------
# 4. 模拟设置与运行
# --------------------------
# 设定环境常量 (假设在短时间内相对稳定)
env_vars = {
    'S_O': 2.0,           # 溶解氧 (mg/L)
    'S_ND': 5.0,          # 有机氮 (mg/L)
    'X_BA': 100.0,        # 自养菌浓度 (mg/L) - 恒定
    'X_BH': 2000.0        # 异养菌浓度 (mg/L) - 恒定
}

# 初始状态
S_NH_init = 20.0       # 初始氨氮 (mg/L)

# 模拟时间 (模拟 24 小时)
t_span = (0, 24)
t_eval = np.linspace(0, 24, 200)

# 求解微分方程
sol = solve_ivp(
    fun=lambda t, y: nh_dynamics(t, y, env_vars, params),
    t_span=t_span,
    y0=[S_NH_init],
    t_eval=t_eval,
    method='RK45'
)

# --------------------------
# 5. 结果可视化
# --------------------------
plt.figure(figsize=(10, 6))
# 氨氮变化曲线
plt.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='氨氮浓度 (S_NH)')
plt.axhline(y=params['S_NH_in'], color='r', linestyle='--', alpha=0.6, label='进水浓度 (30 mg/L)')

# 标注关键信息
plt.xlabel('时间 (小时)', fontsize=12)
plt.ylabel('氨氮浓度 (mg N/L)', fontsize=12)
plt.title(f'好氧池氨氮动态模拟 (DO={env_vars["S_O"]} mg/L)', fontsize=14)
plt.grid(alpha=0.3)
plt.legend(fontsize=10)

# 标注最终平衡浓度
final_conc = sol.y[0][-1]
plt.text(0.6, 0.5, f'最终平衡浓度：{final_conc:.2f} mg/L', transform=plt.gca().transAxes, fontsize=10, bbox=dict(facecolor='white', alpha=0.8, edgecolor='gray'))
plt.tight_layout()
plt.show()

# --------------------------
# 6. 关键速率输出 (t=0 时刻)
# --------------------------
rho3_init, rho6_init, r_NH_init = calc_rates_nh(S_NH_init, env_vars, params)
print("--- 初始时刻 (t=0) 速率分析 ---")
print(f"1. 硝化反应速率 (rho3): {rho3_init:.4f} mg VSS/(L·h)")
print(f" -> 导致氨氮消耗：{(params['i_XB'] + 1/params['Y_A']) * rho3_init:.4f} mg N/(L·h)")
print(f"2. 氨化反应速率 (rho6): {rho6_init:.4f} mg N/(L·h)")
print(f"3. 净生化反应速率 (r_NH): {r_NH_init:.4f} mg N/(L·h) (正值代表生成，负值代表消耗)")
print(f"4. 物理稀释/积累速率：{(params['q_in']*(params['S_NH_in']-S_NH_init))/params['V']:.4f} mg N/(L·h)")