河间地下水非稳定流演示（Python）
如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。
需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions
用 pip 安装：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
程序可在 Jupyter notebook 中运行。
mpl_interactions.ipyplot 模块可实现交互绘图，与 ipywidgets 相比，可变参数设定比较简单。
以下内容源自地下水动力学课程教学内容。
承压水
数学模型
式中，.
记 ，
，
，
。
河渠水位
式中
为河渠水位函数，
为河渠水位余函数，且
计算河渠水位的函数
import numpy as np
np.set_printoptions(precision=4)
# 河渠水位函数
def river_head(t, x, nmax):
if t <=0:
f = 0.0
else:
f = 1.0 - x
for n in range(1, nmax + 1):
alpha = n * np.pi
term = -2.0 * np.sin(alpha*x) * np.exp(-t*alpha**2) / alpha
f = f + term
return f
# 将函数向量化
river_head = np.vectorize(river_head)
水位函数与流量函数计算有限位数即可满足精度要求。
x = np.linspace(0,1,11)
t = 0.4
nmax = 20
print('水位 (t = {},'.format(t), 'nmax = {}):'.format(nmax))
print(river_head(t, x, nmax))
nmax = 50
print('水位 (t = {},'.format(t), 'nmax = {}):'.format(nmax))
print(river_head(t, x, nmax))
单宽流量
式中,
为河渠流量函数，
为河渠流量余函数，
.
计算河渠流量的函数
import numpy as np
# np.set_printoptions( ) 控制显示的小数位, 默认 8 位小数位. 可用 formatter 参数自定义显示的格式
# 例如: np.set_printoptions(formatter={'float': '{:.2f}'.format})
np.set_printoptions(precision=4)
# 河渠流量函数
def river_flow(t, x, nmax):
if t <=0:
g = 0.0
else:
g = 1.0
for n in range(1, nmax + 1):
alpha = n * np.pi
term = 2.0 * np.cos(alpha*x) * np.exp(-t*alpha**2)
g = g + term
return g
# 将函数向量化
river_flow = np.vectorize(river_flow)
x = np.linspace(0,1,11)
t = 0.4
nmax = 20
print('水位 (t = {},'.format(t), 'nmax = {}):'.format(nmax))
print(river_flow(t, x, nmax))
nmax = 50
print('水位 (t = {},'.format(t), 'nmax = {}):'.format(nmax))
print(river_flow(t, x, nmax))
潜水
数学模型
记 ，
，
河渠水位
单宽流量
式中
绘图演示 - 潜水
库文件导入与设置
%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
import mpl_interactions.ipyplot as iplt
设置一些初始参数
L = 1000 # aquifer length, m
K = 10 # hydraulic conductivity, m/d
H = 40 # aquifer top
zt = 30 # water level top, m
zb = 10 # aquifer bottom, m
mu = 0.15 # specific yield
a = K * 30 / mu
计算水头的函数
# hydraulic head
def head(x, h0_, hl_, h0t, hlt, t):
h02_ = h0_**2 - (h0_**2 - hl_**2) * x /L
t_ = a * t /L /L
x_ = x /L
dh0 = h0t*h0t-h0_*h0_
dhl = hlt*hlt-hl_*hl_
ht = h02_ + dh0*river_head(t_, x_, 20) + dhl*river_head(t_, 1 - x_, 20)
return np.sqrt(ht)
可交互的绘图程序
# make the interactive figure
x = np.linspace(0, L, 101)
t = np.linspace(0, 1000, 1001)
h0_ = np.linspace(zb, zt, 201)
hl_ = h0_
h0t = np.linspace(zb, H, 301)
hlt = h0t
fig, ax = plt.subplots(figsize=[6, 4])
ax.set_xlim(0, L)
ax.set_ylim(zb, H)
controls = iplt.plot(x, head, h0_=h0_, hl_=hl_,
h0t=h0t, hlt=hlt, t=t,
label="潜水位(m)")
plt.legend(prop={'family': 'KaiTi'},
handlelength=2, loc=2,
title="图例",
title_fontproperties={'family': 'SimHei'})
image.png
保存为矢量图
plt.savefig('plot.svg', format='svg')
欢迎点评！
