配线法的程序实现（Python）

欢迎点评！如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。

需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions

用 pip 安装：

pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions

ipywidgets 模块可实现交互绘图。

以下内容源自地下水动力学课程教学内容，计算井函数的程序要独立保存为 wellfunction.py，程序可在 Jupyter notebook 中运行。

配线法原理

Theis 配线法是 Hantush-Jacob 配线法的特例，在此只介绍后一种。

对

两边取对数

思路：

在标准曲线坐标系中移动 散点图，根据平移量计算 ：

记坐标平移量为

通过输入确定。

参数计算公式：

代码中给出一些关于初始参数的解释。

%matplotlib widget

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
import math as math
from wellfunction import *

# 控制小数的显示精度
np.set_printoptions(precision=4)

plt.rcParams['axes.unicode_minus']=False  #用来正常显示负号

# 准备数据, 根据配线法的种类，需要专门准备数据
Q = 528/1440  # m^3/min
r = 90 # m
t = np.array([1, 2, 4, 6, 9, 20, 30, 40, 50, 60, 90, 120, 150, 360, 550, 720]) # min
s = np.array([2.5, 3.9, 6.1, 8.0, 10.6, 16.8, 20.0, 22.6, 24.7, 26.4, 30.4, 33.0, 35.0, 42.6, 44.0, 44.5])/100 # m

ymin = 1.0e-2
ymax = 10.0

xmin = 1.0e-1
xmax = 1.0e4

x = np.linspace(-1, 4, 51)
x = np.float_power(10,x)

# 初始参数，最简单方法是取中间相邻两点用 Jacob 公式计算
# 取数组长度
n = len(t)
i1 = int(n/2)
i2 = i1 + 1
t1 = t[i1]
t2 = t[i2]
s1 = s[i1]
s2 = s[i2]
kk = (s1 - s2)/np.log10(t1/t2)

T = 0.183*Q/kk
S = 2.25*T*t1/r**2/np.float_power(10, s1/kk) 
  
B = 1e06

# Plot the data
def Hantush_Jacob_fit(dlogs, dlogt, beta):
    global T, S, B # 设置全局变量，值可以穿透程序
    
    plt.style.use('default')
    fig, ax = plt.subplots()
    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_aspect(1)
    plt.xlabel('$\log 1/u$')
    plt.ylabel('$W(u)$')
    ax.grid(True)
    
    T = Q/4/np.pi*np.float_power(10, dlogs)  # 配线法不同，公式也不一样
    S = 4*T/r**2*np.float_power(10, -dlogt)  # 配线法不同，公式也不一样
    u = r**2*S/4/T/x
    u0 = r**2*S/4/T/t
    ax.plot(x, hantush_jacob(u, beta),label="Hantush-Jacob曲线")
    ax.plot(x, theis(u),label="Theis曲线")
    ax.plot(t*np.float_power(10, dlogt), s*np.float_power(10, dlogs), '*', label="观测值")
    
    plt.legend(
        prop={'family': 'Simsun', 'size': 10}, handlelength=2,
        title="图例",
        title_fontproperties={'family': 'KaiTi', 'size': 12})
    ax.set_title("配线法", fontproperties={'family': 'KaiTi', 'size': 12}) # 指定显示中文字体

    plt.show()

    RSS = np.sum((hantush_jacob(u0, beta) - s)**2)

    print('             T(m^2/min): ', '{:.4f}'.format(T))
    print('                      S: ', '{:.4e}'.format(S))
    if beta>0 :
        print('                      B: ', '{:.4e}'.format(r/beta))
    print('Residual Sum of Squares: ', '{:.4e}'.format(RSS))

# ipywidgets 小部件控制参数, ipywidgets 有缓冲功能，同一个 Notebook 复制的代码得不到所期望的结果
dlogs_v=np.log10(T*4*np.pi/Q)
dlogt_v=np.log10(4*T/S/r/r)

widgets.interact(
    Hantush_Jacob_fit,
    dlogs = widgets.FloatSlider(
        value=dlogs_v, min=dlogs_v-1.0, max=dlogs_v+1.0, step=0.05,
        description=r"$\Delta\lg s$ [-]:", continuous_update=False,
        readout_format='.2f', disabled=False),
    dlogt = widgets.FloatSlider(
        value=dlogt_v, min=dlogt_v-1.0, max=dlogt_v+1.0, step=0.05,
        description=r'$\Delta\lg t$ [-]:', continuous_update=False,
        readout_format='.2f', disabled=False),
    beta = widgets.FloatSlider(
        value=r/B, min=0.0, max=+3.0, step=0.01,
        description=r'$\beta$ [-]:', continuous_update=False,
        readout_format='.2f', disabled=False),
    );   

image.png