实验1-配线法求参

《地下水动力学》实验实验一：配线法求参一、实验名称配线法求参。二、实验性质必做。三、实验类型验证、研究。四、实验目的掌握配线法求参的原理，并利用EXCEL或者VisualBasic等编程工具进行配线法的编程分，根据抽水试验资料进行水文地质参数计算。五、实验内容5.1.原理配线法是《地下水动力学》中介绍的求取水文地质参数的一大类方法。5.1.1越流含水层完整井稳定流定流量抽水Hantush-Jacob公式：，为Bessel函数，可用Excel工程函数库计算，也可利用附录程序BessK0(x)计算。配线法公式：5.1.2无补给无限承压含水层完整井非稳定流定流量抽水Theis公式：，，Theis井函数计算方法很多，其中以多项式逼近最为常用。Allen(1954)andHastings(1955):式中系数如下：i0-0.577220.267773.9585010.999998.6347621.099652-0.2499118.0590225.6329630.055198.573339.573324-0.0097650.00108程序见附录Exp1(x)函数。配线法：5.1.3越流含水层完整井非稳定流定流量抽水Hantush-Jacob公式：，可利用附录程序w(x,y)函数计算。配线法：5.2配线法的Excel表本实验以Theis与Hantush-Jacob井函数为例，介绍用Excel表实现配线法的方法。5.2.1Theis井函数求参文件名：Theis_test.xlsTheis表中包含双对数坐标的配线法图形，A、B列为抽水实验观测数据(单位：min、m)。Q、R、S列分别为标准曲线的1/u、u、W(u)，图中的标准曲线即是用Q、S列生成的。滚动条ScrollBar1、ScrollBar2控制坐标偏移量。属性如下：Min=0，Max=100，LargeChange=10，SmallChange=1滚动条ScrollBar1、ScrollBar2的事件过程为：Submain()ScrollBar1_ChangeScrollBar2_ChangeEndSubPrivateSubScrollBar1_Change()WithScrollBar1delta=Range(H24).Value-Range(G24).ValueRange(I24).Value=Range(G24).Value+delta*.Value/100EndWithEndSubPrivateSubScrollBar1_Scroll()ScrollBar1_ChangeEndSubPrivateSubScrollBar2_Change()WithScrollBar2delta=Range(H26).Value-Range(G26).ValueRange(I26).Value=Range(G26).Value+delta*.Value/100EndWithEndSubPrivateSubScrollBar2_Scroll()ScrollBar2_ChangeEndSub横轴、纵轴的偏移量Min、Max分别为双对数坐标的偏移量，向左取负值。Offset为偏移量的对数值，单元格为I24与I26（记为）。单元L23抽水量（）,单元L24为观测孔距抽水井径向距离，L25为计算的T值（），L26为计算的值。$L$23为取单元L23的绝对引用。C、D列为A、B列偏移后的数值，E列为计算误差，公式如下：，参数计算公式：曲线：，曲线：，曲线：，Question表为抽水实验观测数据。Excel图形如下：5.2.2Hantush-Jacob井函数求参文件名：Hantush-Jacob_test.xlsData表含有不同值的Hantush-Jacob井函数数值表。Hantush-Jacob表中包含双对数坐标的配线法图形，A、B列为抽水实验观测数据，单位为min、m。标准曲线的绘图数据在Data表中。滚动条ScrollBar1、ScrollBar2控制坐标偏移量。属性如下：Min=0，Max=100，LargeChange=10，SmallChange=1滚动条ScrollBar1、ScrollBar2的事件过程为：Submain()ScrollBar1_ChangeScrollBar2_ChangeEndSubPrivateSubScrollBar1_Change()WithScrollBar1Delta=Range(H27).Value-Range(G27).ValueRange(I27).Value=Range(G27).Value+Delta*.Value/100EndWithEndSubPrivateSubScrollBar1_Scroll()ScrollBar1_ChangeEndSubPrivateSubScrollBar2_Change()WithScrollBar2Delta=Range(H29).Value-Range(G29).ValueRange(I29).Value=Range(G29).Value+Delta*.Value/100EndWithEndSubPrivateSubScrollBar2_Scroll()ScrollBar2_ChangeEndSub横轴、纵轴的偏移量Min、Max分别为双对数坐标的偏移量，向左取负值。Offset为偏移量的对数值，单元格为I27与I29（记为）。先移动图形，选择合适的，填入单元L25。单元L26抽水量（）,单元L27为观测孔距抽水井径向距离，L28为计算的T值（），L29为计算的值。$L$25为取单元L25的绝对引用。C、D列为A、B列偏移后的数值，E列为计算误差，公式如下：，参数计算公式：曲线：，曲线：，曲线：，Question表为抽水实验观测数据。Excel图形如下：附录：井函数计算程序'ThiscalculatesdrawdownsforflowtoawellfromtheTheissolution.FunctionExp1(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=1Ifx=1ThenExp1=-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)/xEndIfEndFunction'TocomputetheHantushleaky-aquiferfunctionW(x,y).Appendroutinesto'computeExp1(x),ExpInt(n,x),BessK0(x),BessI0(x)andBessI1(x).Functionw(x,y)Ifx=0Thenw=2*BessK0(y)Elser=1t=y^2/(4*x)b=2*xIfy=bThenw=0n=0Doterm=r*ExpInt(n+1,x)w=w+termn=n+1r=r*(-t)/nLoopUntilAbs(term)0.0000000001Elsew=2*BessK0(y)n=0Doterm=r*ExpInt(n+1,t)w=w-termn=n+1r=r*(-x)/nLoopUntilAbs(term)0.0000000001EndIfEndIfEndFunction'Tocomputetheexponentialintegralofordernforn=0.FunctionExpInt(n,x)Ifn=0ThenExpInt=Exp(-x)/xElseIfn=1ThenExpInt=Exp1(x)ElseIf(n1)And(x=5)Thena=Exp1(x)Fori=2Tona=(Exp(-x)-x*a)/(i-1)NextiExpInt=aElseIf(n1)And(x5)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)/tIfn=N1Theni=N1DoWhileini=i-1a=(Exp(-x)-i*a)/xLoopExpInt=aElsei=N1DoWhileini=i+1a=(Exp(-x)-x*a)/(i-1)LoopExpInt=aEndIfEndIfEndFunction'TocalculatethemodifiedBesselfunctionK0(x)for0xinfinity.FunctionBessK0(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.00053208Ifx=2Thent=(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)EndIfEndFunction'TocomputethemodifiedBesselfunctionI0(x)for0xinfinity.FunctionBessI0(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.00392377Ifx=3.75Thent=(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)EndIfEndFunction