多场耦合优化-主题054-水处理与环境工程仿真

kkchenkx

148人浏览 · 2026-03-18 21:12:39

kkchenkx · 2026-03-18 21:12:39 发布

主题054：水处理与环境工程仿真

一、水处理与环境工程概述

水资源是人类生存和发展的重要基础资源。随着工业化进程加快和人口增长，水环境污染问题日益严峻，水处理技术的重要性愈发凸显。环境工程仿真通过数学建模和数值计算，为水处理工艺设计、优化和运行管理提供科学依据，是实现水资源可持续利用的重要技术手段。

1.1 水处理的分类与目标

给水处理：

原水取水与预处理
混凝沉淀
过滤消毒
深度处理（膜分离、活性炭吸附）
出水水质达标

污水处理：

物理处理（格栅、沉砂、沉淀）
化学处理（中和、氧化还原、混凝）
生物处理（活性污泥法、生物膜法）
深度处理（脱氮除磷、膜生物反应器）
污泥处理与处置

工业废水处理：

重金属去除
有机污染物降解
高盐废水处理
零排放技术

水环境修复：

河流生态修复
湖泊富营养化控制
地下水污染修复
人工湿地构建

1.2 环境工程仿真的意义

工艺设计优化：

反应器尺寸确定
工艺参数优化
能耗分析
成本估算

运行管理支持：

水质预测
故障诊断
自动控制
应急预案

环境风险评估：

污染物迁移转化
生态影响评价
健康风险分析
环境容量计算

二、反应器水力学

2.1 反应器类型与流态

理想反应器模型：

完全混合反应器（CSTR）：
$\frac{dC}{dt} = Q(C_{in} - C) + rV$

其中， $V$ 为反应器体积， $Q$ 为流量， $C$ 为浓度， $r$ 为反应速率。

推流式反应器（PFR）：
$∂C∂t+u∂C∂x=r\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} = r$

实际反应器：
实际反应器介于CSTR和PFR之间，存在：

短路流
死区
回流混合
非理想流动

2.2 停留时间分布（RTD）

停留时间分布函数：

F曲线（累积分布）：
$\frac{C(t)}{C_0}$

E曲线（概率密度）：
$\frac{dF(t)}{dt}$

平均停留时间：
$tˉ=∫0∞tE(t)dt=VQ\bar{t} = \int_0^\infty t E(t) dt = \frac{V}{Q}$

方差：
$σt2=∫0∞(t−tˉ)2E(t)dt\sigma_t^2 = \int_0^\infty (t - \bar{t})^2 E(t) dt$

无量纲方差：
$σθ2=σt2tˉ2\sigma_\theta^2 = \frac{\sigma_t^2}{\bar{t}^2}$

$σθ2=0\sigma_\theta^2 = 0$ ：理想PFR
$σθ2=1\sigma_\theta^2 = 1$ ：理想CSTR
$\sigma_\theta^2 < 1$ ：实际反应器

2.3 流动模型

串联CSTR模型：
$\frac{n^n t^{n-1}}{(n-1)! \bar{t}^n} \exp\left(-\frac{nt}{\bar{t}}\right)$

其中， $n$ 为串联级数， $\to \infty$ 时趋近于PFR。

分散模型：
$∂C∂t+u∂C∂x=DL∂2C∂x2\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} = D_L \frac{\partial^2 C}{\partial x^2}$

其中， $D_L$ 为纵向分散系数。

Peclet数：
$\frac{uL}{D_L}$

$\to \infty$ ：理想PFR
$\to 0$ ：理想CSTR

三、活性污泥法仿真

3.1 活性污泥工艺原理

活性污泥法是利用微生物代谢作用去除水中有机污染物和氮磷的生物处理技术。

主要微生物：

异养菌：降解有机物
自养硝化菌：氨氮氧化
反硝化菌：硝酸盐还原
聚磷菌：过量吸磷

工艺参数：

污泥龄（SRT）：微生物在系统中的平均停留时间
水力停留时间（HRT）：污水在反应器中的停留时间
污泥回流比：回流污泥与进水流量之比
溶解氧（DO）：好氧区溶解氧浓度

3.2 ASM模型族

ASM1模型（活性污泥1号模型）：

包括13个组分和8个过程：

组分：

$S_I$ ：可溶性惰性有机物
$S_S$ ：易生物降解底物
$X_I$ ：颗粒性惰性有机物
$X_S$ ：慢速可生物降解底物
$X_{BH}$ ：异养菌
$X_{BA}$ ：自养菌
$X_P$ ：微生物残骸
$S_O$ ：溶解氧
$S_{NO}$ ：硝酸盐和亚硝酸盐氮
$S_{NH}$ ：氨氮
$S_{ND}$ ：可溶性可生物降解有机氮
$X_{ND}$ ：颗粒性可生物降解有机氮
$S_{ALK}$ ：碱度

典型过程速率方程：

异养菌好氧生长：
$ρ1=μHSSKS+SSSOKOH+SOXBH\rho_1 = \mu_H \frac{S_S}{K_S + S_S} \frac{S_O}{K_{OH} + S_O} X_{BH}$

异养菌缺氧生长：
$ρ2=μHSSKS+SSKOHKOH+SOSNOKNO+SNOηgXBH\rho_2 = \mu_H \frac{S_S}{K_S + S_S} \frac{K_{OH}}{K_{OH} + S_O} \frac{S_{NO}}{K_{NO} + S_{NO}} \eta_g X_{BH}$

自养菌好氧生长：
$ρ3=μASNHKNH+SNHSOKOA+SOXBA\rho_3 = \mu_A \frac{S_{NH}}{K_{NH} + S_{NH}} \frac{S_O}{K_{OA} + S_O} X_{BA}$

3.3 生物膜模型

生物膜特性：

微生物固定在载体表面
传质限制显著
分层结构明显

一维生物膜模型：
$∂Ci∂t=Deff∂2Ci∂z2+ri\frac{\partial C_i}{\partial t} = D_{eff} \frac{\partial^2 C_i}{\partial z^2} + r_i$

其中， $D_{eff}$ 为有效扩散系数， $z$ 为生物膜深度。

传质边界层：
$N = k_L (C_b - C_s)$

其中， $k_L$ 为液膜传质系数， $C_b$ 为主体浓度， $C_s$ 为表面浓度。

四、膜分离过程仿真

4.1 膜分离基本原理

膜分离类型：

微滤（MF）：0.1-10 μm
超滤（UF）：1-100 nm
纳滤（NF）：0.5-2 nm
反渗透（RO）：< 0.5 nm

传质机理：

筛分效应
溶解-扩散
道南效应

4.2 膜通量模型

阻力模型：
$\frac{\Delta P}{\mu (R_m + R_c + R_f)}$

其中：

$J$ ：膜通量
$ΔP\Delta P$ ：跨膜压差
$μ\mu$ ：粘度
$R_m$ ：膜阻力
$R_c$ ：浓差极化层阻力
$R_f$ ：膜污染阻力

浓差极化：
$\ln\left(\frac{C_w - C_p}{C_b - C_p}\right)$

其中：

$k$ ：传质系数
$C_w$ ：膜表面浓度
$C_b$ ：主体浓度
$C_p$ ：透过液浓度

4.3 膜污染模型

标准堵塞模型：
$d2tdV2=kbdtdV\frac{d^2 t}{dV^2} = k_b \frac{dt}{dV}$

中间堵塞模型：
$d2tdV2=ki(dtdV)2\frac{d^2 t}{dV^2} = k_i \left(\frac{dt}{dV}\right)^2$

滤饼过滤模型：
$dtdV=μRmAΔP+μαCb2A2ΔPV\frac{dt}{dV} = \frac{\mu R_m}{A \Delta P} + \frac{\mu \alpha C_b}{2 A^2 \Delta P} V$

完全堵塞模型：
$dtdV=μRmAΔPexp⁡(ksV)\frac{dt}{dV} = \frac{\mu R_m}{A \Delta P} \exp(k_s V)$

五、污染物迁移转化

5.1 对流-扩散-反应方程

基本方程：
$∂C∂t+∇⋅(uC)=∇⋅(D∇C)+r\frac{\partial C}{\partial t} + \nabla \cdot (\mathbf{u} C) = \nabla \cdot (D \nabla C) + r$

一维形式：
$∂C∂t+u∂C∂x=D∂2C∂x2+r\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} = D \frac{\partial^2 C}{\partial x^2} + r$

各项物理意义：

累积项：浓度随时间变化
对流项：随流体运动
扩散项：分子扩散和湍流扩散
反应项：生化反应

5.2 吸附过程

Langmuir等温线：
$\frac{q_{max} K_L C}{1 + K_L C}$

Freundlich等温线：
$q = K_F C^{1/n}$

吸附动力学：

准一级动力学：
$dqdt=k1(qe−q)\frac{dq}{dt} = k_1 (q_e - q)$

准二级动力学：
$dqdt=k2(qe−q)2\frac{dq}{dt} = k_2 (q_e - q)^2$

5.3 降解动力学

一级降解：
$dCdt=−kC\frac{dC}{dt} = -kC$
$C = C_0 \exp(-kt)$

二级降解：
$dCdt=−kC2\frac{dC}{dt} = -kC^2$
$1C=1C0+kt\frac{1}{C} = \frac{1}{C_0} + kt$

Monod动力学（微生物降解）：
$μ=μmaxSKS+S\mu = \mu_{max} \frac{S}{K_S + S}$

六、人工湿地仿真

6.1 人工湿地类型

表面流人工湿地：

水流在地表流动
类似自然湿地
氧传递效率低

潜流人工湿地：

