当前位置：首页 > news >正文

别再死记硬背达西定律了！用Python模拟地下水流动，直观理解渗流速度与达西速度的区别

news 2026/6/1 8:23:30

用Python模拟地下水流动：可视化理解达西定律与渗流速度差异

从抽象公式到动态模型

当环境工程专业的学生第一次接触达西定律时，往往会被各种速度概念搞得晕头转向。达西速度、渗流速度、实际流速...这些术语在教材中通常以数学公式呈现，缺乏直观的物理图像。传统教学方法依赖静态示意图和理论推导，难以帮助学生建立清晰的概念框架。

我在教授地下水动力学课程时发现，即使是成绩优秀的学生，也常常混淆这些基本概念。直到我开始将Python编程引入教学，情况才发生根本改变。通过构建可交互的砂管水流模型，学生们能够亲眼看到水如何流过多孔介质，以及不同"速度"之间的本质区别。

1. 理论基础与模型构建

1.1 达西定律的物理意义

达西定律描述多孔介质中流体流动的基本规律，其核心表达式为：

Q = -K * A * (h2 - h1)/L

其中：

Q：体积流量（m³/s）
K：渗透系数（m/s）
A：横截面积（m²）
h1,h2：上下游水头（m）
L：砂段长度（m）

这个看似简单的公式背后隐藏着几个关键概念：

达西速度(q)：也称为比流量，计算公式为：

q = Q/A

它表示单位截面积上的流量，而非水的实际运动速度。

渗流速度(v)：水在孔隙中运动的表观速度，与有效孔隙度(θ)相关：

v = q/θ

1.2 有效孔隙度的关键作用

有效孔隙度θ是理解速度差异的核心参数。它表示参与流动的孔隙体积占总体积的比例：

参数	物理意义	典型值范围
总孔隙度	孔隙总体积/样品总体积	0.2-0.4
有效孔隙度	流动孔隙体积/样品总体积	0.1-0.3

注意：死端孔隙和不连通孔隙中的水不参与流动，因此有效孔隙度通常小于总孔隙度

2. Python实现砂管流动模型

2.1 环境配置与基础库

我们需要以下Python库构建模型：

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation

2.2 砂管模型的参数化表示

首先定义模型的基本参数：

# 砂管参数 L = 1.0 # 砂段长度 (m) A = 0.01 # 横截面积 (m²) K = 1e-4 # 渗透系数 (m/s) theta = 0.2 # 有效孔隙度 # 水头条件 h1 = 10.0 # 上游水头 (m) h2 = 9.0 # 下游水头 (m) # 计算基本量 Q = -K * A * (h2 - h1)/L # 体积流量 q = Q/A # 达西速度 v = q/theta # 渗流速度

2.3 流动可视化实现

创建动态展示渗流过程的动画：

def create_flow_animation(): fig, ax = plt.subplots(figsize=(10, 4)) # 初始化砂管和粒子 tube = plt.Rectangle((0, 0.4), L, 0.2, fc='tan', alpha=0.3) particles = [plt.Circle((0.1*i, 0.5), 0.01, fc='blue') for i in range(1, 10)] ax.add_patch(tube) for p in particles: ax.add_patch(p) ax.set_xlim(0, L) ax.set_ylim(0, 1) ax.set_aspect('equal') ax.set_title("砂管中的水流与粒子运动") def update(frame): for i, p in enumerate(particles): # 粒子位置更新 x, y = p.center new_x = x + v * 0.05 # 按渗流速度移动 if new_x > L: new_x = 0.1 # 循环边界 p.center = (new_x, y) return particles ani = FuncAnimation(fig, update, frames=100, interval=50, blit=True) plt.close() return ani

这段代码创建了一个直观的动画：

棕色半透明矩形代表砂管
蓝色圆点代表水分子或示踪剂粒子
粒子以渗流速度(v)运动，而非达西速度(q)

3. 溶质运移模拟

3.1 对流迁移的基本原理

溶质随水流迁移的过程可以用质量守恒描述：

∂(θC)/∂t = -∇·(qC) + 源汇项

其中C为溶质浓度。在稳定流条件下简化为：

∂C/∂t = -v·∇C

3.2 Python实现溶质锋面推进

模拟浓度锋面在砂管中的运动：

def simulate_solute_transport(): # 空间和时间离散 nx = 100 dx = L/nx nt = 200 dt = 0.005 # 初始条件 C = np.zeros(nx) C[:10] = 1.0 # 初始浓度锋面 # 存储结果 C_history = [C.copy()] for _ in range(nt): # 迎风格式求解 C[1:] = C[1:] - v * dt/dx * (C[1:] - C[:-1]) C_history.append(C.copy()) return np.array(C_history)

3.3 结果可视化

def plot_solute_transport(C_history): plt.figure(figsize=(10, 5)) for i in range(0, len(C_history), 20): plt.plot(np.linspace(0, L, 100), C_history[i], label=f't={i*0.005:.2f}s') plt.xlabel('距离 (m)') plt.ylabel('相对浓度') plt.title('溶质锋面随时间推进') plt.legend() plt.grid(True)

4. 参数影响分析

4.1 渗透系数的影响

保持其他参数不变，改变渗透系数K：

K (m/s)	达西速度q (m/s)	渗流速度v (m/s)
1e-5	1e-5	5e-5
1e-4	1e-4	5e-4
1e-3	1e-3	5e-3

