简介面向环境科学与水利工程领域的MATLAB源码包针对河流、渠道等一维流动水体水质变化模拟而编写采用一维差分法求解连续性方程与质量守恒原理。压缩包内共2个.m脚本文件整体仅1KB体量轻、易阅读适合作为一维水质建模与数值计算入门的教学示例或科研模板。已有285人学习浏览简洁直接的实现思路对同类需求有参考价值。程序覆盖离散网格划分、多种差分格式选择、边界条件设定、时间步长控制与逐断面迭代求解等关键环节能够输出水质参数随时间和空间的变化曲线帮助用户完整掌握一维水质模型从方程离散到结果可视化的过程。在此基础上可结合溶解氧、氨氮等实际指标修改参数用于污染预测、水环境容量评估及水资源管理对策制定。1. 一维水质模型什么时候值得把河段简化成一维一个沿江排污口连续排放含氨氮废水下游三个监测断面浓度持续超标。你想知道削减多少负荷、在哪里削减最有效一维水质模型是环评和水环境容量计算里最常用的工具。所谓一维是指污染物浓度在河流断面上已混合均匀只沿河长方向变化——对中小河段而言这比二维、三维模型简单得多数据需求也低得多。一维差分法的任务就是把描述对流和弥散的控制方程在网格上拆成代数方程让浓度分布不再是解析失真的曲线而是能贴近真实边界条件的数值解。这篇文章写给环境工程、水利水务和做模型开发的从业者会从方程讲到显式和隐式差分再到网格步长的取舍和代码验证。2. 一维水质模型的控制方程与一维差分法离散化2.1 对流-弥散-衰减方程里每一项的实际含义常用的一个一维水质模型控制方程是∂C/∂t -u·∂C/∂x D·∂²C/∂x² - k·C S其中 C 是断面平均浓度单位常用 mg/Lu 是断面平均流速m/sD 是纵向弥散系数m²/sk 是一级衰减系数s⁻¹S 是源汇项代表单位时间、单位体积内由排污口、支流汇入或取水口引起的浓度变化单位 mg/(L·s)。对流项描述水流拉着污染物向下游走弥散项描述流速非均匀和湍流造成的纵向展开衰减项挂钩生化降解S 则体现外部输入输出。实际建模中最容易翻车的是单位换算衰减系数从 d⁻¹ 换成 s⁻¹ 要除以 86400源项写成 kg/d 时要换算成 mg/s再折算进单个网格的水体体积。每个新项目我都会重新查一遍这套换算而不是沿用旧脚本。横向扩散项被略掉背后的工程假设是“断面已充分混合”。工程上常用完全混合长度 L 0.4·u·B²/D_t 估算这个前提是否成立其中 B 是河宽D_t 是横向扩散系数。排污口下游的建模断面离源头的距离不够这个长度时简化成一维会让 D 吸收掉横向未混合的效应标定出的弥散系数虚高。我一般要求建模断面距离排污口不小于完全混合长度的 1.5 倍再谈一维建模。2.2 空间导数怎么差分迎风与中心差分怎么组合把空间导数替换成差商是一维差分法区别于解析解的关键一步。对流项用中心差分 (C_{i1} - C_{i-1})/(2Δx) 有二阶精度但在对流占优区域容易振荡工程代码里更常见的是迎风差分上游差分∂C/∂x ≈ (C_i - C_{i-1})/Δx迎风差分只取上游方向的信息牺牲了一点精度但换来单调无振荡的解。弥散项固定用二阶中心差分∂²C/∂x² ≈ (C_{i1} - 2C_i C_{i-1})/Δx²两者组合后半离散方程变成一组常微分方程dC_i/dt -u·(C_i - C_{i-1})/Δx D·(C_{i1} - 2C_i C_{i-1})/Δx² - k·C_i S_i对流项用迎风、弥散项用中心这种混合离散的精度并不统一但我在一维河段模型里很少用中心差分处理对流项网格 Peclet 数偏大时中心差分产生的负浓度比那一点精度损失难解释得多。2.3 显式、隐式、Crank-Nicolson:三套时间推进的取舍时间方向用 θ 加权统一描述C^{n1} C^n Δt·[θ·F(C^{n1}) (1-θ)·F(C^n)]。θ0 是显式向前欧拉θ1 是全隐式向后欧拉θ0.5 是 Crank-Nicolson。显式格式每个网格独立推进写起来最顺手但稳定性条件苛刻对流项要求 Courant 数 Cr uΔt/Δx ≤ 1弥散项要求 DΔt/Δx² ≤ 0.5。实际河段模型中弥散项的限制往往更严导致时间步长比直觉小一个量级。全隐式无稳定条件限制时间步可以放大但一阶时间精度在大步长下有明显数值耗散浓度峰会被人为拉平。Crank-Nicolson 精度和稳定性折中最好代价是实现多一层矩阵处理。格式时间精度稳定性限制数值耗散计算量显式一阶Cr ≤ 1 且 DΔt/Δx² ≤ 0.5小最小全隐式一阶无条件稳定较大解三对角方程组Crank-Nicolson二阶极轻微振荡可能小解三对角方程组我自己做主长期稳态水环境容量计算时会用隐式事故性瞬时泄漏场景需要精确刻画浓度峰则偏向 Crank-Nicolson或者把隐式的时间步压到峰值时间尺度的十分之一以下。2.4 一维差分法的三对角系数矩阵怎么组装无论隐式还是 Crank-Nicolson最终都落在线性方程组 A·C^{n1} b 上。A 的主体是三对角形式主对角线来自时间项、弥散项的对角部分和衰减项上对角线来自弥散项的 i1 系数下对角线来自迎风对流项和弥散项的 i-1 系数。以下用全隐式格式演示import numpy as np N 400 # 网格数 dx 25.0 # 空间步长m dt 120.0 # 时间步长s u 0.5 # 流速m/s D 15.0 # 纵向弥散系数m²/s k 1.0 / 86400 # 衰减系数s⁻¹ A np.zeros((N, N)) r D * dt / dx**2 # 弥散数 p u * dt / dx # Courant 数 for i in range(1, N - 1): A[i, i - 1] p r A[i, i] 1.0 2.0 * r k * dt A[i, i 1] -r # 上游 Dirichlet 边界第 0 个网格浓度固定为入流值 A[0, 0] 1.0 # 下游 Neumann 边界C_N C_{N-1}实现为零梯度出流 A[N - 1, N - 2] -1.0 A[N - 1, N - 1] 1.0弥散数 r 是 DΔt/Δx²Courant 数 p 是 uΔt/Δx。主对角线上的 1.0 来自时间项2r 是弥散项的中心对角贡献kΔt 是衰减项。对流项迎风离散在下对角线贡献 p位置取决于水流方向水流反向时这一项要挪到上对角线。下游边界写成 C_{N} C_{N-1}相当于零浓度梯度出流比直接设 0 合理得多后面章节会专门讲这个坑。3. 用 Python 实现一维水质模型显式与隐式两版代码3.1 模型场景、单位统一和网格尺寸为一个 10 km 长的河段建模流速 0.5 m/s纵向弥散系数 15 m²/s氨氮一级衰减系数按 0.3 d⁻¹ 取值。x3 km 处有连续排污口排放浓度 20 mg/L排放流量 0.5 m³/s上游来流和本底浓度均为 1.5 mg/L。先统一单位k 0.3 / 86400 3.47×10⁻⁶ s⁻¹。河流断面积按流量 20 m³/s 反推取 40 m²。网格尺寸按 Peclet 约束来定。网格 Peclet 数 Pe_grid uΔx/D 要小于 1D15 m²/s 时 Δx 30 m取 25 m网格数 N400。显式格式的时间步要同时满足 Cr ≤ 1 和 DΔt/Δx² ≤ 0.5Cr 条件给出 Δt 50 s弥散条件给出 Δt 20.8 s所以显式取 Δt10 s。隐式格式没有稳定性限制取 Δt120 s模拟同样时长耗时只有显式的约 1/12。参数值说明河长 L10000 m计算域长度网格数 N400对应 Δx25 m网格步长 Δx25 m由 Pe_grid 1 约束显式时间步10 s受 DΔt/Δx² ≤ 0.5 限制隐式时间步120 s无稳定限制看精度需求流速 u0.5 m/s断面平均弥散系数 D15 m²/s经验估算起点衰减系数 k3.47×10⁻⁶ s⁻¹0.3 d⁻¹ 换算排口位置x3000 m网格序号 120背景浓度1.5 mg/L初始与上游边界值3.2 显式差分的最小可运行代码显式格式不需要解方程组直接从旧浓度推出新浓度。边界处理和源项注入写清楚代码就能跑通import numpy as np L, N 10000.0, 400 dx L / N u, D 0.5, 15.0 k 0.3 / 86400.0 # 单位统一为 s^-1 dt 10.0 # 显式步长满足稳定性条件 cr u * dt / dx # Courant 数 dr D * dt / dx**2 # 弥散条件参数 assert cr 1.0 and dr 0.5 C np.full(N, 1.5) # 初始浓度 mg/L src_node 120 # 排口对应网格序号 src_rate 0.5 * 20.0 / 40.0 * (dt / dx) # 每步浓度抬升量 for step in range(7200): # 显式模拟 20 小时 C_new C.copy() for i in range(1, N - 1): advect -u * (C[i] - C[i - 1]) / dx diff D * (C[i 1] - 2.0 * C[i] C[i - 1]) / dx**2 C_new[i] C[i] dt * (advect diff - k * C[i]) C_new[0] 1.5 # 上游入流浓度固定 C_new[-1] C_new[-2] # 下游零梯度出流 C_new[src_node] src_rate C C_new np.save(concentration_explicit.npy, C)代码里 src_rate 的折算逻辑是排污负荷 0.5 m³/s×20 mg/L10 mg/s除过水断面积 40 m²再乘时间步 dt就是该网格在一个时间步内应当抬高的浓度。这种折算把单网格水体体积隐含在断面积和步长里适合均匀网格网格尺度差异大的工程模型要按每个网格的实际水体体积 V A·Δx 重新折算。这段代码跑完排口下游会看到浓度从 1.5 mg/L 抬升后逐渐向下游推进并衰减。3.3 隐式差分与 Thomas 追赶法隐式格式每个时间步要解三对角方程组。Thomas 追赶法复杂度 O(N)适合这类一维问题不依赖额外的求解器import numpy as np a np.full(N, u * dt / dx D * dt / dx**2) # 下对角线 b np.full(N, 1.0 2.0 * D * dt / dx**2 k * dt) # 主对角线 c np.full(N, -D * dt / dx**2) # 上对角线 b[0] 1.0; c[0] 0.0 # 上游 Dirichlet a[-1] -1.0; b[-1] 1.0 # 下游 Neumann for step in range(3000): # 隐式模拟 100 小时dt120s R C.copy() R[0] 1.5 # 上游边界值 R[src_node] 0.5 * 20.0 / 40.0 * (dt / dx) # Thomas 追赶法前代 cp np.zeros(N); dp np.zeros(N) cp[0] -c[0] / b[0] dp[0] R[0] / b[0] for i in range(1, N): denom b[i] a[i] * cp[i - 1] cp[i] -c[i] / denom dp[i] (R[i] - a[i] * dp[i - 1]) / denom # 后代回代 C_new np.zeros(N) C_new[-1] dp[-1] for i in range(N - 2, -1, -1): C_new[i] dp[i] cp[i] * C_new[i 1] C C_new np.save(concentration_implicit.npy, C)a、b、c 分别是三对角矩阵的下对角线、主对角线和上对角线。迎风对流项放在下对角线决定了信息从上游向下游传递。下游 Neumann 边界在 Thomas 算法里表现为 R 的最后一项保持 0因为 C_{N} C_{N-1} 时该行方程右端为零。把这套代码和显式版跑同样的 Δt 对比显式会出现数值振荡因为 Δt 超出了弥散稳定条件隐式解则稳定吸收掉了高频扰动但峰值会被略微压低这就是隐式格式数值耗散的直接体现。4. 参数取值、网格约束与边界条件一维水质模型落地细节4.1 弥散系数 D 和衰减系数 k 没有实测时怎么估D 最常见的初值来自 D α·uα 是纵向弥散度。山区陡坡小河流 α 约 2-10 m中下游平原河流可到 30-60 m。取 α30 m、u0.5 m/s得到 D15 m²/s这是我在无实测资料时的标准起点。衰减系数按污染物类型先给初值再靠实测修正污染物k 典型范围 (d⁻¹)备注BOD₅0.15-0.35水温越高衰减越快COD0.02-0.10难降解组分占主导NH₃-N0.05-0.50硝化作用为主总磷0.005-0.05沉降与底泥吸附为主用上下游两个实测断面校准时我一般先用简化稳态对流衰减式 C(x)C₀·exp(-kx/u) 粗标定 k再通过浓度曲线的纵向展宽调 D。标定顺序必须是先衰减后弥散否则 k 和 D 高度耦合参数不唯一。校准迭代两轮后再看排口下游细节偏差微调源项。4.2 网格 Peclet 数和 Courant 数动手前先算的两个数网格 Peclet 数 Pe_grid uΔx/D 决定空间离散的成败。Pe_grid 2 时即使隐式格式时间上稳定空间上也会出现非物理振荡或过度数值耗散浓度峰后跟着一个下冲。Courant 数 Cr uΔt/Δx 则掌控时间推进中对流信息的传播速度。显式格式 Cr 1 会直接不收敛隐式格式虽然没有这个硬约束但 Cr 太大时解的相位偏差会变大峰值位置偏移以公里计。选网格的顺序是先根据河段长度和需要分辨的浓度梯度定初始 Δx算 Pe_grid不满足就加密再用格式类型定 Δt。比如 u1 m/s、D10 m²/s 时Δx 只能取到 10 m 才能使 Pe_grid ≤ 1。如果不想加密网格单纯靠隐式格式放大时间步结果是解看起来稳定但浓度峰位置和幅度都失真。隐式格式解决的是时间步长稳定性解决不了空间离散的震荡问题。这两个数我每次建模都会写在脚本第一段的注释里方便复查。4.3 点源、支流汇入和边界条件的 3 个常见坑第一个坑是下游边界设成零浓度。有些初版代码把 A[-1,-1]1.0 且右端项 R[-1]0.0等于强制边界节点浓度为零边界瞬间变成一个吸收阱上游整段浓度都会被拉低。正确做法是零梯度出流即 C_N C_{N-1}在矩阵里表现为最后一行的系数为 -1 和 1。第二个坑是点源加的时间位置不对。连续排污应作为源项 S 叠加在每个时间步的右端向量上而不是写进初始条件事故性瞬时排放才应该在初始浓度里一次性置入置入量按总排放质量除以网格水体体积计算。二者搞反排口浓度会随时间线性累积或瞬间流失。第三个坑是时间单位混用。代码里全部用秒单位是唯一不会出错的方案。衰减系数、流量、排放速率、模拟总时长任何一项用了天或小时都要在进入方程前统一换算。现象是浓度偏置好几倍但数值解形态完全正常这类 bug 最难定位。常见坑现象处理下游边界零浓度边界附近浓度骤降上游偏低换成零梯度出流连续源写成初值排口浓度逐时步累积连续源加在右端向量瞬时源加初值单位混用 day 与 s衰减快慢差 86400 倍全模型统一用 SI 秒单位Pe_grid 过大浓度出现负值或锯齿加密 Δx 或改高阶格式5. 用解析解和质量守恒给差分代码“验算”5.1 稳态解析解做基准有衰减、忽略弥散的稳态连续源条件下浓度沿程满足 C(x) C₁ (C₀ - C₁)·exp(-k·x/u)这是一个可以直接手算的解析基准。我会把模型跑到数百小时后取稳态剖面与解析解比较。相对误差控制在 5% 以内算通过误差集中在排口附近多是源项折算或网格编码偏差误差整体沿程放大几乎可以肯定是 k 或单位用错。x np.linspace(0, L, N) c_analytic 1.5 (20.0 - 1.5) * np.exp(-k * (x - 3000.0) / u) mask x 3000.0 err np.max(np.abs(C[mask] - c_analytic[mask]) / c_analytic[mask]) print(f最大相对误差: {err:.2%})这段代码只适用于排口下游且弥散影响较小的区段。如果整段误差都大先检查 k 的秒/天换算如果只在峰后出现局部误差再回头检查网格 Peclet 数。5.2 质量守恒核算每一步都该稳得住无论代码来自现成的参考包还是自己重写质量守恒都是一维水质模型的第一道验收。核算式是入流质量 源项质量 出流质量 衰减质量 系统增量。稳定运行后最后一项应趋近零。核算方法累计每个时刻的 u·C_in·A·dt 与 u·C_out·A·dt再累计 k·ΣC·A·dx·dt。三者误差超过 2% 就不该继续调参先回头检查矩阵边界行和源项折算。比如你在资料包里拿到 shuizhi.zip 之类的现成实现第一件事不是换参数而是跑一组已知解析解的工况确认守恒误差落在 2% 以内。5.3 浓度出现负值或锯齿状振荡怎么办按固定顺序排查避免乱改参数计算 Pe_grid确认空间步长是否满足 uΔx/D ≤ 1显式格式核对 Cr 和 DΔt/Δx² 两个稳定条件隐式格式检查 Δt 是否大到使单步浓度变化超过峰值的 10%观察振荡位置若只出现在排口下游两三个网格优先加密排口附近的 Δx。修复操作与排查顺序严格对应Pe_grid 超限就加密网格显式超稳定条件就缩小 Δt隐式步长过大就按峰值时间尺度重新选 Δt排口局部振荡就在源项附近做局部网格加密。修改完成后重新运行 5.1 的解析解对比确认相对误差回到 5% 以内再看浓度曲线是否保持单调下降。本文还有配套的精品资源点击获取