水流通过基质
水平潜流或垂直潜流
氧传递效率较高

6.2 湿地水力学

Darcy定律：
$\frac{dh}{dl}$

其中：

$Q$ ：流量
$K$ ：水力传导系数
$A$ ：横截面积
$d h / d l$ ：水力梯度

停留时间分布：
人工湿地存在明显的非理想流动，常用 tanks-in-series 模型描述。

6.3 污染物去除机理

物理过程：

沉淀
过滤
吸附

化学过程：

沉淀反应
氧化还原
络合反应

生物过程：

微生物降解
植物吸收
硝化反硝化

一级去除模型：
$CoutCin=exp⁡(−kq)\frac{C_{out}}{C_{in}} = \exp\left(-\frac{k}{q}\right)$

其中， $k$ 为去除速率常数， $q$ 为水力负荷。

k-C*模型：
$dCdt=−k(C−C∗)\frac{dC}{dt} = -k(C - C^*)$

其中， $C^*$ 为背景浓度。

七、地下水污染修复

7.1 地下水流动

Darcy定律：
$q=−K∇h\mathbf{q} = -K \nabla h$

其中， $q\mathbf{q}$ 为达西流速， $K$ 为水力传导系数， $h$ 为水头。

连续性方程：
$Ss∂h∂t=∇⋅(K∇h)+WS_s \frac{\partial h}{\partial t} = \nabla \cdot (K \nabla h) + W$

其中， $S_s$ 为比储水系数， $W$ 为源汇项。

7.2 溶质运移

对流-弥散方程（ADE）：
$∂C∂t=∇⋅(D∇C)−∇⋅(vC)+r\frac{\partial C}{\partial t} = \nabla \cdot (D \nabla C) - \nabla \cdot (\mathbf{v} C) + r$

水动力弥散：
$D_m + \alpha_L v$

其中：

$D_m$ ：分子扩散系数
$αL\alpha_L$ ：纵向弥散度
$v$ ：孔隙流速

7.3 多相流

非水相液体（NAPL）：

轻非水相液体（LNAPL）：密度小于水
重非水相液体（DNAPL）：密度大于水

多相流方程：
$∂(ϕραSα)∂t=−∇⋅(ραqα)+rα\frac{\partial (\phi \rho_\alpha S_\alpha)}{\partial t} = -\nabla \cdot (\rho_\alpha \mathbf{q}_\alpha) + r_\alpha$

其中， $α\alpha$ 为相（水、NAPL、气）， $ϕ\phi$ 为孔隙度， $S$ 为饱和度。

八、数值方法

8.1 有限差分法

对流项离散：

迎风格式：
$\frac{\partial C}{\partial x} \approx u \frac{C_i - C_{i-1}}{\Delta x} \quad (u > 0)$

中心差分：
$\frac{\partial C}{\partial x} \approx u \frac{C_{i+1} - C_{i-1}}{2\Delta x}$

QUICK格式：
$\frac{\partial C}{\partial x} \approx u \frac{3C_{i+1} + 3C_i - 7C_{i-1} + C_{i-2}}{8\Delta x}$

扩散项离散：
$\frac{\partial^2 C}{\partial x^2} \approx D \frac{C_{i+1} - 2C_i + C_{i-1}}{\Delta x^2}$

8.2 有限体积法

通量计算：
$∂∂t∫VCdV=−∮SJ⋅dS+∫VrdV\frac{\partial}{\partial t} \int_V C dV = -\oint_S \mathbf{J} \cdot d\mathbf{S} + \int_V r dV$

界面通量：
$Ji+1/2=ui+1/2Ci+1/2−Di+1/2Ci+1−CiΔxJ_{i+1/2} = u_{i+1/2} C_{i+1/2} - D_{i+1/2} \frac{C_{i+1} - C_i}{\Delta x}$

8.3 特征线法

特征方程：
$dxdt=u\frac{dx}{dt} = u$
$dCdt=D∂2C∂x2+r\frac{dC}{dt} = D \frac{\partial^2 C}{\partial x^2} + r$

随机行走法：
$\Delta t) = x(t) + u \Delta t + Z \sqrt{2D\Delta t}$

其中， $Z$ 为标准正态分布随机数。

九、工程案例分析

9.1 污水处理厂优化设计

问题描述：
某城市污水处理厂设计处理能力10万吨/日，进水COD 400 mg/L，氨氮35 mg/L，总磷5 mg/L。要求出水达到一级A标准。

仿真方法：

建立ASM2d模型
模拟不同工艺配置（A²/O、氧化沟、SBR）
优化运行参数（DO、SRT、内回流比）

优化结果：

推荐工艺：A²/O
最佳SRT：15天
好氧区DO：2 mg/L
内回流比：300%
出水COD < 50 mg/L，氨氮 < 5 mg/L，总磷 < 0.5 mg/L

9.2 人工湿地设计优化

问题描述：
设计处理农村生活污水的人工湿地，流量100 m³/d，进水COD 200 mg/L，总氮40 mg/L。

仿真方法：

建立湿地水力学模型
模拟污染物去除过程
优化湿地尺寸和填料

设计结果：

类型：水平潜流人工湿地
面积：500 m²
水力负荷：0.2 m/d
填料：砾石+沸石
出水COD < 50 mg/L，总氮 < 15 mg/L

9.3 地下水污染修复

问题描述：
某工业区地下水受到氯代烃污染，需要设计修复方案。

仿真方法：

建立地下水流动模型
模拟DNAPL运移
评估不同修复技术（抽出-处理、原位生物修复、PRB）

修复方案：

抽出-处理系统：5口抽水井
原位生物修复：注入电子供体
PRB：渗透反应墙
预计修复时间：5年

十、前沿技术与发展趋势

10.1 智慧水务

数字孪生：

实时监测数据驱动
虚拟-现实交互
预测性维护

人工智能应用：

水质预测模型
异常检测
智能加药控制
能耗优化

10.2 新型处理技术

高级氧化：

臭氧氧化
Fenton反应
光催化
电化学氧化

膜生物反应器（MBR）：

高效固液分离
占地面积小
出水水质好

厌氧氨氧化：

节能脱氮
污泥产量低
温室气体减排

10.3 资源回收

磷回收：

鸟粪石结晶
磷酸钙沉淀
生物聚磷

能源回收：

厌氧消化产甲烷
微生物燃料电池
热水解

水资源回用：

再生水利用
零液体排放
海水淡化

10.4 气候变化适应

暴雨管理：

低影响开发（LID）
海绵城市
雨水调蓄

极端天气应对：

洪涝风险评估
应急处理
韧性设计

十一、总结

水处理与环境工程仿真是解决水环境问题的重要技术手段。通过建立准确的数学模型和采用高效的数值方法，可以深入理解水处理过程的物理化学机制，优化工艺设计，提高处理效率，降低运行成本。

随着计算能力的提升和人工智能技术的发展，水处理仿真正朝着多尺度、多物理场、实时化和智能化方向发展。数字孪生技术将虚拟仿真与物理系统深度融合，为智慧水务建设提供了强有力的技术支撑。

未来，水处理仿真将在新型污染物控制、资源能源回收、气候变化适应等方面发挥更加重要的作用，推动水处理技术向高效、低碳、循环的方向发展，为实现水资源的可持续利用做出贡献。

6 深入理论分析

6.1 数学模型推导

流体力学问题的数学描述基于守恒定律和本构关系。控制方程的一般形式可以表示为：

$∂(ρϕ)∂t+∇⋅(ρuϕ)=∇⋅(Γ∇ϕ)+Sϕ\frac{\partial (\rho \phi)}{\partial t} + \nabla \cdot (\rho \mathbf{u} \phi) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi$

其中， $ϕ\phi$ 表示通用变量， $ρ\rho$ 是密度， $u\mathbf{u}$ 是速度矢量， $Γ\Gamma$ 是扩散系数， $SϕS_\phi$ 是源项。

对于稳态问题，控制方程简化为：

$∇⋅(ρuϕ)=∇⋅(Γ∇ϕ)+Sϕ\nabla \cdot (\rho \mathbf{u} \phi) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi$

边界条件的数学表达：

Dirichlet边界条件：
$ϕ=ϕspecified在ΓD上\phi = \phi_{specified} \quad \text{在} \quad \Gamma_D \text{上}$

Neumann边界条件：
$−Γ∂ϕ∂n=qspecified在ΓN上-\Gamma \frac{\partial \phi}{\partial n} = q_{specified} \quad \text{在} \quad \Gamma_N \text{上}$

Robin边界条件：
$−Γ∂ϕ∂n=h(ϕ−ϕ∞)在ΓR上-\Gamma \frac{\partial \phi}{\partial n} = h(\phi - \phi_\infty) \quad \text{在} \quad \Gamma_R \text{上}$

6.2 数值离散方法

有限差分法的离散格式：

时间离散：

显式格式： $ϕn+1=ϕn+Δt⋅f(ϕn)\phi^{n+1} = \phi^n + \Delta t \cdot f(\phi^n)$
隐式格式： $ϕn+1=ϕn+Δt⋅f(ϕn+1)\phi^{n+1} = \phi^n + \Delta t \cdot f(\phi^{n+1})$
Crank-Nicolson格式： $ϕn+1=ϕn+Δt2[f(ϕn)+f(ϕn+1)]\phi^{n+1} = \phi^n + \frac{\Delta t}2[f(\phi^n) + f(\phi^{n+1})]$

空间离散：

