企业官网建设流程全解析

简介：面向环境科学与水利工程领域的MATLAB源码包，针对河流、渠道等一维流动水体水质变化模拟而编写，采用一维差分法求解连续性方程与质量守恒原理。压缩包内共2个.m脚本文件，整体仅1KB，体量轻、易阅读，适合作为一维水质建模与数值计算入门的教学示例或科研模板。已有285人学习浏览，简洁直接的实现思路对同类需求有参考价值。程序覆盖离散网格划分、多种差分格式选择、边界条件设定、时间步长控制与逐断面迭代求解等关键环节，能够输出水质参数随时间和空间的变化曲线，帮助用户完整掌握一维水质模型从方程离散到结果可视化的过程。在此基础上，可结合溶解氧、氨氮等实际指标修改参数，用于污染预测、水环境容量评估及水资源管理对策制定。

1. 一维水质模型：什么时候值得把河段简化成一维

一个沿江排污口连续排放含氨氮废水，下游三个监测断面浓度持续超标。你想知道削减多少负荷、在哪里削减最有效，一维水质模型是环评和水环境容量计算里最常用的工具。所谓一维，是指污染物浓度在河流断面上已混合均匀，只沿河长方向变化——对中小河段而言，这比二维、三维模型简单得多，数据需求也低得多。一维差分法的任务，就是把描述对流和弥散的控制方程在网格上拆成代数方程，让浓度分布不再是解析失真的曲线，而是能贴近真实边界条件的数值解。这篇文章写给环境工程、水利水务和做模型开发的从业者，会从方程讲到显式和隐式差分，再到网格步长的取舍和代码验证。

2. 一维水质模型的控制方程与一维差分法离散化

2.1 对流-弥散-衰减方程里每一项的实际含义

常用的一个一维水质模型控制方程是：

∂C/∂t = -u·∂C/∂x + D·∂²C/∂x² - k·C + S

其中 C 是断面平均浓度，单位常用 mg/L；u 是断面平均流速，m/s；D 是纵向弥散系数，m²/s；k 是一级衰减系数，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_{i+1} - C_{i-1})/(2Δx) 有二阶精度，但在对流占优区域容易振荡；工程代码里更常见的是迎风差分（上游差分）：

∂C/∂x ≈ (C_i - C_{i-1})/Δx

迎风差分只取上游方向的信息，牺牲了一点精度，但换来单调无振荡的解。弥散项固定用二阶中心差分：

∂²C/∂x² ≈ (C_{i+1} - 2C_i + C_{i-1})/Δx²

两者组合后，半离散方程变成一组常微分方程：

dC_i/dt = -u·(C_i - C_{i-1})/Δx + D·(C_{i+1} - 2C_i + C_{i-1})/Δx² - k·C_i + S_i

对流项用迎风、弥散项用中心，这种混合离散的精度并不统一，但我在一维河段模型里很少用中心差分处理对流项：网格 Peclet 数偏大时，中心差分产生的负浓度比那一点精度损失难解释得多。

2.3 显式、隐式、Crank-Nicolson:三套时间推进的取舍

时间方向用 θ 加权统一描述：C^{n+1} = C^n + Δt·[θ·F(C^{n+1}) + (1-θ)·F(C^n)]。θ=0 是显式（向前欧拉），θ=1 是全隐式（向后欧拉），θ=0.5 是 Crank-Nicolson。显式格式每个网格独立推进，写起来最顺手，但稳定性条件苛刻：对流项要求 Courant 数 Cr = uΔt/Δx ≤ 1，弥散项要求 DΔt/Δx² ≤ 0.5。实际河段模型中，弥散项的限制往往更严，导致时间步长比直觉小一个量级。全隐式无稳定条件限制，时间步可以放大，但一阶时间精度在大步长下有明显数值耗散，浓度峰会被人为拉平。Crank-Nicolson 精度和稳定性折中最好，代价是实现多一层矩阵处理。

我自己做主长期稳态水环境容量计算时会用隐式；事故性瞬时泄漏场景需要精确刻画浓度峰，则偏向 Crank-Nicolson，或者把隐式的时间步压到峰值时间尺度的十分之一以下。

2.4 一维差分法的三对角系数矩阵怎么组装

无论隐式还是 Crank-Nicolson，最终都落在线性方程组 A·C^{n+1} = b 上。A 的主体是三对角形式：主对角线来自时间项、弥散项的对角部分和衰减项；上对角线来自弥散项的 i+1 系数；下对角线来自迎风对流项和弥散项的 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⁻¹ 取值。x=3 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 要小于 1，D=15 m²/s 时 Δx < 30 m，取 25 m，网格数 N=400。显式格式的时间步要同时满足 Cr ≤ 1 和 DΔt/Δx² ≤ 0.5：Cr 条件给出 Δt < 50 s，弥散条件给出 Δt < 20.8 s，所以显式取 Δt=10 s。隐式格式没有稳定性限制，取 Δt=120 s，模拟同样时长耗时只有显式的约 1/12。

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/L=10 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 小时，dt=120s 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、u=0.5 m/s，得到 D=15 m²/s，这是我在无实测资料时的标准起点。衰减系数按污染物类型先给初值，再靠实测修正：

用上下游两个实测断面校准时，我一般先用简化稳态对流衰减式 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。比如 u=1 m/s、D=10 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 最难定位。

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 浓度出现负值或锯齿状振荡怎么办

按固定顺序排查，避免乱改参数：

1. 计算 Pe_grid，确认空间步长是否满足 uΔx/D ≤ 1；
2. 显式格式核对 Cr 和 DΔt/Δx² 两个稳定条件；
3. 隐式格式检查 Δt 是否大到使单步浓度变化超过峰值的 10%；
4. 观察振荡位置，若只出现在排口下游两三个网格，优先加密排口附近的 Δx。

修复操作与排查顺序严格对应：Pe_grid 超限就加密网格；显式超稳定条件就缩小 Δt；隐式步长过大就按峰值时间尺度重新选 Δt；排口局部振荡就在源项附近做局部网格加密。修改完成后，重新运行 5.1 的解析解对比，确认相对误差回到 5% 以内，再看浓度曲线是否保持单调下降。

本文还有配套的精品资源，点击获取