可视化代码：

def plot_K_effect(): K_values = [1e-5, 1e-4, 1e-3] q_values = [-K * (h2 - h1)/L for K in K_values] v_values = [q/theta for q in q_values] plt.figure(figsize=(8, 4)) plt.plot(K_values, q_values, 'o-', label='达西速度 q') plt.plot(K_values, v_values, 's-', label='渗流速度 v') plt.xscale('log') plt.yscale('log') plt.xlabel('渗透系数 K (m/s)') plt.ylabel('速度 (m/s)') plt.legend() plt.grid(True)

4.2 有效孔隙度的影响

固定K=1e-4 m/s，变化θ：

theta_range = np.linspace(0.1, 0.3, 5) v_values = [q/theta for theta in theta_range] plt.figure(figsize=(8, 4)) plt.plot(theta_range, v_values, 'o-') plt.xlabel('有效孔隙度 θ') plt.ylabel('渗流速度 v (m/s)') plt.title('渗流速度随有效孔隙度变化') plt.grid(True)

5. 教学应用与扩展

5.1 交互式教学工具

使用Jupyter Widgets创建交互界面：

from ipywidgets import interact, FloatSlider def interactive_darcy(K=1e-4, theta=0.2, h1=10.0, h2=9.0): Q = -K * A * (h2 - h1)/L q = Q/A v = q/theta print(f"体积流量 Q = {Q:.2e} m³/s") print(f"达西速度 q = {q:.2e} m/s") print(f"渗流速度 v = {v:.2e} m/s") interact(interactive_darcy, K=FloatSlider(min=1e-5, max=1e-3, step=1e-5, value=1e-4), theta=FloatSlider(min=0.1, max=0.4, step=0.01, value=0.2), h1=FloatSlider(min=5, max=15, step=0.1, value=10), h2=FloatSlider(min=5, max=15, step=0.1, value=9))

5.2 三维扩展与复杂边界

将模型扩展到三维情况：

def darcy_3D(Kx, Ky, Kz, grad_h, theta): """ 三维达西定律实现 grad_h: 水头梯度向量 (dh/dx, dh/dy, dh/dz) """ q = -np.array([Kx, Ky, Kz]) * grad_h v = q/theta return q, v

6. 常见误区与验证方法

6.1 典型概念混淆

学生常见的理解误区包括：

认为达西速度就是水的实际流速
忽略有效孔隙度的影响
混淆质量流量与体积流量
不理解水头的能量意义

6.2 模型验证技巧

验证模拟结果的合理性：

质量守恒检查：流入量=流出量+储存变化
量纲一致性验证
极限情况测试（如K→0或θ→1）
与解析解对比（简单条件下）

def mass_balance_check(C_history, v, dx, dt): """检查溶质质量是否守恒""" total_mass = np.sum(C_history * dx, axis=1) rel_change = (total_mass - total_mass[0])/total_mass[0] plt.plot(rel_change) plt.title('相对质量变化') plt.xlabel('时间步') plt.ylabel('质量变化率')

7. 工程应用实例

7.1 污染羽运移预测

将模型应用于污染物迁移预测：

def contaminant_transport(): # 污染源参数 source_x = 0.2 # 污染源位置 C0 = 100.0 # 污染源浓度 (mg/L) # 模拟参数 nx = 200 x = np.linspace(0, L, nx) C = np.zeros(nx) C[(x >= source_x - 0.05) & (x <= source_x + 0.05)] = C0 # 模拟100天运移 for day in range(100): C[1:] = C[1:] - v * dt/dx * (C[1:] - C[:-1]) plt.plot(x, C) plt.title('100天后污染物分布') plt.xlabel('距离 (m)') plt.ylabel('浓度 (mg/L)')

7.2 地下水修复设计

基于模拟优化修复方案：

方案	抽水速率 (m³/d)	预计修复时间 (年)
1	10	5.2
2	20	3.8
3	30	2.5

8. 高级主题与扩展方向

8.1 非均质介质模拟

考虑渗透系数的空间变化：

K_field = np.ones(nx) * 1e-4 K_field[50:70] = 5e-5 # 低渗透带 def heterogeneous_flow(K_field): q = np.zeros_like(K_field) for i in range(1, nx): q[i] = -K_field[i] * (h2 - h1)/L return q

8.2 动态边界条件

模拟季节性水位变化：

h1_t = lambda t: 10 + 2*np.sin(2*np.pi*t/365) # 年周期变化 def time_dependent_simulation(days=365): C_history = [] C = np.zeros(nx) for day in range(days): h1 = h1_t(day) q = -K * (h2 - h1)/L v = q/theta # 更新浓度场 C[1:] = C[1:] - v * dt/dx * (C[1:] - C[:-1]) C_history.append(C.copy()) return np.array(C_history)

9. 性能优化技巧

9.1 向量化计算

替换循环为向量运算：

# 优化前 for i in range(1, nx): C[i] = C[i] - v * dt/dx * (C[i] - C[i-1]) # 优化后 C[1:] = C[1:] - v * dt/dx * (C[1:] - C[:-1])

9.2 使用Numba加速

from numba import jit @jit(nopython=True) def simulate_transport_numba(C, v, dt, dx, nt): for _ in range(nt): C[1:] = C[1:] - v * dt/dx * (C[1:] - C[:-1]) return C