中心差分： $∂ϕ∂x≈ϕi+1−ϕi−12Δx\frac{\partial \phi}{\partial x} \approx \frac{\phi_{i+1} - \phi_{i-1}}{2\Delta x}$
迎风格式： $∂ϕ∂x≈ϕi−ϕi−1Δx\frac{\partial \phi}{\partial x} \approx \frac{\phi_i - \phi_{i-1}}{\Delta x}$ （ $u > 0$ 时）

6.3 稳定性与收敛性分析

Von Neumann稳定性分析：

对于一维热传导方程的显式格式，稳定性条件为：

$\frac{\alpha \Delta t}{(\Delta x)^2} \leq \frac12$

其中 $F o$ 是傅里叶数， $α\alpha$ 是热扩散系数。

收敛性条件：

根据Lax等价定理，对于适定的线性初值问题，一致性和稳定性是收敛性的充分必要条件。

误差分析：

截断误差： $τ=O(Δt)+O((Δx)2)\tau = O(\Delta t) + O((\Delta x)^2)$

累积误差： $ϵ=O(Δt)+O((Δx)2)\epsilon = O(\Delta t) + O((\Delta x)^2)$

7 工程案例分析

7.1 案例背景

本案例研究流体力学在管道流动优化中的应用。该问题涉及复杂的物理过程和边界条件，具有典型的工程实际意义。

问题描述：

几何尺寸：根据实际工程确定
材料属性：标准工程材料
边界条件：典型工况条件
初始条件：稳态或瞬态初始场

物理过程分析：

该问题涉及以下物理过程：

场变量的空间分布与演化
边界上的能量/质量交换
多物理场之间的耦合作用

这些过程之间存在强烈的耦合关系，需要采用多场耦合方法求解。

7.2 数值建模

网格划分：

采用结构化或非结构化网格，根据问题复杂度确定网格数量。在关键区域进行局部加密，确保计算精度。

网格质量指标：

最小单元尺寸：根据梯度确定
最大单元尺寸：保证整体精度
长宽比：小于5-10
扭曲度：小于0.5-0.6

求解器设置：

时间步长：根据稳定性条件确定
收敛准则：残差下降3-6个数量级
迭代次数：根据问题复杂度
计算时间：取决于网格规模和硬件

7.3 结果分析

稳态结果：

场变量分布云图显示了物理量的空间分布特征。最大值出现在关键区域，最小值出现在边界区域。

瞬态演化：

随时间演化，系统逐渐趋于稳态。在早期阶段，场变量变化剧烈；在后期阶段，变化趋于平缓。

参数敏感性分析：

研究了关键参数对结果的影响：

当参数增加20%时，输出响应变化
当参数增加50%时，输出响应显著变化

验证与确认：

将数值结果与理论解析解或实验数据对比，验证模型的准确性。相对误差通常在1-5%范围内。

7.4 优化设计

基于仿真结果，提出了以下优化建议：

几何优化：改进几何形状，减小不利因素
材料选择：采用性能更优的材料
工艺参数：优化操作条件，提高效率

优化后的设计方案相比原方案，性能显著提升，满足设计要求。

8 高级主题与前沿研究

8.1 多尺度建模方法

多尺度问题涉及从微观到宏观的多个尺度，需要采用特殊的方法处理：

尺度分离方法：

均质化方法（Homogenization）
代表性体积单元（RVE）
多尺度有限元法

耦合策略：

顺序多尺度：微观→介观→宏观
并发多尺度：同时求解多个尺度
自适应多尺度：根据局部特征动态选择尺度

8.2 不确定性量化

工程问题中存在各种不确定性来源：

不确定性类型：

偶然不确定性（Aleatory）：固有的随机性
认知不确定性（Epistemic）：知识不足导致的

分析方法：

蒙特卡洛模拟（MCS）
多项式混沌展开（PCE）
随机有限元法（SFEM）
贝叶斯推断

灵敏度分析：

Sobol指数用于量化各输入参数对输出的影响：

$Si=VXi(EX∼i[Y∣Xi])V(Y)S_i = \frac{V_{X_i}(E_{X_{\sim i}}[Y|X_i])}{V(Y)}$

8.3 高性能计算

并行计算策略：

域分解：将计算域划分为子域，分配给不同处理器
任务并行：不同任务并行执行
数据并行：同一任务在不同数据上并行

加速技术：

GPU加速：利用CUDA/OpenCL实现大规模并行
向量化：利用SIMD指令加速
混合精度：平衡精度和速度

可扩展性分析：

加速比： $\frac{T_1}{T_p}$

效率： $\frac{S(p)}{p}$

8.4 机器学习融合

数据驱动建模：

神经网络代理模型
高斯过程回归
支持向量机

物理信息神经网络（PINN）：

将物理约束嵌入损失函数：

$L=Ldata+λLPDE\mathcal{L} = \mathcal{L}_{data} + \lambda \mathcal{L}_{PDE}$

强化学习优化：

智能体通过与环境交互学习最优策略，用于：

拓扑优化
参数优化
控制优化

8.5 发展趋势

本领域的未来发展方向：

数字孪生：虚实融合，实时仿真
云仿真：基于云计算的仿真服务
边缘计算：嵌入式实时仿真
量子计算：量子算法加速
自主仿真：AI驱动的自动化仿真

附录

附录A：常用物理常数

常数	符号	数值	单位
重力加速度	$g$	9.80665	m/s²
标准大气压	$p_{atm}$	101325	Pa
水的密度（4°C）	$ρwater\rho_{water}$	1000	kg/m³
空气密度（标准状态）	$ρair\rho_{air}$	1.225	kg/m³
水的比热容	$c_{p,water}$	4186	J/(kg·K)
空气比热容	$c_{p,air}$	1005	J/(kg·K)
水的导热系数	$k_{water}$	0.6	W/(m·K)
空气导热系数	$k_{air}$	0.026	W/(m·K)
斯特藩-玻尔兹曼常数	$σ\sigma$	5.67×10⁻⁸	W/(m²·K⁴)

附录B：单位换算

长度：

1 m = 100 cm = 1000 mm
1 inch = 25.4 mm
1 ft = 0.3048 m

压力：

1 Pa = 1 N/m²
1 bar = 10⁵ Pa
1 atm = 101325 Pa
1 psi = 6894.76 Pa

能量：

1 J = 1 N·m
1 cal = 4.184 J
1 kWh = 3.6×10⁶ J

功率：

1 W = 1 J/s
1 hp = 745.7 W

附录C：数值方法参数选择指南

时间步长选择：

显式格式： $Δt≤(Δx)22α\Delta t \leq \frac{(\Delta x)^2}{2\alpha}$ （扩散问题）
显式格式： $Δt≤Δx∣u∣\Delta t \leq \frac{\Delta x}{|u|}$ （对流问题，CFL条件）
隐式格式：稳定性较好，可适当增大

网格尺寸选择：

根据梯度大小局部加密
边界层区域需要细化
网格无关性验证

收敛准则：

残差下降3-6个数量级
监测关键物理量的变化
质量/能量守恒检查

附录D：编程技巧与最佳实践

代码结构：

# 1. 参数设置
# 2. 网格生成
# 3. 初始化
# 4. 时间推进
# 5. 后处理
# 6. 结果输出

性能优化：

使用NumPy数组运算替代循环
预分配数组内存
向量化计算
使用JIT编译（Numba）

调试技巧：

从小规模问题开始
与解析解对比
网格收敛性验证
物理合理性检查

附录E：扩展阅读资源

经典教材：

《数值传热学》— 陶文铨
《有限元法原理与应用》— 朱伯芳
《计算流体力学基础》— 刘儒勋
《多物理场耦合分析》— 相关专著

在线资源：

FEniCS Project（开源有限元软件）
OpenFOAM（开源CFD软件）
deal.II（自适应有限元库）

学术期刊：

Journal of Computational Physics
Computer Methods in Applied Mechanics and Engineering
International Journal for Numerical Methods in Engineering
Applied Mathematical Modelling

9 扩展理论深化

9.1 守恒方程的推导

从基本物理定律出发，可以推导出控制方程。

质量守恒方程：

考虑一个控制体，质量守恒要求：

$∂ρ∂t+∇⋅(ρu)=0\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$

对于不可压缩流体，密度 $ρ\rho$ 为常数，方程简化为：

$∇⋅u=0\nabla \cdot \mathbf{u} = 0$

动量守恒方程（Navier-Stokes方程）：

$ρ(∂u∂t+u⋅∇u)=−∇p+∇⋅τ+ρg\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g}$

其中， $p$ 是压力， $τ\boldsymbol{\tau}$ 是粘性应力张量， $g\mathbf{g}$ 是重力加速度。

对于牛顿流体，粘性应力与应变率成正比：

$τ=μ(∇u+(∇u)T)−23μ(∇⋅u)I\boldsymbol{\tau} = \mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right) - \frac{2}{3}\mu (\nabla \cdot \mathbf{u})\mathbf{I}$

能量守恒方程：

$ρcp(∂T∂t+u⋅∇T)=∇⋅(k∇T)+Φ+q˙\rho c_p \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) = \nabla \cdot (k \nabla T) + \Phi + \dot{q}$

其中， $c_p$ 是比热容， $k$ 是导热系数， $Φ\Phi$ 是粘性耗散函数， $q˙\dot{q}$ 是体积热源。

9.2 无量纲化与相似准则

通过无量纲化可以简化问题并提取关键参数。

特征尺度：

长度尺度： $L$
速度尺度： $U$
时间尺度： $t_0 = L/U$
压力尺度： $p0=ρU2p_0 = \rho U^2$
温度尺度： $ΔT\Delta T$

无量纲变量：

$x∗=xL,u∗=uU,t∗=tt0,p∗=pp0x^* = \frac{x}{L}, \quad u^* = \frac{u}{U}, \quad t^* = \frac{t}{t_0}, \quad p^* = \frac{p}{p_0}$

