《地下水动力学》实验报告班级：地质10-4班姓名：顾春生学号：05102064实验一：thise 曲线配线法求参一、实验名称配线法求参。
二、实验性质必做。
三、实验类型验证、研究。
四、实验目的掌握配线法求参的原理，并利用EXCEL 或者Visual Basic 等编程工具进行配线法的编程分，根据抽水试验资料进行水文地质参数计算。
五、实验内容5.1. 原理：：配线法是《地下水动力学》中介绍的求取水文地质参数的一大类方法。
5.1.2 无补给无限承压含水层完整井非稳定流定流量抽水Theis公式：，，配线法：参数计算公式：(曲线未作处理)曲线：，曲线：，Question表为抽水实验观测数据。
5.2.1 Theis井函数求参Theis公式：，，文件名：实验一（泰斯曲线.xls）Theis表中包含双对数坐标的配线法图形，A、B列为抽水实验观测数据(单位：min、m)。
Q、R、S列分别为标准曲线的1/u、u、W(u)，图中的标准曲线即是用Q、S列生成的。
横轴、纵轴的偏移量Min、Max分别为双对数坐标的偏移量，向左取负值。
Offset为偏移量的对数值，单元格为I24与I26（记为）。
单元L23抽水量（）,单元L24为观测孔距抽水井径向距离，L25为计算的T值（），L26为计算的值。
$L$23为取单元L23的绝对引用。
C、D列为A、B列偏移后的数值，E列为计算误差，公式如下：，在山西某地一承压含水层中进行抽水，含水层的顶、底板绝对隔水，抽水时侧向边界尚无影响。
14号孔为完整抽水井，抽水量为60m3/h。
2号观测孔距抽水孔43m，15号观测孔距抽水孔125m。
求导水系数和贮水系数。
(参考值：T=180m2/d，S=3.7×10-4)#2 s-t数据：#2 s-t曲线：分析：#15 s-t数据：#15 s-t曲线：分析：s-t /r*r数据：s-t /r*r曲线：分析：实验二：Hantush-Jacob井函数配线法求参一、实验名称配线法求参。
二、实验性质必做。
三、实验类型验证、研究。
四、实验目的掌握配线法求参的原理，并利用EXCEL或者Visual Basic等编程工具进行配线法的编程分，根据抽水试验资料进行水文地质参数计算。
五、实验内容2.1. 原理配线法是《地下水动力学》中介绍的求取水文地质参数的一大类方法。
2.2.1 越流含水层完整井稳定流定流量抽水Hantush-Jacob公式：，为Bessel函数，可用Excel工程函数库计算，也可利用附录程序BessK0(x)计算。
配线法公式：2.2.2 Hantush-Jacob井函数求参2.2.3 越流含水层完整井非稳定流定流量抽水Hantush-Jacob公式：可利用附录程序w(x, y)函数计算（成品见下页）。
，2.2.4题目：2.2.5配线需要掌握的过程：井函数计算程序(见末页附录)横轴、纵轴的偏移量Min、Max分别为双对数坐标的偏移量，向左取负值。
Offset为偏移量的对数值，单元格为I27与I29（记为）。
先移动图形，选择合适的，填入单元L25。
单元L26抽水量（）,单元L27为观测孔距抽水井径向距离，L28为计算的T值（），L29为计算的值。
$L$25为取单元L25的绝对引用。
C、D列为A、B列偏移后的数值，E列为计算误差，公式如下：，2.2.6配线结果：曲线：Hantush-Jacob表中包含双对数坐标的配线法图形，A、B列为抽水实验观测数据，单位为min、m。
标准曲线的绘图数据在Data 表中。
滚动条ScrollBar1、ScrollBar2控制坐标偏移量。
2.2.7结果分析：附录：井函数计算程序' This calculates drawdowns for flow to a well from the Theis solution.Function Exp1(x)A0 = -0.57721566a1 = 0.99999193a2 = -0.24991055A3 = 0.05519968A4 = -0.00976004A5 = 0.00107857B0 = 0.2677737343B1 = 8.6347608925B2 = 18.059016973B3 = 8.5733287401B4 = 1c0 = 3.9584969228c1 = 21.0996530827c2 = 25.6329561486c3 = 9.5733223454C4 = 1If x <= 1 ThenExp1 = -Log(x) + A0 + x * (a1 + x * (a2 + x * (A3 + x * (A4 + x * A5)))) Elsep1 = B0 + x * (B1 + x * (B2 + x * (B3 + x * B4)))P2 = c0 + x * (c1 + x * (c2 + x * (c3 + x * C4)))Exp1 = (p1 / P2) * Exp(-x) / xEnd IfEnd Function' 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 Thenw = 2 * BessK0(y)Elser = 1t = y ^ 2 / (4 * x)b = 2 * xIf y <= b Thenw = 0n = 0Doterm = r * ExpInt(n + 1, x)w = w + termn = n + 1r = r * (-t) / nLoop Until Abs(term) < 0.0000000001Elsew = 2 * BessK0(y)n = 0Doterm = r * ExpInt(n + 1, t)w = w - termn = n + 1r = r * (-x) / nLoop Until Abs(term) < 0.0000000001End IfEnd IfEnd Function' To compute the exponential integral of order n for n >=0.Function ExpInt(n, x)If n = 0 ThenExpInt = Exp(-x) / xElseIf n = 1 ThenExpInt = Exp1(x)ElseIf (n > 1) And (x <= 5) Thena = Exp1(x)For i = 2 To na = (Exp(-x) - x * a) / (i - 1)Next iExpInt = aElseIf (n > 1) And (x > 5) ThenN1 = Int(x)t = x + N1a = 1 + N1 / t ^ 2 + N1 * (N1 - 2 * x) / t ^ 4 + N1 * (6 * x ^ 2 - 8 * N1 * x + N1 ^ 2) / t ^ 6a = a * Exp(-x) / tIf n <= N1 Theni = N1Do While i > ni = i - 1a = (Exp(-x) - i * a) / xLoopExpInt = aElsei = N1Do While i < ni = i + 1a = (Exp(-x) - x * a) / (i - 1)LoopExpInt = aEnd IfEnd IfEnd Function' To calculate the modified Bessel function K0(x) for 0<x<infinity.Function BessK0(x)A0 = -0.57721566a1 = 0.4227842a2 = 0.23069756A3 = 0.0348859A4 = 0.00262698A5 = 0.0001075A6 = 0.0000074B0 = 1.25331414B1 = -0.07832358B2 = 0.02189568B3 = -0.01062446B4 = 0.00587872B5 = -0.0025154B6 = 0.00053208If x <= 2 Thent = (x / 2) ^ 2BessK0 = A0 + t * (a1 + t * (a2 + t * (A3 + t * (A4 + t * (A5 + t * A6)))))BessK0 = BessK0 - Application.Ln(x / 2) * BessI0(x)Elset = 2 / xBessK0 = B0 + t * (B1 + t * (B2 + t * (B3 + t * (B4 + t * (B5 + t * B6)))))BessK0 = BessK0 * Exp(-x) / Sqr(x)End IfEnd Function' To compute the modified Bessel function I0(x) for 0<x<infinity.Function BessI0(x)A0 = 1a1 = 3.5156229a2 = 3.0899424A3 = 1.2067492A4 = 0.2659732A5 = 0.0360768A6 = 0.0045813B0 = 0.39894228B1 = 0.01328592B2 = 0.00225319B3 = -0.00157565B4 = 0.00916281B5 = -0.02057706B6 = 0.02635537B7 = -0.01647633B8 = 0.00392377If x <= 3.75 Thent = (x / 3.75) ^ 2BessI0 = A0 + t * (a1 + t * (a2 + t * (A3 + t * (A4 + t * (A5 + t * A6)))))Elset = 3.75 / xBessI0 = B0 + t * (B1 + t * (B2 + t * (B3 + t * (B4 + t * (B5 + t * (B6 + t * (B7 + t * B8))))))) BessI0 = BessI0 * Exp(x) / Sqr(x)End IfEnd Function。