4. 井函数W(u)的数值计算方法

Theis公式的核心是井函数 \(W(u)\)，其中 \(u=\dfrac{r^2 S}{4Tt}\)。降深 \(s=\dfrac{Q}{4\pi T}W(u)\)。实际计算中，\(W(u)\) 没有初等闭式表达，必须依靠数值方法。工程上常用的三条路径是：级数求和、渐近展开、查表/近似公式。理解它们的适用区间与精度差异，是避免“公式用错区间”这类系统性错误的关键。

一、级数求和：小 u 区间的精确解

井函数的标准定义为指数积分 \(W(u)=\displaystyle\int_u^{\infty}\frac{e^{-x}}{x}dx\)，其级数展开为：

W(u) = -0.5772156649 - ln(u) + u - u²/(2·2!) + u³/(3·3!) - u⁴/(4·4!) + …
     = -γ - ln(u) + Σ_{n=1}^{∞} (-1)^{n+1} u^n / (n·n!)

其中 γ≈0.5772156649 为欧拉常数。该级数对任意 u>0 收敛，但收敛速度随 u 增大迅速变慢：u 较小时，前几项即可达到很高精度；u 接近 1 时需十余项；u 大于 2 后出现明显的“大数相减”抵消误差，双精度下甚至可能失真。因此级数求和的实际适用区间是 u ≤ 0.1（工程上常取 u ≤ 0.05 更稳妥），此时取前 8～10 项即可保证 1e-8 量级精度。

需要注意：抽水初期（t 小）或近井（r 小）时 u 反而较大，级数法恰恰在这些时刻不适用；抽水时间较长、u 变小时，级数法才成为主力。

二、渐近展开：大 u 区间的有效工具

当 u 较大时，指数积分可用渐近展开：

W(u) ≈ e^{-u}/u · [1 - 1!/u + 2!/u² - 3!/u³ + 4!/u⁴ - …]

这是发散的渐近级数，并非收敛级数，其意义是“截断到最小项”时误差最小。实践中的判据是：当 u > 10 时，仅取首项 \(W(u)\approx e^{-u}/u\) 的相对误差已小于约 10%；u > 20 时误差约 5% 以内；若取到 \(1-1/u+2/u^2\) 三项，u > 10 时相对误差可降到 1% 以下。因此渐近展开的适用区间是 u ≥ 10，且项数不宜贪多——超过最小项后继续加项反而使误差增大。

物理上，大 u 对应抽水极早期或距井极远的位置，此时降深很小，渐近展开正好覆盖这一“信号微弱”的区域。

三、常用近似公式与精度对比

工程手册中流传最广的是 Jacob 近似（又称 Cooper-Jacob 近似）：当 u < 0.01 时，忽略级数中 u 及更高次项，得

W(u) ≈ -0.5772 - ln(u) = ln(2.25 T t / (r² S))

其误差来源就是被丢弃的 u 项：u=0.01 时相对误差约 0.5%，u=0.05 时约 2.5%，u=0.1 时约 5%，u=0.3 时可达 15% 以上。所以 Jacob 公式必须严格守住 u ≤ 0.01～0.05 的门槛，否则会系统性高估降深。

另一类是高精度有理逼近（如 Abramowitz & Stegun 的指数积分近似），在全区间内相对误差可控制在 1e-7 量级，适合编程实现，但公式较长、不便手算。查表法则介于两者之间：以 u 的对数为横坐标、W(u) 为纵坐标制成数表，配合对数插值，精度约 0.1%～1%，适合无计算器时代的现场估算，如今更多作为校核手段。

四、编程实现与查表方法

下面给出一个可直接使用的 Python 实现，按 u 自动切换算法，兼顾精度与稳健性：

import math

def well_function(u):
    if u <= 0:
        raise ValueError("u 必须为正")
    if u < 0.05:                      # 级数求和
        s = -0.5772156649015329 - math.log(u)
        term, n = u, 1
        while abs(term) > 1e-12:
            s += term
            n += 1
            term *= -u * (n - 1) / (n * n)
        return s
    if u >= 10:                       # 渐近展开，取到最小项
        s, term, n = 1.0, 1.0, 1
        while True:
            term *= -n / u
            if abs(term) < 1e-12 or abs(term) > abs(prev):
                break
            prev = term
            s += term
            n += 1
        return math.exp(-u) / u * s
    # 过渡带：用高精度有理逼近或数值积分兜底
    return exp1_rational(u)           # 相对误差 < 1e-7

实现要点有三：一是级数循环用递推项 \(t_n=t_{n-1}\cdot(-u)(n-1)/n^2\) 避免重复求阶乘；二是渐近展开必须在“项开始变大”时立即截断，否则发散级数会越加越错；三是过渡带不要用 Jacob 硬撑，应调用有理逼近或直接数值积分。

查表法的操作流程是：先算 \(u=r^2S/(4Tt)\)，取对数 \(\lg u\)，在 W(u) 表中定位相邻两行，按线性（或对数）插值得到 W(u)，再乘 \(Q/(4\pi T)\) 得降深。查表时务必确认表的 u 列是对数刻度还是线性刻度，插值方式不同会带来 0.5% 左右的差异。

最后强调一个易错点：无论用哪种方法，都要先统一单位再算 u。r 用 m、S 无量纲、T 用 m²/d、t 用 d，则 u 无量纲；若 T 用 m²/s 而 t 用 d，u 会差 86400 倍，直接导致方法选错区间。建议在代码入口处做一次单位归一化，并在输出中同时打印 u 值与所选算法，便于复核。