重要无量纲数：

雷诺数（Reynolds Number）：
$\frac{\rho U L}{\mu} = \frac{UL}{\nu}$

表示惯性力与粘性力之比。 $\ll 1$ 时为蠕动流， $\gg 1$ 时为湍流。

普朗特数（Prandtl Number）：
$\frac{\nu}{\alpha} = \frac{c_p \mu}{k}$

表示动量扩散与热扩散之比。对于空气， $\approx 0.7$ ；对于水， $\approx 7$ 。

格拉晓夫数（Grashof Number）：
$\frac{g \beta \Delta T L^3}{\nu^2}$

表示浮升力与粘性力之比，用于自然对流。

瑞利数（Rayleigh Number）：
$\cdot Pr = \frac{g \beta \Delta T L^3}{\nu \alpha}$

努塞尔数（Nusselt Number）：
$\frac{h L}{k}$

表示对流换热与纯导热之比。

毕渥数（Biot Number）：
$\frac{h L_c}{k_s}$

表示表面换热热阻与内部导热热阻之比。

傅里叶数（Fourier Number）：
$\frac{\alpha t}{L^2}$

表示热传导速率与储热速率之比。

9.3 边界层理论

速度边界层：

在固体壁面附近，流体速度从壁面的无滑移条件（ $u = 0$ ）逐渐过渡到主流速度。边界层厚度 $δ\delta$ 定义为速度达到主流速度99%的位置。

对于层流平板边界层，Blasius解给出：
$δ=5.0xRex\delta = \frac{5.0 x}{\sqrt{Re_x}}$

热边界层：

类似地，温度从壁面温度逐渐过渡到主流温度。热边界层厚度 $δT\delta_T$ 与速度边界层厚度之比取决于普朗特数：

$δTδ≈Pr−1/3\frac{\delta_T}{\delta} \approx Pr^{-1/3}$

边界层方程：

通过量级分析，可以简化Navier-Stokes方程得到边界层方程：

$\frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = -\frac{1}{\rho} \frac{dp}{dx} + \nu \frac{\partial^2 u}{\partial y^2}$

$∂u∂x+∂v∂y=0\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0$

9.4 湍流模型

当 $Re_{crit} \approx 2300$ （管道流动）时，流动转变为湍流。

雷诺平均Navier-Stokes（RANS）方程：

将瞬时量分解为平均量和脉动量：
$\bar{u} + u'$

对NS方程进行雷诺平均：

$ρ(∂uˉi∂t+uˉj∂uˉi∂xj)=−∂pˉ∂xi+∂∂xj(μ∂uˉi∂xj−ρui′uj′‾)\rho \left( \frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j \frac{\partial \bar{u}_i}{\partial x_j} \right) = -\frac{\partial \bar{p}}{\partial x_i} + \frac{\partial}{\partial x_j} \left( \mu \frac{\partial \bar{u}_i}{\partial x_j} - \rho \overline{u'_i u'_j} \right)$

其中，$ - \rho \overline{u’_i u’_j}$是雷诺应力，需要模型封闭。

k-\varepsilon湍流模型：

引入湍动能 $k$ 和湍流耗散率 $ε\varepsilon$ ：

$∂(ρk)∂t+∂(ρuˉjk)∂xj=∂∂xj[(μ+μtσk)∂k∂xj]+Pk−ρε\frac{\partial (\rho k)}{\partial t} + \frac{\partial (\rho \bar{u}_j k)}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right] + P_k - \rho \varepsilon$

$∂(ρε)∂t+∂(ρuˉjε)∂xj=∂∂xj[(μ+μtσε)∂ε∂xj]+Cε1εkPk−Cε2ρε2k\frac{\partial (\rho \varepsilon)}{\partial t} + \frac{\partial (\rho \bar{u}_j \varepsilon)}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \mu + \frac{\mu_t}{\sigma_\varepsilon} \right) \frac{\partial \varepsilon}{\partial x_j} \right] + C_{\varepsilon 1} \frac{\varepsilon}{k} P_k - C_{\varepsilon 2} \rho \frac{\varepsilon^2}{k}$

其中， $μt=Cμρk2ε\mu_t = C_\mu \rho \frac{k^2}{\varepsilon}$ 是湍流粘度。

标准模型常数： $Cμ=0.09C_\mu = 0.09$ ， $Cε1=1.44C_{\varepsilon 1} = 1.44$ ， $Cε2=1.92C_{\varepsilon 2} = 1.92$ ， $σk=1.0\sigma_k = 1.0$ ， $σε=1.3\sigma_\varepsilon = 1.3$ 。

大涡模拟（LES）：

直接求解大尺度涡，对小尺度涡进行建模：

$∂uˉi∂t+∂(uˉiuˉj)∂xj=−1ρ∂pˉ∂xi+ν∂2uˉi∂xj∂xj−∂τijsgs∂xj\frac{\partial \bar{u}_i}{\partial t} + \frac{\partial (\bar{u}_i \bar{u}_j)}{\partial x_j} = -\frac{1}{\rho} \frac{\partial \bar{p}}{\partial x_i} + \nu \frac{\partial^2 \bar{u}_i}{\partial x_j \partial x_j} - \frac{\partial \tau_{ij}^{sgs}}{\partial x_j}$

其中， $τijsgs\tau_{ij}^{sgs}$ 是亚格子应力。

9.5 数值方法进阶

高阶格式：

QUICK格式： $ϕe=38ϕP+34ϕE−18ϕW\phi_e = \frac{3}{8}\phi_P + \frac{3}{4}\phi_E - \frac{1}{8}\phi_W$
MUSCL格式：单调迎风格式
TVD格式：总变差减小格式
ENO/WENO格式：本质无振荡/加权本质无振荡格式

多重网格方法：

利用不同尺度的网格加速收敛：

在细网格上迭代，得到近似解
将残差限制到粗网格
在粗网格上求解误差方程
将误差延拓到细网格，修正解
重复直到收敛

收敛速度： $O (N)$ ，其中 $N$ 是网格节点数。

预处理技术：

对于大型稀疏线性方程组 $A x = b$ ，预处理可以改善系数矩阵的条件数：

$M^{-1}Ax = M^{-1}b$

常用预处理器：

Jacobi预处理： $\text{diag}(A)$
ILU预处理：不完全LU分解
多重网格预处理

9.6 误差估计与自适应

后验误差估计：

基于计算结果估计误差：

$η=∥∇ϕh−gh∥\eta = \left\| \nabla \phi_h - \mathbf{g}_h \right\|$

其中， $gh\mathbf{g}_h$ 是恢复梯度（如ZZ恢复）。

自适应网格细化：

根据误差指示子局部加密网格：

计算误差指示子
标记误差大的单元
细化标记单元
重新求解
重复直到满足精度要求

h-细化：减小单元尺寸
p-细化：提高多项式阶数
r-细化：重新分布节点

9.7 验证与确认

代码验证（Verification）：

验证数值解是否正确实现了数学模型：

网格收敛性研究
与解析解对比
与基准解对比
守恒律检查

模型确认（Validation）：

验证数学模型是否准确描述物理现象：

与实验数据对比
与现场数据对比
不确定性量化

理查森外推：

利用不同网格上的解估计精确解：

$fexact≈f1+f1−f2rp−1f_{exact} \approx f_1 + \frac{f_1 - f_2}{r^p - 1}$

其中， $r$ 是网格细化比， $p$ 是收敛阶数。

9.8 工程应用实例

换热器设计：

管壳式换热器的传热计算：

$\Delta T_{lm}$

其中， $ΔTlm=ΔT1−ΔT2ln⁡(ΔT1/ΔT2)\Delta T_{lm} = \frac{\Delta T_1 - \Delta T_2}{\ln(\Delta T_1/\Delta T_2)}$ 是对数平均温差。

总传热系数：
$1UA=1hiAi+ln⁡(do/di)2πkL+1hoAo\frac{1}{UA} = \frac{1}{h_i A_i} + \frac{\ln(d_o/d_i)}{2\pi k L} + \frac{1}{h_o A_o}$

电子冷却：

芯片结温计算：

$Tj=Ta+P⋅θjaT_j = T_a + P \cdot \theta_{ja}$

其中， $θja\theta_{ja}$ 是结到环境的热阻。

散热器的优化设计需要考虑：

鳍片几何形状
鳍片间距
材料选择
冷却方式（自然对流/强制对流/液冷）

建筑能耗：

围护结构传热：

$Q = U A (T_{in} - T_{out})$

整体传热系数：
$\frac{1}{R_{si} + \sum R_i + R_{se}}$

其中， $R_{si}$ 和 $R_{se}$ 是内外表面热阻， $R_i$ 是各层材料热阻。

