地下水动力学中井函数的计算
地下动力学中的定流量抽水的井流模型，可以写成统一的形式
井函数的计算方法散见有关文献。以下介绍 Theis、Hantush - Jacob、Neuman 三种函数的计算方法。
以下内容源自地下水动力学课程教学内容。
1. Theis井函数
这是一个指数积分，有多种计算方法，如教材中的级数求和法。在此介绍另外两种计算方法。
1.1 SciPy 科学计算库
Python 的 SciPy 库有 scipy.special.exp1 函数可以计算井函数。
%matplotlib inline
# 导入一些库
import matplotlib.pyplot as plt
import numpy as np
import scipy.special as sps
np.set_printoptions(precision=4)
# 定义一个测试函数
def timer(func):
def func_wrapper(*args, **kwargs):
from time import time
time_start = time()
result = func(*args, **kwargs)
time_end = time()
time_spend = time_end - time_start
print("%s cost time: %.6f s" % (func.__name__, time_spend))
return result
return func_wrapper
用 sps.exp1 计算井函数
@timer
def theis1(u):
return sps.exp1(u)
u = np.array([10**x for x in range(-10, 1)])
print(theis1(u))
theis1 cost time: 0.007292 s
[22.4486 20.1461 17.8435 15.5409 13.2383 10.9357 8.6332 6.3315 4.0379
1.8229 0.2194]
1.2 多项式逼近方法
Theis 井函数可以分段用多项式逼近。公式如下
-
时
式中
-
时
式中
[1]Abramowitz, Milton, Ed. "Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards Applied Mathematics Series 55. Tenth Printing. " Engineering (1972):1076.
@timer
def theis2(u):
"""
多项式逼近方法计算 Theis 井函数.
u = r^2 S/(4Tt), u 为数组时返回数组。
"""
def wellfunc(u):
a = [-0.57721566, 0.99999193, -0.24991055, 0.05519968, -0.00976004, 0.00107857]
b = [0.2677737343, 8.6347608925, 18.059016973, 8.5733287401]
c = [3.9584969228, 21.0996530827, 25.6329561486, 9.5733223454]
if u <= 1:
w = (
-np.log(u) + a[0]
+ u * (a[1] + u * (a[2] + u * (a[3] + u * (a[4] + u * a[5]))))
)
else:
w = c[0] + u * (c[1] + u * (c[2] + u * (c[3] + u)))
w = (b[0] + u * (b[1] + u * (b[2] + u * (b[3] + u)))) / w
w = w * np.exp(-u) / u
return w
well = np.vectorize(wellfunc) # 向量化函数
return 1.0 * well(u)
u = np.array([10**x for x in range(-10, 1)])
print(theis2(u))
theis2 cost time: 0.000000 s
[22.4486 20.1461 17.8435 15.5409 13.2383 10.9357 8.6332 6.3315 4.0379
1.8229 0.2194]
2. Hantush-Jacob 井函数
式中：.
Hantush-Jacob 井函数还有如下的性质
该性质也为求参的拐点法所用。
Hunt（1977）给出了 Hantush-Jacob 井函数的级数形式：
级数形式（Hunt，1977）
式中，
为指数积分，当 时级数快速收敛。
[1] Hunt, Bruce . "Calculation of the leaky aquifer function." Journal of Hydrology 33.1-2(1977):179-183.
科学计算库 scipy.special.expn 可用于计算该指数积分，scipy.special.k0 是计算 0 阶第二类修正 Bessel 函数 的程序。
下面的代码是 Hunt 方法的实现。
@timer
def hantush_jacob(u, beta):
"""
指数积分级数方法计算 Hantush-Jacob 井函数.
u = r^2 S/(4Tt);
beta = r/B;
u, beta 可以为数组并返回数组。
"""
def wellfunc(u, beta):
if u < 0:
print("Negative are not allowed")
return np.nan
if u == 0:
return 2.0 * sps.k0(beta)
r = 1
t = beta**2 / (4 * u)
b = 2 * u
if beta <= b: # beta<2u
W = 0
n = 0
term = r * sps.expn(n + 1, u)
while np.abs(term) > 1e-10:
W = W + term
n = n + 1
r = r * (-t) / n
term = r * sps.expn(n + 1, u)
else:
W = 2.0 * sps.k0(beta)
n = 0
term = r * sps.expn(n + 1, t)
while np.abs(term) > 1e-10:
W = W - term
n = n + 1
r = r * (-u) / n
term = r * sps.expn(n + 1, t)
return W
well = np.vectorize(wellfunc)
return 1.0 * well(u, beta)
beta = 0.05
u = np.array([10**x for x in range(-4, 1)])
print(hantush_jacob(u, beta))
hantush_jacob cost time: 0.001002 s
[6.2282 5.7965 3.9795 1.8184 0.2193]
beta = np.array([0.01 * x for x in range(1, 6)])
print(beta)
u = np.array([10**x for x in range(-4, 1)])
print(hantush_jacob(u, beta))
[0.01 0.02 0.03 0.04 0.05]
hantush_jacob cost time: 0.000000 s
[8.3983 6.2347 4.0167 1.82 0.2193]
3. Neuman 井函数
Neuman 井函数计算比较复杂，Cheng, A.H.-D. 的著作中有详细介绍
- Cheng, A.H.-D., Multilayered Aquifer Systems-Fundamentals and Applications, Marcel Dekker, New York/Basel, 384 p., 2000.
对于较短的时间，有
对于较长的时间，有
式中
由于 是振荡衰减的，先将积分转换为
零点区间积分的级数和，每项用 Gaussian 正交多项式插值计算，由于级数和是衰减收敛的，计算该级数和用到了 epsilon 算法。
from math import *
import numpy as np
import scipy.special as sp
import time
# Modified Bessel function K_0
def besselk0(x):
if x <= 2.0:
y = x*x/4.0
besselk0 = (-log(x/2.)*besseli0(x) - 0.57721566 \
+ y*(0.4227842 + y*(0.23069756 \
+ y*(0.348859e-1 + y*(0.262698e-2 \
+ y*(0.1075e-3 + y*0.74e-5))))))
else:
y = (2.0/x)
besselk0 = (exp(-x)/sqrt(x)*(1.25331414 \
+ y*(-0.7832358e-1 + y*(0.2189568e-1 \
+ y*(-0.1062446e-1 + y*(0.587872e-2 \
+ y*(-0.25154e-2 + y*0.53208e-3)))))))
return besselk0
# Modified Bessel function K_1
def besselk1(x):
if x <= 2.0:
y = x*x/4.0
besselk1 = (log(x/2.)*besseli1(x) \
+ (1./x)*(1. + y*(0.15443144 + y*(-0.67278579 \
+ y*(-0.18156897 + y*(-0.1919402e-1 \
+ y*(-0.110404e-2 + y*(-0.4686e-4))))))))
else:
y = 2.0/x
besselk1 = (exp(-x)/sqrt(x)*(1.25331414 \
+ y*(0.23498619 + y*(-0.3655620e-1 \
+ y*(0.1504268e-1 + y*(-0.780353e-2 \
+ y*(0.325614e-2 + y*(-0.68245e-3))))))))
return besselk1
# Modified Bessel function I_0
def besseli0(x):
if x < 3.75:
y = (x/3.75)**2
besseli0 = (1.+ y*(3.5156229 + y*(3.0899424 \
+ y*(1.2067492 + y*(0.2659732 \
+ y*(0.360768e-1 + y*0.45813e-2))))))
else:
y = 3.75/x
besseli0 = (exp(x)/sqrt(x)*(0.39894228 \
+ y*(0.1328592e-1 + y*(0.225319e-2 \
+ y*(-0.157565e-2 + y*(0.916281e-2 \
+ y*(-0.2057706e-1 + y*(0.2635537e-1 \
+ y*(-0.1647633e-1 + y*0.392377e-2)))))))))
return besseli0
# Modified Bessel function I_1
def besseli1(x):
if x < 3.75:
y = (x/3.75)**2
besseli1 = (x*(0.5e0 + y*(0.87890594 \
+ y*(0.51498869 + y*(0.15084934 \
+ y*(0.2658733e-1 + y*(0.301532e-2 \
+ y*0.32411e-3)))))))
else:
y = 3.75/x
besseli1 = (exp(x)/sqrt(x)*(0.39894228 \
+ y*(-0.3988024e-1 + y*(-0.362018e-2 \
+ y*(0.163801e-2 + y*(-0.1031555e-1 \
+ y*(0.2282967e-1 + y*(-0.2895312e-1 \
+ y*(0.1787654e-1 + y*(-0.420059e-2))))))))))
return besseli1
# Bessel function J_0
def besselj0(x):
if abs(x) < 8.:
y = x**2
besselj0 = (
(57568490574.+y*(-13362590354.+y*(651619640.7 +
y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))))) /
(57568490411.+y*(1029532985.+y*(9494680.718+y *
(59272.64853+y*(267.8532712+y)))))
)
else:
ax = abs(x)
z = 8.0/ax
y = z**2
xx = ax-.785398164
besselj0 = (
sqrt(.636619772/ax)*(cos(xx)*(1.+y *
(-.1098628627e-2+y*(.2734510407e-4+y * (-.2073370639e-5 +
y*.2093887211e-6))))-z*sin(xx)*(-.1562499995e-1+y *
(.1430488765e-3+y*(-.6911147651e-5+y*(.7621095161e-6+y
* (-.934945152e-7)))))) )
return besselj0
# Bessel function J_1
def besselj1(x):
if abs(x) < 8.:
y = x**2
besselj1 = (x*(72362614232. \
+ y*(-7895059235. + y*(242396853.1 \
+ y*(-2972611.439 + y*(15704.48260 \
+ y*(-30.16036606))))))/(144725228442. \
+ y*(2300535178. + y*(18583304.74 \
+ y*(99447.43394+y*(376.9991397+y))))))
else:
ax = abs(x)
z = 8.0/ax
y = z**2
xx = ax-2.356194491
besselj1 = (sqrt(.636619772/ax)*(cos(xx)*(1. \
+ y*(.183105e-2 + y*(-.3516396496e-4 \
+ y*(.2457520174e-5 + y*(-.240337019e-6))))) \
- z*sin(xx)*(.04687499995 + y*(-.2002690873e-3 \
+ y*(.8449199096e-5 + y*(-.88228987e-6 \
+ y*.105787412e-6)))))*np.sign(complex(1.0, x)).real)
return besselj1
# Find a root of BesselK0 near rstart by Newton-Raphson method
def bj0root(rstart):
xacc = 1.0e-8
bj0root = rstart
for j in range(100):
dx = -besselj0(bj0root) / besselj1(bj0root)
bj0root = bj0root - dx
if abs(dx) < xacc:
break
return bj0root
# Calculate Gaussian quadrature nodes and weights
# Source: Numerical Recipes, Press, et al. 1992
def gauleg(x1, x2, n):
x, w = np.polynomial.legendre.leggauss(n)
xm = 0.5*(x2 + x1)
xl = 0.5*(x2 - x1)
return xm + xl*x[:], xl*w[:]
# Extrapolate a series by epsilon algorithm
def epsilonn(n, psum):
eps = np.zeros((202, 201))
for m in range(n + 1):
eps[0, m] = psum[m]
for i in range(1, n + 1):
for m in range(n - i, -1, -1):
eps[i, m] = eps[i - 2, m + 1] + 1.0 / \
(eps[i - 1, m + 1] - eps[i - 1, m])
return eps[n, 0]
def func(u, gamma, x, ncase):
if ncase == 1:
func = 64.0*x*besselj0(x*np.sqrt(gamma))*summ(x, u, gamma, ncase)
if ncase == 2:
func = (4.0*x*besselj0(x*np.sqrt(gamma))*((1.0 \
- np.exp(-gamma*x*np.tanh(x)/(4.0*u)))*np.tanh(x)/(2.0*x**3) \
+ summ(x, u, gamma, ncase)))
return func
def summ(x, u, gamma, ncase):
nsum = 12
psum = np.zeros(201)
for m in range(1, nsum + 1):
a = 4.0*x**2 + (2.0*m - 1.0)**2*np.pi**2
denom = (2.0*m - 1.0)**2*np.pi**2*a
if ncase == 1:
dsum = (1 - np.exp(-gamma/16.0/u*a))/denom
if ncase == 2:
dsum = 16.0/denom
psum[m] = psum[m - 1] + dsum
# return epsilonn(nsum, psum)
return psum[nsum]
# Neuman's type A and B unconfined aquifer well function
def neumanw(u, gamma, ncase):
r = np.zeros(201)
# Define integration ranges using even roots of J0
nr = 200
r[0] = 0.0
r[1] = 5.520078110080565
for i in range(2, nr+1):
r[i] = bj0root(r[i - 1] + 2*np.pi)
ng = 30
x = np.zeros(30)
w = np.zeros(30)
wsum = np.zeros(201)
gammasq = np.sqrt(gamma)
wsum[0] = 0.0
# Integrate by subintervals
for j in range(1, nr+1):
x, w = gauleg(r[j - 1] / gammasq, r[j] / gammasq, ng)
dwf = 0.0
# Perform Gaussian quadrature
for i in range(ng):
dwf = dwf + w[i] * func(u, gamma, x[i], ncase)
wsum[j] = wsum[j - 1] + dwf
# Use epsilon algorithm to extrapolate
return epsilonn(nr, wsum)
定义 Neuman 井函数
@timer
def Neumann(ua, uy, gamma):
n = len(gamma)
m = len(ua)
value = np.zeros((m, n))
for i in range(m):
for j in range(n):
value[i, j] = neumanw(ua[i], gamma[j], 1)
np.savetxt("w_ua1.csv", value, delimiter=",")
m = len(uy)
value = np.zeros((m, n))
for i in range(m):
for j in range(n):
value[i, j] = neumanw(uy[i], gamma[j], 2)
np.savetxt("w_uy1.csv", value, delimiter=",")
return "Data has been written to file: w_ua.csv and w_uy.csv."
根据 ua, uy, gamma 计算井函数，并将结果保存到 csv 文件。下面程序中 的取值，转换为
后，计算结果可以与教材中的表直接对比。
ua = np.array([
2.50E+00, 1.25E+00, 7.14E-01, 4.17E-01, 2.50E-01,
1.25E-01, 7.14E-02, 4.17E-02, 2.50E-02, 1.25E-02,
7.14E-03, 4.17E-03, 2.50E-03, 1.25E-03, 7.14E-04,
4.17E-04, 2.50E-04, 1.25E-04, 7.14E-05, 4.17E-05])
uy = np.array([
2.50E+03, 1.25E+03, 7.14E+02, 4.17E+02, 2.50E+02,
1.25E+02, 7.14E+01, 4.17E+01, 2.50E+01, 1.25E+01,
7.14E+00, 4.17E+00, 2.50E+00, 1.25E+00, 7.14E-01,
4.17E-01, 2.50E-01, 1.25E-01, 7.14E-02, 4.17E-02,
2.50E-02, 1.25E-02, 7.14E-03, 4.17E-03, 2.50E-03])
gamma = np.array([0.001, 0.004, 0.010, 0.030, 0.060,
0.100, 0.200, 0.400, 0.600, 0.800,
1.000, 1.500, 2.000, 2.500, 3.000,
4.000, 5.000, 6.000, 7.000])
Neumann(ua, uy, gamma)
从计算过程可以看出，计算一张井函数表用时还是比较长的（429.6 s），可以考虑用 numba 加速，用笔者的计算机配置，速度可以提高近 5 倍。
上述代码没有用 scipy.special 中的 Bessel 函数程序就是为了使程序可用 numba 加速。
