地下水動力学中Matlab的运用(井函数与贝塞尔函数).docxVIP

地下水动力学中Matlab的运用 越流含水层中贝塞尔函数的实现 越流含水层中地下水向承压水井运动的问题中，贝塞尔函数大量运用，其中精确解中运用了零阶第二类虚宗量Bessel函数K0(rB)，一阶第二类 s= 经简化后的Hantush-Jacob公式中也零阶第二类虚宗量Bessel函数K0 s≈ 在配线法中使用的是Hantush-Jacob公式，需要在双对数纸上绘制K0rB-rB曲线，而这在Matlab中很容易实现。Matlab中内置大量函数，其中包括五类Bessel函数，即besselj(nu,Z)、bessely(nu,Z)、 besselh(nu,Z)、 besseli(nu,Z)、 besselk( x=0:0.01:10; y0=besselk(0,x); y1=besselk(1,x); loglog(x,y0,x,y1); grid on; 井函数的实现 地下水向完整井的非稳定运动中需要运用井函数W(u)，其指数积分式为 W 在Matlab中利用quad或quad8等积命令可实现求其近似值。但Matlab中内置的Maple函数库中包含Ei函数，但不可直接显示其函数值，可直接利用mfun函数调用Matlab中的Maple函数库，以达到求值的要求。相应的语句及图像如下： for i=1:160 u(i)=10^(-15+i/10); %生成等比数列，便于画双对数坐标图像 end for i=1:160 w(i)=mfun(Ei,1,u(i)); end loglog(u,w); grid on; 井函数数值的验证： U 10-15 10-14 10-13 10-12 10-11 10-10 W(u) 33.9615607 31.658976 29.356391 27.053805 24.751220 22.448635 U 10-9 10-8 10-7 10-6 10-5 10-4 W(u) 20.146050 17.843465 15.540880 13.238296 10.935720 8.633225 U 10-3 10-2 10-1 100 101 W(u) 6.331539 4.037930 1.822924 0.219384 0.000004 利用非线性规化获得Theis解析模型数值解 Theis公式既可以用于水位预测，也可以用于求含水层参数。而在解决后者这一逆问题时，传统的方法是配线法或Jacob直线图解法，其精度不高且随意性较大。利用Matlab中非线性规化函数fmincon的功能实现统一的评判标准，保证解的唯一性与可靠性。由于我们在上面已经解决了井函数的求值问题，故只需编写非线性规划的目标函数即可。 此处最优解的原则是使目标函数值最小，也即目标函数决定了整体的优化评判标准。常见的方法是使理论值与实际值的方差最小，并认为此时函数曲线最契合实际数值点。所以Theis模型的目标函数如下 E 其中wij表示权数，代表数据的可信度，一般设为1。Hjti与 function E =fune(x) t=T; %调用时间点的M文件 s=S; %调用降深值的M文件 r=R; %调用观测距离的M文件 Q=q; %调用流量的M文件 E=0; for n=1:length(r) for m=1:length(t) a=r(n)^2*x(2)/(4*x(1)*t(m)); E=E+(Q/(4*pi*x(1))*mfun(Ei,1,a)-s(n,m))^2; end end 目标函数实现后即可在非线性规划函数中调用，fmincon中变量很多，但一般Theis模型没有对应的线性不等约束或者线性等式约束，故只需限定上下限即可，而初值可根据经验取得相应的数量级即可。具体实现代码如下： function [x fval]=funm [x,fval] = fmincon(fune,[300;0.0001],[],[],[],[],[1;0.00001],[10000;0.01]); 验证与计算 利用上述函数可根据所测得的数据计算得含水层的导水系数T与贮水系数μ*，而所需要提供的数据包括观测时间，观测点与井的距离以及观测点的降深。由于编写的fune （1）利用《地下动力学》课本中100页中观2与观15的两个孔的数据进行测试，结果如下： T=195.94m2/d μ* 书上利用配线法的计算结果如下： T=197.67m2/d μ （2）利用所提供的PDF文件中观测数据，上述Matlab程序计算结果如下： T=84.82m2/d μ*=1.46x PDF文件中自身程序计算结果为： T=77.37m2/d μ* 两者的区别源于井函数的求值部