"""
主题054：水处理与环境工程仿真
Water Treatment and Environmental Engineering Simulation

本程序演示水处理过程的仿真方法，包含以下4个案例：
1. 活性污泥反应器仿真（ASM简化模型）
2. 膜分离过程与膜污染分析
3. 人工湿地污染物去除仿真
4. 地下水污染物迁移扩散

作者：仿真领域专家
日期：2026-03-03
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
from scipy.interpolate import interp1d
import matplotlib.patches as patches
from matplotlib.patches import FancyBboxPatch, Circle, Rectangle
from matplotlib.collections import LineCollection
import math
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 输出目录
OUTPUT_DIR = 'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题054\\output\\'

# ==================== 案例1: 活性污泥反应器仿真 ====================
def case1_activated_sludge():
    """案例1: 活性污泥反应器仿真（ASM简化模型）"""
    print("\n" + "="*60)
    print("案例1: 活性污泥反应器仿真")
    print("="*60)
    
    # 反应器参数
    reactor_params = {
        'V': 10000,      # 反应器体积 (m³)
        'Q': 2000,       # 进水流量 (m³/d)
        'X_recycle': 0.5, # 污泥回流比
        'DO_set': 2.0,   # 溶解氧设定值 (mg/L)
    }
    
    # 进水水质
    influent = {
        'COD': 400,      # 化学需氧量 (mg/L)
        'NH4': 35,       # 氨氮 (mg/L)
        'NO3': 0,        # 硝酸盐 (mg/L)
        'PO4': 5,        # 磷酸盐 (mg/L)
        'TSS': 200,      # 总悬浮固体 (mg/L)
    }
    
    # ASM1简化模型参数
    asm_params = {
        'mu_H': 4.0,     # 异养菌最大比生长速率 (1/d)
        'K_S': 10,       # 半饱和常数 (mg COD/L)
        'K_OH': 0.2,     # 氧半饱和常数 (mg O2/L)
        'Y_H': 0.67,     # 异养菌产率系数
        'b_H': 0.3,      # 异养菌衰减系数 (1/d)
        'mu_A': 0.5,     # 自养菌最大比生长速率 (1/d)
        'K_NH': 1.0,     # 氨氮半饱和常数 (mg N/L)
        'K_OA': 0.4,     # 自养菌氧半饱和常数 (mg O2/L)
        'Y_A': 0.24,     # 自养菌产率系数
        'b_A': 0.05,     # 自养菌衰减系数 (1/d)
        'f_P': 0.08,     # 生物量转化为颗粒产物的比例
    }
    
    # 定义ODE系统
    def asm_model(y, t, reactor, influent, asm):
        """ASM简化模型"""
        S_S, X_BH, S_NH, S_NO, X_BA, S_O = y
        
        # 水力停留时间
        HRT = reactor['V'] / reactor['Q']
        
        # 异养菌好氧生长
        rho_1 = asm['mu_H'] * S_S / (asm['K_S'] + S_S) * S_O / (asm['K_OH'] + S_O) * X_BH
        
        # 异养菌缺氧生长（简化）
        eta = 0.8  # 缺氧生长校正因子
        rho_2 = asm['mu_H'] * S_S / (asm['K_S'] + S_S) * asm['K_OH'] / (asm['K_OH'] + S_O) * \
                S_NO / (0.5 + S_NO) * eta * X_BH
        
        # 自养菌好氧生长
        rho_3 = asm['mu_A'] * S_NH / (asm['K_NH'] + S_NH) * S_O / (asm['K_OA'] + S_O) * X_BA
        
        # 衰减过程
        rho_4 = asm['b_H'] * X_BH
        rho_5 = asm['b_A'] * X_BA
        
        # 微分方程
        dS_S = (influent['COD'] - S_S) / HRT - rho_1 / asm['Y_H'] - rho_2 / asm['Y_H']
        dX_BH = -X_BH / HRT + rho_1 + rho_2 - rho_4
        dS_NH = (influent['NH4'] - S_NH) / HRT - asm['Y_A'] * rho_3 + 0.08 * (rho_4 + rho_5)
        dS_NO = -S_NO / HRT + (1.0 - asm['Y_H']) / (2.86 * asm['Y_H']) * rho_2 - rho_3 / asm['Y_A']
        dX_BA = -X_BA / HRT + rho_3 - rho_5
        dS_O = (8.0 - S_O) / HRT - (1.0 - asm['Y_H']) / asm['Y_H'] * rho_1 - (4.57 - asm['Y_A']) / asm['Y_A'] * rho_3
        
        # 溶解氧控制
        dS_O += max(0, (reactor['DO_set'] - S_O) * 10)
        
        return [dS_S, dX_BH, dS_NH, dS_NO, dX_BA, dS_O]
    
    # 初始条件
    y0 = [50, 100, 10, 5, 20, 2]
    
    # 时间范围
    t = np.linspace(0, 30, 300)  # 30天
    
    # 求解ODE
    solution = odeint(asm_model, y0, t, args=(reactor_params, influent, asm_params))
    
    S_S, X_BH, S_NH, S_NO, X_BA, S_O = solution.T
    
    # 计算出水COD和总氮
    effluent_COD = S_S + 0.5 * X_BH
    total_N = S_NH + S_NO
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    fig.suptitle('案例1: 活性污泥反应器仿真 (ASM简化模型)', fontsize=14, fontweight='bold')
    
    # 1. 有机物降解
    ax1 = axes[0, 0]
    ax1.plot(t, S_S, 'b-', linewidth=2, label='易降解有机物 $S_S$')
    ax1.plot(t, effluent_COD, 'r--', linewidth=2, label='出水COD')
    ax1.axhline(y=50, color='g', linestyle=':', label='一级A标准')
    ax1.set_xlabel('时间 (d)')
    ax1.set_ylabel('浓度 (mg/L)')
    ax1.set_title('有机物降解过程')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 微生物浓度
    ax2 = axes[0, 1]
    ax2.plot(t, X_BH, 'b-', linewidth=2, label='异养菌 $X_{BH}$')
    ax2.plot(t, X_BA, 'r-', linewidth=2, label='自养菌 $X_{BA}$')
    ax2.plot(t, X_BH + X_BA, 'g--', linewidth=2, label='总生物量')
    ax2.set_xlabel('时间 (d)')
    ax2.set_ylabel('浓度 (mg/L)')
    ax2.set_title('微生物浓度变化')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 氮素转化
    ax3 = axes[0, 2]
    ax3.plot(t, S_NH, 'b-', linewidth=2, label='氨氮 $S_{NH}$')
    ax3.plot(t, S_NO, 'r-', linewidth=2, label='硝酸盐 $S_{NO}$')
    ax3.plot(t, total_N, 'g--', linewidth=2, label='总氮')
    ax3.axhline(y=5, color='orange', linestyle=':', label='氨氮标准')
    ax3.axhline(y=15, color='purple', linestyle=':', label='总氮标准')
    ax3.set_xlabel('时间 (d)')
    ax3.set_ylabel('浓度 (mg/L)')
    ax3.set_title('氮素转化过程')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 溶解氧
    ax4 = axes[1, 0]
    ax4.plot(t, S_O, 'b-', linewidth=2)
    ax4.axhline(y=reactor_params['DO_set'], color='r', linestyle='--', label='设定值')
    ax4.set_xlabel('时间 (d)')
    ax4.set_ylabel('溶解氧 (mg/L)')
    ax4.set_title('溶解氧浓度')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. 去除效率
    ax5 = axes[1, 1]
    COD_removal = (influent['COD'] - effluent_COD) / influent['COD'] * 100
    N_removal = (influent['NH4'] - total_N) / influent['NH4'] * 100
    ax5.plot(t, COD_removal, 'b-', linewidth=2, label='COD去除率')
    ax5.plot(t, N_removal, 'r-', linewidth=2, label='总氮去除率')
    ax5.set_xlabel('时间 (d)')
    ax5.set_ylabel('去除效率 (%)')
    ax5.set_title('污染物去除效率')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. 反应器示意图
    ax6 = axes[1, 2]
    ax6.set_xlim(0, 10)
    ax6.set_ylim(0, 10)
    ax6.set_aspect('equal')
    ax6.axis('off')
    ax6.set_title('活性污泥工艺示意图')
    
    # 绘制反应器
    reactor = FancyBboxPatch((2, 3), 6, 4, boxstyle="round,pad=0.1", 
                              edgecolor='black', facecolor='lightblue', linewidth=2)
    ax6.add_patch(reactor)
    ax6.text(5, 5, '曝气池\nV=10000m³', ha='center', va='center', fontsize=10)
    
    # 进水箭头
    ax6.arrow(0, 5, 1.5, 0, head_width=0.3, head_length=0.3, fc='blue', ec='blue')
    ax6.text(0.5, 5.5, '进水', fontsize=9)
    
    # 出水箭头
    ax6.arrow(8.5, 5, 1.5, 0, head_width=0.3, head_length=0.3, fc='green', ec='green')
    ax6.text(9, 5.5, '出水', fontsize=9)
    
    # 污泥回流
    ax6.annotate('', xy=(5, 3), xytext=(5, 2),
                arrowprops=dict(arrowstyle='->', color='brown', lw=2))
    ax6.text(5.5, 2.3, '污泥回流', fontsize=9)
    
    # 曝气符号
    for i in range(3):
        circle = Circle((3.5 + i*1.5, 3.5), 0.2, facecolor='white', edgecolor='black')
        ax6.add_patch(circle)
        ax6.plot([3.5 + i*1.5, 3.5 + i*1.5], [3.3, 3.1], 'k-', linewidth=1)
    ax6.text(5, 2.7, '曝气', fontsize=9, ha='center')
    
    plt.tight_layout()
    plt.savefig(f'{OUTPUT_DIR}case1_activated_sludge.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"案例1完成！图片已保存至 {OUTPUT_DIR}case1_activated_sludge.png")
    print(f"稳态出水COD: {effluent_COD[-1]:.1f} mg/L")
    print(f"稳态出水氨氮: {S_NH[-1]:.1f} mg/L")
    print(f"稳态出水总氮: {total_N[-1]:.1f} mg/L")
    print(f"COD去除率: {COD_removal[-1]:.1f}%")
    print(f"总氮去除率: {N_removal[-1]:.1f}%")

# ==================== 案例2: 膜分离过程与膜污染分析 ====================
def case2_membrane_filtration():
    """案例2: 膜分离过程与膜污染分析"""
    print("\n" + "="*60)
    print("案例2: 膜分离过程与膜污染分析")
    print("="*60)
    
    # 膜参数
    membrane = {
        'A': 0.1,        # 膜面积 (m²)
        'R_m': 1e11,     # 膜阻力 (1/m)
        'TMP': 0.5,      # 跨膜压差 (bar)
        'mu': 0.001,     # 粘度 (Pa·s)
    }
    
    # 料液参数
    feed = {
        'C_b': 10,       # 主体浓度 (g/L)
        'k': 5e-5,       # 传质系数 (m/s)
    }
    
    # 污染参数
    fouling = {
        'k_c': 1e15,     # 滤饼阻力系数
        'k_s': 0.1,      # 完全堵塞系数
        'alpha': 1e12,   # 滤饼比阻 (m/kg)
    }
    
    # 过滤时间
    t = np.linspace(0, 3600, 1000)  # 1小时
    
    # 1. 清洁膜通量（恒定压力）
    J_clean = membrane['TMP'] * 1e5 / (membrane['mu'] * membrane['R_m']) * 1000 * 3600  # L/(m²·h)
    
    # 2. 浓差极化模型
    def concentration_polarization(t, J0, C_b, k):
        """浓差极化导致通量下降"""
        # 简化模型：通量随时间指数衰减
        J = J0 * np.exp(-t / 1000)
        C_w = C_b * np.exp(J / (k * 3600 * 1000))  # 膜表面浓度
        return J, C_w
    
    J_cp, C_w = concentration_polarization(t, J_clean, feed['C_b'], feed['k'])
    
    # 3. 滤饼过滤模型
    def cake_filtration(t, J0, C_b, alpha, A, mu, TMP):
        """滤饼过滤模型"""
        V = np.zeros_like(t)
        J = np.zeros_like(t)
        J[0] = J0
        
        for i in range(1, len(t)):
            dt = t[i] - t[i-1]
            # 累积滤液体积 (L/m² -> m³/m²)
            V[i] = V[i-1] + J[i-1] * dt / 3600 / 1000
            # 滤饼阻力
            R_c = alpha * C_b * V[i] * 1000  # 转换为kg/m³
            # 通量
            J[i] = TMP * 1e5 / (mu * (membrane['R_m'] + R_c)) * 1000 * 3600
        
        return V, J
    
    V_cake, J_cake = cake_filtration(t, J_clean, feed['C_b'], fouling['alpha'], 
                                      membrane['A'], membrane['mu'], membrane['TMP'])
    
    # 4. 完全堵塞模型
    def complete_blocking(t, J0, k_s):
        """完全堵塞模型"""
        J = J0 * np.exp(-k_s * J0 * t / 3600 / 1000)
        return J
    
    J_blocking = complete_blocking(t, J_clean, fouling['k_s'])
    
    # 5. 中间堵塞模型
    def intermediate_blocking(t, J0, k_i):
        """中间堵塞模型"""
        J = J0 / (1 + k_i * J0 * t / 3600 / 1000)
        return J
    
    J_intermediate = intermediate_blocking(t, J_clean, 0.05)
    
    # 6. 标准堵塞模型
    def standard_blocking(t, J0, k_b):
        """标准堵塞模型"""
        J = J0 / (1 + k_b * J0 * t / 3600 / 1000)**2
        return J
    
    J_standard = standard_blocking(t, J_clean, 0.02)
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    fig.suptitle('案例2: 膜分离过程与膜污染分析', fontsize=14, fontweight='bold')
    
    # 1. 通量衰减曲线
    ax1 = axes[0, 0]
    ax1.plot(t/60, J_clean * np.ones_like(t), 'k--', linewidth=2, label='清洁膜')
    ax1.plot(t/60, J_cp, 'b-', linewidth=2, label='浓差极化')
    ax1.plot(t/60, J_cake, 'r-', linewidth=2, label='滤饼过滤')
    ax1.plot(t/60, J_blocking, 'g-', linewidth=2, label='完全堵塞')
    ax1.set_xlabel('时间 (min)')
    ax1.set_ylabel('膜通量 (L/(m²·h))')
    ax1.set_title('膜通量衰减曲线')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 不同污染模型比较
    ax2 = axes[0, 1]
    ax2.plot(t/60, J_cake/J_clean, 'r-', linewidth=2, label='滤饼过滤')
    ax2.plot(t/60, J_blocking/J_clean, 'g-', linewidth=2, label='完全堵塞')
    ax2.plot(t/60, J_intermediate/J_clean, 'b-', linewidth=2, label='中间堵塞')
    ax2.plot(t/60, J_standard/J_clean, 'm-', linewidth=2, label='标准堵塞')
    ax2.set_xlabel('时间 (min)')
    ax2.set_ylabel('归一化通量 J/J₀')
    ax2.set_title('不同污染模型比较')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 累积滤液体积
    ax3 = axes[0, 2]
    ax3.plot(t/60, V_cake * membrane['A'] * 1000, 'b-', linewidth=2)
    ax3.set_xlabel('时间 (min)')
    ax3.set_ylabel('累积滤液体积 (L)')
    ax3.set_title('滤饼过滤累积滤液量')
    ax3.grid(True, alpha=0.3)
    
    # 4. 膜表面浓度分布
    ax4 = axes[1, 0]
    y = np.linspace(0, 0.1, 100)  # 边界层厚度 (mm)
    # 稳态浓度分布
    C_profile = feed['C_b'] * np.exp((J_cake[-1] / 3600 / 1000) * (y/1000) / feed['k'])
    ax4.plot(C_profile, y*1000, 'b-', linewidth=2)
    ax4.axvline(x=feed['C_b'], color='r', linestyle='--', label='主体浓度')
    ax4.set_xlabel('浓度 (g/L)')
    ax4.set_ylabel('距膜表面距离 (μm)')
    ax4.set_title('浓差极化层浓度分布')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. 阻力分析
    ax5 = axes[1, 1]
    R_cake = fouling['alpha'] * feed['C_b'] * V_cake * 1000
    R_total = membrane['R_m'] + R_cake
    ax5.plot(t/60, membrane['R_m'] * np.ones_like(t)/1e11, 'b-', linewidth=2, label='膜阻力')
    ax5.plot(t/60, R_cake/1e11, 'r-', linewidth=2, label='滤饼阻力')
    ax5.plot(t/60, R_total/1e11, 'g--', linewidth=2, label='总阻力')
    ax5.set_xlabel('时间 (min)')
    ax5.set_ylabel('阻力 (×10¹¹ 1/m)')
    ax5.set_title('过滤阻力变化')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. 膜污染机理示意图
    ax6 = axes[1, 2]
    ax6.set_xlim(0, 10)
    ax6.set_ylim(0, 10)
    ax6.set_aspect('equal')
    ax6.axis('off')
    ax6.set_title('膜污染机理示意图')
    
    # 绘制膜
    membrane_rect = Rectangle((1, 2), 8, 0.3, facecolor='gray', edgecolor='black', linewidth=2)
    ax6.add_patch(membrane_rect)
    ax6.text(5, 1.5, '膜', ha='center', fontsize=10)
    
    # 污染物颗粒
    np.random.seed(42)
    for i in range(15):
        x = np.random.uniform(1.5, 8.5)
        y = np.random.uniform(3, 8)
        circle = Circle((x, y), 0.2, facecolor='brown', edgecolor='black', alpha=0.6)
        ax6.add_patch(circle)
    
    # 滤饼层
    cake = Rectangle((1, 2.3), 8, 0.5, facecolor='brown', edgecolor='black', alpha=0.5)
    ax6.add_patch(cake)
    ax6.text(5, 2.8, '滤饼层', ha='center', fontsize=9)
    
    # 流动方向
    ax6.arrow(5, 8.5, 0, -1.5, head_width=0.3, head_length=0.3, fc='blue', ec='blue')
    ax6.text(5, 9, '进水', ha='center', fontsize=9)
    ax6.arrow(5, 1.2, 0, -0.8, head_width=0.3, head_length=0.3, fc='green', ec='green')
    ax6.text(5, 0.5, '滤液', ha='center', fontsize=9)
    
    plt.tight_layout()
    plt.savefig(f'{OUTPUT_DIR}case2_membrane_filtration.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"案例2完成！图片已保存至 {OUTPUT_DIR}case2_membrane_filtration.png")
    print(f"初始通量: {J_clean:.1f} L/(m²·h)")
    print(f"1小时后滤饼过滤通量: {J_cake[-1]:.1f} L/(m²·h)")
    print(f"通量衰减率: {(1 - J_cake[-1]/J_clean)*100:.1f}%")
    print(f"累积滤液量: {V_cake[-1] * membrane['A'] * 1000:.1f} L")

# ==================== 案例3: 人工湿地污染物去除仿真 ====================
def case3_constructed_wetland():
    """案例3: 人工湿地污染物去除仿真"""
    print("\n" + "="*60)
    print("案例3: 人工湿地污染物去除仿真")
    print("="*60)
    
    # 湿地参数
    wetland = {
        'L': 50,         # 长度 (m)
        'W': 10,         # 宽度 (m)
        'H': 0.6,        # 深度 (m)
        'n': 0.4,        # 孔隙率
        'K': 0.01,       # 水力传导系数 (m/s)
        'Q': 100,        # 流量 (m³/d)
    }
    
    # 进水水质
    influent = {
        'COD': 200,      # mg/L
        'NH4': 40,       # mg/L
        'NO3': 0,        # mg/L
        'TP': 5,         # mg/L
    }
    
    # 去除速率常数 (1/d)
    k_rates = {
        'COD': 0.5,
        'NH4': 0.3,
        'NO3': 0.2,
        'TP': 0.15,
    }
    
    # 背景浓度 (mg/L)
    C_star = {
        'COD': 20,
        'NH4': 1,
        'NO3': 0.5,
        'TP': 0.2,
    }
    
    # 空间离散
    nx = 100
    x = np.linspace(0, wetland['L'], nx)
    dx = x[1] - x[0]
    
    # 水力负荷
    q = wetland['Q'] / (wetland['W'] * wetland['L'])  # m/d
    
    # 孔隙流速
    v = q / wetland['n']  # m/d
    
    # 停留时间
    tau = wetland['L'] / v  # d
    
    # k-C*模型求解
    def kCstar_model(x, C_in, k, C_star, v):
        """k-C*模型解析解"""
        return C_star + (C_in - C_star) * np.exp(-k * x / v)
    
    # 计算各污染物浓度分布
    C_COD = kCstar_model(x, influent['COD'], k_rates['COD'], C_star['COD'], v)
    C_NH4 = kCstar_model(x, influent['NH4'], k_rates['NH4'], C_star['NH4'], v)
    C_NO3 = C_star['NO3'] + (influent['NO3'] - C_star['NO3']) * np.exp(-k_rates['NO3'] * x / v)
    # 考虑硝化产生NO3
    C_NO3 += (influent['NH4'] - C_NH4) * 0.8  # 80%氨氮转化为硝酸盐
    C_TP = kCstar_model(x, influent['TP'], k_rates['TP'], C_star['TP'], v)
    
    # 去除效率
    removal_COD = (influent['COD'] - C_COD[-1]) / influent['COD'] * 100
    removal_NH4 = (influent['NH4'] - C_NH4[-1]) / influent['NH4'] * 100
    removal_TP = (influent['TP'] - C_TP[-1]) / influent['TP'] * 100
    
    # 停留时间分布 (串联CSTR模型)
    n_tanks = 3  # 等效串联级数
    t_rtd = np.linspace(0, 3*tau, 100)
    E = (n_tanks**n_tanks * t_rtd**(n_tanks-1)) / (math.factorial(n_tanks-1) * tau**n_tanks) * \
        np.exp(-n_tanks * t_rtd / tau)
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    fig.suptitle('案例3: 人工湿地污染物去除仿真', fontsize=14, fontweight='bold')
    
    # 1. COD去除
    ax1 = axes[0, 0]
    ax1.plot(x, C_COD, 'b-', linewidth=2)
    ax1.axhline(y=C_star['COD'], color='r', linestyle='--', label='背景浓度')
    ax1.set_xlabel('距离 (m)')
    ax1.set_ylabel('COD浓度 (mg/L)')
    ax1.set_title(f'COD去除 (去除率: {removal_COD:.1f}%)')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 氮素转化
    ax2 = axes[0, 1]
    ax2.plot(x, C_NH4, 'b-', linewidth=2, label='氨氮')
    ax2.plot(x, C_NO3, 'r-', linewidth=2, label='硝酸盐')
    ax2.plot(x, C_NH4 + C_NO3, 'g--', linewidth=2, label='总氮')
    ax2.set_xlabel('距离 (m)')
    ax2.set_ylabel('浓度 (mg/L)')
    ax2.set_title(f'氮素转化 (氨氮去除率: {removal_NH4:.1f}%)')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 总磷去除
    ax3 = axes[0, 2]
    ax3.plot(x, C_TP, 'b-', linewidth=2)
    ax3.axhline(y=C_star['TP'], color='r', linestyle='--', label='背景浓度')
    ax3.set_xlabel('距离 (m)')
    ax3.set_ylabel('TP浓度 (mg/L)')
    ax3.set_title(f'总磷去除 (去除率: {removal_TP:.1f}%)')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 停留时间分布
    ax4 = axes[1, 0]
    ax4.plot(t_rtd, E, 'b-', linewidth=2)
    ax4.axvline(x=tau, color='r', linestyle='--', label=f'理论停留时间 ({tau:.1f} d)')
    ax4.set_xlabel('停留时间 (d)')
    ax4.set_ylabel('E(t)')
    ax4.set_title('停留时间分布')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. 去除效率对比
    ax5 = axes[1, 1]
    pollutants = ['COD', 'NH4-N', 'TP']
    removals = [removal_COD, removal_NH4, removal_TP]
    colors = ['blue', 'green', 'orange']
    bars = ax5.bar(pollutants, removals, color=colors, alpha=0.7, edgecolor='black')
    ax5.set_ylabel('去除效率 (%)')
    ax5.set_title('污染物去除效率对比')
    ax5.set_ylim(0, 100)
    for bar, rem in zip(bars, removals):
        ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
                f'{rem:.1f}%', ha='center', fontsize=10)
    ax5.grid(True, alpha=0.3, axis='y')
    
    # 6. 人工湿地剖面示意图
    ax6 = axes[1, 2]
    ax6.set_xlim(0, 10)
    ax6.set_ylim(0, 10)
    ax6.set_aspect('equal')
    ax6.axis('off')
    ax6.set_title('水平潜流人工湿地示意图')
    
    # 湿地床体
    wetland_bed = Rectangle((1, 3), 8, 3, facecolor='sandybrown', edgecolor='black', linewidth=2)
    ax6.add_patch(wetland_bed)
    
    # 填料层
    for i in range(20):
        for j in range(6):
            circle = Circle((1.2 + i*0.38, 3.2 + j*0.45), 0.15, 
                          facecolor='gray', edgecolor='black', alpha=0.5)
            ax6.add_patch(circle)
    
    # 水面
    water = Rectangle((1, 4.5), 8, 1.5, facecolor='lightblue', alpha=0.5)
    ax6.add_patch(water)
    
    # 进水管
    ax6.plot([0.5, 1], [5.25, 5.25], 'b-', linewidth=4)
    ax6.arrow(0.5, 5.25, 0.3, 0, head_width=0.2, head_length=0.1, fc='blue', ec='blue')
    ax6.text(0.3, 5.6, '进水', fontsize=9)
    
    # 出水管
    ax6.plot([9, 9.5], [5.25, 5.25], 'g-', linewidth=4)
    ax6.arrow(9.2, 5.25, 0.3, 0, head_width=0.2, head_length=0.1, fc='green', ec='green')
    ax6.text(9.2, 5.6, '出水', fontsize=9)
    
    # 植物
    for i in range(5):
        x_plant = 2 + i * 1.5
        # 茎
        ax6.plot([x_plant, x_plant], [6, 7.5], 'g-', linewidth=3)
        # 叶
        ax6.plot([x_plant, x_plant-0.3], [7, 7.3], 'g-', linewidth=2)
        ax6.plot([x_plant, x_plant+0.3], [7, 7.3], 'g-', linewidth=2)
        ax6.plot([x_plant, x_plant], [7, 7.5], 'g-', linewidth=2)
    
    ax6.text(5, 2.5, '砾石填料', ha='center', fontsize=9)
    
    plt.tight_layout()
    plt.savefig(f'{OUTPUT_DIR}case3_constructed_wetland.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"案例3完成！图片已保存至 {OUTPUT_DIR}case3_constructed_wetland.png")
    print(f"水力负荷: {q*1000:.1f} mm/d")
    print(f"停留时间: {tau:.1f} d")
    print(f"出水COD: {C_COD[-1]:.1f} mg/L (去除率{removal_COD:.1f}%)")
    print(f"出水氨氮: {C_NH4[-1]:.1f} mg/L (去除率{removal_NH4:.1f}%)")
    print(f"出水总磷: {C_TP[-1]:.1f} mg/L (去除率{removal_TP:.1f}%)")

# ==================== 案例4: 地下水污染物迁移扩散 ====================
def case4_groundwater_contamination():
    """案例4: 地下水污染物迁移扩散"""
    print("\n" + "="*60)
    print("案例4: 地下水污染物迁移扩散")
    print("="*60)
    
    # 含水层参数
    aquifer = {
        'K': 1e-4,       # 水力传导系数 (m/s)
        'n': 0.3,        # 孔隙度
        'alpha_L': 10,   # 纵向弥散度 (m)
        'alpha_T': 1,    # 横向弥散度 (m)
        'D_m': 1e-9,     # 分子扩散系数 (m²/s)
    }
    
    # 水力梯度
    dh_dx = 0.001
    
    # 达西流速
    q = aquifer['K'] * dh_dx  # m/s
    
    # 孔隙流速
    v = q / aquifer['n']  # m/s
    v_day = v * 86400  # m/d
    
    print(f"达西流速: {q*1000*86400:.2f} mm/d")
    print(f"孔隙流速: {v_day:.2f} m/d")
    
    # 弥散系数
    D_L = aquifer['D_m'] + aquifer['alpha_L'] * v  # 纵向
    D_T = aquifer['D_m'] + aquifer['alpha_T'] * v  # 横向
    
    print(f"纵向弥散系数: {D_L:.2e} m²/s")
    print(f"横向弥散系数: {D_T:.2e} m²/s")
    
    # 污染源参数
    source = {
        'x0': 0,         # 污染源位置
        'y0': 0,
        'C0': 1000,      # 初始浓度 (mg/L)
        'M': 1000,       # 污染物质量 (kg)
        't_release': 30, # 释放时间 (d)
    }
    
    # 计算网格
    nx, ny = 100, 50
    x = np.linspace(-50, 450, nx)  # m
    y = np.linspace(-50, 50, ny)   # m
    X, Y = np.meshgrid(x, y)
    
    # 时间步
    times = [10, 30, 90, 365]  # d
    
    # 解析解：二维对流-弥散方程（点源）
    def advection_dispersion_2D(X, Y, t, C0, x0, y0, v, D_L, D_T, n):
        """二维对流-弥散方程解析解"""
        t_sec = t * 86400  # 转换为秒
        
        # 浓度分布
        C = np.zeros_like(X)
        
        # 只计算t>0的区域
        if t > 0:
            # 解析解公式
            exponent = -((X - x0 - v*t_sec)**2) / (4*D_L*t_sec) - (Y - y0)**2 / (4*D_T*t_sec)
            C = C0 / (4 * np.pi * t_sec * np.sqrt(D_L * D_T)) * np.exp(exponent)
            
            # 考虑孔隙度
            C = C / n
        
        return C
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    fig.suptitle('案例4: 地下水污染物迁移扩散', fontsize=14, fontweight='bold')
    
    # 浓度分布随时间变化
    for idx, (ax, t) in enumerate(zip(axes.flat, times)):
        C = advection_dispersion_2D(X, Y, t, source['C0'], source['x0'], source['y0'], 
                                     v, D_L, D_T, aquifer['n'])
        
        # 限制最大浓度以便显示
        C_plot = np.clip(C, 0, 100)
        
        # 绘制浓度分布
        contour = ax.contourf(X, Y, C_plot, levels=20, cmap='jet', alpha=0.8)
        ax.contour(X, Y, C_plot, levels=5, colors='black', alpha=0.5, linewidths=0.5)
        
        # 标记污染源
        ax.plot(source['x0'], source['y0'], 'r*', markersize=15, label='污染源')
        
        # 添加流线
        y_stream = np.linspace(-40, 40, 9)
        for y_s in y_stream:
            ax.annotate('', xy=(400, y_s), xytext=(0, y_s),
                       arrowprops=dict(arrowstyle='->', color='white', alpha=0.5, lw=1))
        
        ax.set_xlabel('距离 (m)')
        ax.set_ylabel('距离 (m)')
        ax.set_title(f't = {t} d')
        ax.set_aspect('equal')
        ax.legend()
        
        # 添加颜色条
        cbar = plt.colorbar(contour, ax=ax, shrink=0.6)
        cbar.set_label('浓度 (mg/L)')
    
    plt.tight_layout()
    plt.savefig(f'{OUTPUT_DIR}case4_groundwater_contamination.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    # 创建浓度随时间和距离变化图
    fig2, axes2 = plt.subplots(1, 2, figsize=(14, 5))
    fig2.suptitle('案例4: 污染物浓度时空分布', fontsize=14, fontweight='bold')
    
    # 1. 中心线浓度随时间变化
    ax1 = axes2[0]
    x_center = np.linspace(0, 400, 200)
    colors = plt.cm.viridis(np.linspace(0, 1, len(times)))
    
    for t, color in zip(times, colors):
        C_center = advection_dispersion_2D(x_center, 0, t, source['C0'], source['x0'], source['y0'], 
                                           v, D_L, D_T, aquifer['n'])
        ax1.plot(x_center, C_center, color=color, linewidth=2, label=f't = {t} d')
    
    ax1.set_xlabel('距离 (m)')
    ax1.set_ylabel('浓度 (mg/L)')
    ax1.set_title('中心线浓度随时间变化')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim(0, 100)
    
    # 2. 固定位置浓度随时间变化
    ax2 = axes2[1]
    t_range = np.linspace(1, 365, 100)
    x_obs = [50, 100, 200, 300]  # 观测井位置
    
    for x_o in x_obs:
        C_time = []
        for t in t_range:
            C = advection_dispersion_2D(x_o, 0, t, source['C0'], source['x0'], source['y0'], 
                                        v, D_L, D_T, aquifer['n'])
            C_time.append(C)
        ax2.plot(t_range, C_time, linewidth=2, label=f'x = {x_o} m')
    
    ax2.set_xlabel('时间 (d)')
    ax2.set_ylabel('浓度 (mg/L)')
    ax2.set_title('固定位置浓度随时间变化')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_ylim(0, 100)
    
    plt.tight_layout()
    plt.savefig(f'{OUTPUT_DIR}case4_groundwater_temporal.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"案例4完成！图片已保存至 {OUTPUT_DIR}case4_groundwater_contamination.png")
    print(f"365天时最大浓度: {np.max(C):.1f} mg/L")
    print(f"污染羽前锋位置: 约 {v_day * 365:.0f} m")

# ==================== 创建动画 ====================
def create_wetland_animation():
    """创建人工湿地污染物去除动画"""
    print("\n" + "="*60)
    print("创建人工湿地污染物去除动画...")
    print("="*60)
    
    # 湿地参数
    wetland = {
        'L': 50,         # 长度 (m)
        'W': 10,         # 宽度 (m)
        'n': 0.4,        # 孔隙率
        'Q': 100,        # 流量 (m³/d)
    }
    
    influent = {'COD': 200, 'NH4': 40, 'TP': 5}
    k_rates = {'COD': 0.5, 'NH4': 0.3, 'TP': 0.15}
    C_star = {'COD': 20, 'NH4': 1, 'TP': 0.2}
    
    q = wetland['Q'] / (wetland['W'] * wetland['L'])
    v = q / wetland['n']
    
    # 空间网格
    nx = 100
    x = np.linspace(0, wetland['L'], nx)
    
    # 时间序列（模拟动态过程）
    nt = 50
    times = np.linspace(0, 10, nt)  # 10天
    
    # 创建图形
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    def update(frame):
        for ax in axes:
            ax.clear()
        
        t = times[frame]
        
        # 计算当前时刻的浓度分布（考虑时间延迟）
        x_eff = np.maximum(0, x - v * t)
        
        C_COD = C_star['COD'] + (influent['COD'] - C_star['COD']) * np.exp(-k_rates['COD'] * x_eff / v)
        C_NH4 = C_star['NH4'] + (influent['NH4'] - C_star['NH4']) * np.exp(-k_rates['NH4'] * x_eff / v)
        C_TP = C_star['TP'] + (influent['TP'] - C_star['TP']) * np.exp(-k_rates['TP'] * x_eff / v)
        
        # 污染物前锋位置
        front_pos = min(v * t, wetland['L'])
        
        # 1. COD
        axes[0].fill_between(x, 0, C_COD, alpha=0.5, color='blue')
        axes[0].plot(x, C_COD, 'b-', linewidth=2)
        axes[0].axvline(x=front_pos, color='r', linestyle='--', alpha=0.7)
        axes[0].set_xlabel('距离 (m)')
        axes[0].set_ylabel('COD浓度 (mg/L)')
        axes[0].set_title(f'COD去除 (t = {t:.1f} d)')
        axes[0].set_ylim(0, 220)
        axes[0].grid(True, alpha=0.3)
        
        # 2. 氨氮
        axes[1].fill_between(x, 0, C_NH4, alpha=0.5, color='green')
        axes[1].plot(x, C_NH4, 'g-', linewidth=2)
        axes[1].axvline(x=front_pos, color='r', linestyle='--', alpha=0.7)
        axes[1].set_xlabel('距离 (m)')
        axes[1].set_ylabel('氨氮浓度 (mg/L)')
        axes[1].set_title(f'氨氮去除 (t = {t:.1f} d)')
        axes[1].set_ylim(0, 45)
        axes[1].grid(True, alpha=0.3)
        
        # 3. 总磷
        axes[2].fill_between(x, 0, C_TP, alpha=0.5, color='orange')
        axes[2].plot(x, C_TP, 'orange', linewidth=2)
        axes[2].axvline(x=front_pos, color='r', linestyle='--', alpha=0.7)
        axes[2].set_xlabel('距离 (m)')
        axes[2].set_ylabel('TP浓度 (mg/L)')
        axes[2].set_title(f'总磷去除 (t = {t:.1f} d)')
        axes[2].set_ylim(0, 6)
        axes[2].grid(True, alpha=0.3)
        
        plt.tight_layout()
        
        return axes
    
    # 创建动画
    from matplotlib.animation import FuncAnimation
    
    anim = FuncAnimation(fig, update, frames=nt, interval=200, blit=False)
    
    # 保存为GIF
    anim.save(f'{OUTPUT_DIR}wetland_pollutant_animation.gif', writer='pillow', fps=5, dpi=100)
    plt.close()
    
    print(f"动画创建完成！已保存至 {OUTPUT_DIR}wetland_pollutant_animation.gif")

# ==================== 主程序 ====================
if __name__ == "__main__":
    print("="*60)
    print("主题054：水处理与环境工程仿真")
    print("Water Treatment and Environmental Engineering Simulation")
    print("="*60)
    
    print("\n本程序演示水处理过程的仿真方法，包含以下4个案例：")
    print("1. 活性污泥反应器仿真（ASM简化模型）")
    print("2. 膜分离过程与膜污染分析")
    print("3. 人工湿地污染物去除仿真")
    print("4. 地下水污染物迁移扩散")
    
    # 运行案例
    case1_activated_sludge()
    case2_membrane_filtration()
    case3_constructed_wetland()
    case4_groundwater_contamination()
    
    # 创建动画
    print("\n" + "="*60)
    print("创建水处理过程动画...")
    print("="*60)
    create_wetland_animation()
    
    print("\n" + "="*60)
    print("所有案例完成！")
    print("="*60)
    print("\n输出文件:")
    print("  - case1_activated_sludge.png")
    print("  - case2_membrane_filtration.png")
    print("  - case3_constructed_wetland.png")
    print("  - case4_groundwater_contamination.png")
    print("  - case4_groundwater_temporal.png")
    print("  - wetland_pollutant_animation.gif")

在这里插入图片描述