实验1 配线法求参

单元 L26 抽水量（
）,单元 L27 为观测孔距抽水井径向距离，L28 为计算的 T 值（
），
L29 为计算的 值。
$L$25 为取单元 L25 的绝对引用。
C、D 列为 A、B 列偏移后的数值，E 列为计算误差，公式如下：
，
参数计算公式ຫໍສະໝຸດ Baidu 曲线：
，
曲线：
，
曲线：
， Question 表为抽水实验观测数据。 Excel 图形如下：
delta = Range("H24").Value - Range("G24").Value Range("I24").Value = Range("G24").Value + delta * .Value / 100 End With End Sub Private Sub ScrollBar1_Scroll() ScrollBar1_Change End Sub Private Sub ScrollBar2_Change() With ScrollBar2 delta = Range("H26").Value - Range("G26").Value Range("I26").Value = Range("G26").Value + delta * .Value / 100 End With End Sub Private Sub ScrollBar2_Scroll() ScrollBar2_Change
' To compute the exponential integral of order n for n >=0. Function ExpInt(n, x) If n = 0 Then
ExpInt = Exp(-x) / x ElseIf n = 1 Then
ExpInt = Exp1(x) ElseIf (n > 1) And (x <= 5) Then
End Sub
Private Sub ScrollBar1_Scroll() ScrollBar1_Change
End Sub
Private Sub ScrollBar2_Change() With ScrollBar2 Delta = Range("H29").Value - Range("G29").Value Range("I29").Value = Range("G29").Value + Delta * .Value / 100 End With
a = Exp1(x) For i = 2 To n
a = (Exp(-x) - x * a) / (i - 1) Next i ExpInt = a ElseIf (n > 1) And (x > 5) Then N1 = Int(x) t = x + N1 a = 1 + N1 / t ^ 2 + N1 * (N1 - 2 * x) / t ^ 4 + N1 * (6 * x ^ 2 - 8 * N1 * x + N1 ^ 2) / t ^ 6 a = a * Exp(-x) / t If n <= N1 Then
Min = 0，Max = 100，LargeChange = 10，SmallChange = 1
滚动条 ScrollBar1、ScrollBar2 的事件过程为： Sub main()
ScrollBar1_Change ScrollBar2_Change End Sub Private Sub ScrollBar1_Change() With ScrollBar1
' To compute the Hantush leaky-aquifer function W(x,y). Append routines to ' compute Exp1(x),ExpInt(n,x),BessK0(x),BessI0(x)and BessI1(x). Function w(x, y) If x = 0 Then
C、D 列为 A、B 列偏移后的数值，E 列为计算误差，公式如下：
，
参数计算公式： 曲线：
，
曲线：
，
曲线：
， Question 表为抽水实验观测数据。 Excel 图形如下：
5.2.2 Hantush-Jacob 井函数求参
文件名：Hantush-Jacob_test.xls
Data 表含有不同 值的 Hantush-Jacob 井函数数值表。
5.1.3 越流含水层完整井非稳定流定流量抽水 Hantush-Jacob 公式：
， 可利用附录程序 w(x, y)函数计算。
配线法：
5.2 配线法的 Excel 表
本实验以 Theis 与 Hantush-Jacob 井函数为例，介绍用 Excel 表实现配线法的方法。
5.2.1 Theis 井函数求参
《地下水动力学》实验
实验一：配线法求参
一、实验名称 配线法求参。
二、实验性质 必做。
三、实验类型 验证、研究。
四、实验目的 掌握配线法求参的原理，并利用 EXCEL 或者 Visual Basic 等编程工具进行配线法的编程分，根据抽水
试验资料进行水文地质参数计算。 五、实验内容
5.1. 原理 配线法是《地下水动力学》中介绍的求取水文地质参数的一大类方法。 5.1.1 越流含水层完整井稳定流定流量抽水 Hantush-Jacob 公式：
Exp1 = -Log(x) + A0 + x * (a1 + x * (a2 + x * (A3 + x * (A4 + x * A5)))) Else
p1 = B0 + x * (B1 + x * (B2 + x * (B3 + x * B4))) P2 = c0 + x * (c1 + x * (c2 + x * (c3 + x * C4))) Exp1 = (p1 / P2) * Exp(-x) / x End If End Function
， 计算。
配线法公式：
为 Bessel 函数，可用 Excel 工程函数库计算，也可利用附录程序 BessK0(x)
5.1.2 无补给无限承压含水层完整井非稳定流定流量抽水 Theis 公式：
，
，
Theis 井函数计算方法很多，其中以多项式逼近最为常用。 Allen (1954) and Hastings (1955):
文件名：Theis_test.xls
Theis 表中包含双对数坐标的配线法图形，A、B 列为抽水实验观测数据(单位：min、m)。Q、R、S 列 分别为标准曲线的 1/u、u、W(u)，图中的标准曲线即是用 Q、S 列生成的。
滚动条 ScrollBar1、ScrollBar2 控制坐标偏移量。属性如下：
i = N1 Do While i > n
i=i-1 a = (Exp(-x) - i * a) / x Loop ExpInt = a Else i = N1 Do While i < n i=i+1 a = (Exp(-x) - x * a) / (i - 1) Loop
ExpInt = a End If End If End Function
Hantush-Jacob 表中包含双对数坐标的配线法图形，A、B 列为抽水实验观测数据，单位为 min、m。标
准曲线的绘图数据在 Data 表中。
滚动条 ScrollBar1、ScrollBar2 控制坐标偏移量。属性如下：
Min = 0，Max = 100，LargeChange = 10，SmallChange = 1
End Sub
Private Sub ScrollBar2_Scroll() ScrollBar2_Change
End Sub
横轴、纵轴的偏移量 Min、Max 分别为双对数坐标的偏移量，向左取负值。Offset 为偏移量的对数值，
单元格为 I27 与 I29（记为 ）。
先移动图形，选择合适的 ，填入单元 L25。
滚动条 ScrollBar1、ScrollBar2 的事件过程为：
Sub main() ScrollBar1_Change ScrollBar2_Change
End Sub
Private Sub ScrollBar1_Change() With ScrollBar1 Delta = Range("H27").Value - Range("G27").Value Range("I27").Value = Range("G27").Value + Delta * .Value / 100 End With
End Sub 横轴、纵轴的偏移量 Min、Max 分别为双对数坐标的偏移量，向左取负值。Offset 为偏移量的对数值， 单元格为 I24 与 I26（记为 ）。
单元 L23 抽水量（
）,单元 L24 为观测孔距抽水井径向距离，L25 为计算的 T 值（
），
L26 为计算的 值。
$L$23 为取单元 L23 的绝对引用。
附录：井函数计算程序
' This calculates drawdowns for flow to a well from the Theis solution. Function Exp1(x) A0 = -0.57721566 a1 = 0.99999193 a2 = -0.24991055 A3 = 0.05519968 A4 = -0.00976004 A5 = 0.00107857 B0 = 0.2677737343 B1 = 8.6347608925 B2 = 18.059016973 B3 = 8.5733287401 B4 = 1 c0 = 3.9584969228 c1 = 21.0996530827 c2 = 25.6329561486 c3 = 9.5733223454 C4 = 1 If x <= 1 Then
w = 2 * BessK0(y) Else
r=1 t = y ^ 2 / (4 * x) b=2*x If y <= b Then
w=0 n=0 Do
term = r * ExpInt(n + 1, x) w = w + term n=n+1 r = r * (-t) / n Loop Until Abs(term) < 0.0000000001 Else w = 2 * BessK0(y) n=0 Do term = r * ExpInt(n + 1, t) w = w - term n=n+1 r = r * (-x) / n Loop Until Abs(term) < 0.0000000001 End If End If End Function
t = (x / 2) ^ 2 BessK0 = A0 + t * (a1 + t * (a2 + t * (A3 + t * (A4 + t * (A5 + t * A6))))) BessK0 = BessK0 - Application.Ln(x / 2) * BessI0(x) Else t=2/x BessK0 = B0 + t * (B1 + t * (B2 + t * (B3 + t * (B4 + t * (B5 + t * B6))))) BessK0 = BessK0 * Exp(-x) / Sqr(x) End If End Function
' To calculate the modified Bessel function K0(x) for 0<x<infinity. Function BessK0(x) A0 = -0.57721566 a1 = 0.4227842 a2 = 0.23069756 A3 = 0.0348859 A4 = 0.00262698 A5 = 0.0001075 A6 = 0.0000074 B0 = 1.25331414 B1 = -0.07832358 B2 = 0.02189568 B3 = -0.01062446 B4 = 0.00587872 B5 = -0.0025154 B6 = 0.00053208 If x <= 2 Then
式中系数如下：
i 0 1 2 3 4 5
程序见附录 Exp1(x)函数。
配线法：
-0.57722 0.99999 -0.24991 0.05519 -0.00976 0.00108
0.26777 8.63476 18.05902 8.57333
3.95850 21.09965 25.63296 9.57332