在环境科学中,污染物扩散是一个复杂的过程,涉及到多个因素,如水动力弥散系数、流速、浓度等。为了更直观地理解这个过程,我们可以使用数值方法进行模拟,并通过可视化展示结果。本文将介绍如何使用Python进行这种模拟,并使用matplotlib库进行可视化。
首先,我们需要导入一些必要的库,包括numpy和matplotlib。numpy用于进行数值计算,matplotlib用于绘制图形。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
然后,我们设置一些模拟参数,如空间大小、时间、水动力弥散系数、流速等。
h = 100 # 厚度为2的水层
L = 100 # 设置仿真空间大小
T = 20 # 设置仿真时间
D = 0.38 # 水动力弥散系数
v = 5.01 # 水的流速
k = 436382.0931027513 # 污染物扩散系数
接下来,我们创建空间网格和时间网格,并初始化浓度场。
dz = 5 # 离散化步长
dt = 2.5 # 离散化步长
z = np.arange(0, L, dz) # 空间网格
t = np.arange(0, T, dt) # 时间网格
C = np.zeros((len(z), len(t))) # 初始化浓度场
C[:, 0] = 0.495 # 初始浓度
然后,我们使用二阶中心差分法求解方程,并处理边界条件。
for j in range(1, len(t)):
for i in range(1, len(z) - 1):
C[i, j] = C[i, j - 1] + D * dt / dz ** 2 * (C[i + 1, j - 1] - 2 * C[i, j - 1] + C[i - 1, j - 1]) \
- v * dt / dz * (C[i + 1, j - 1] - C[i - 1, j - 1]) + k * C[i, j - 1] ** 2 * dt
if C[i, j] < 0:
C[i, j] = 0 # 防止出现负浓度
最后,我们使用曲面图可视化浓度场。
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(z, t)
surf = ax.plot_surface(X, Y, C.T, cmap='jet') # C数组需要转置
ax.set_xlabel('空间坐标z')
ax.set_ylabel('时间坐标t')
ax.set_zlabel('污染物浓度')
ax.set_title('污染物浓度随时间和空间的变化')
plt.show()
通过这种方式,我们可以直观地看到污染物浓度随时间和空间的变化情况。