地下水动力学中Matlab的运用(井函数与贝塞尔函数)
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
地下水动力学中Matlab的运用
一、越流含水层中贝塞尔函数的实现
越流含水层中地下水向承压水井运动的问题中，贝塞尔函数大量运用，其中精确解中运用了零阶第二类虚宗量Bessel函数 ，一阶第二类虚宗量Bessel函数 。
经简化后的Hantush配线法中使用的是Hantush-Jacob公式，需要在双对数纸上绘制 曲线，而这在Matlab中很容易实现。Matlab中内置大量函数，其中包括五类Bessel函数，即besselj(nu,Z)、bessely(nu,Z)、besselh(nu,Z)、besseli(nu,Z)、besselk(nu,Z)，分别对应第一类贝塞尔函数，诺依曼函数，汉克尔函数，第一类修正贝塞尔函数以及第二类修正贝塞尔函数。而我们利用的即为第二类修正贝塞尔函数，相应的语句及图像如下：
再代入Theis公式中即可求得降深
S=3.1333m、5.4140m、7.7132m
由此可知同一观测点的水位降深随时间增加而增大，但增大的速率逐渐减小。
（2）将r=3m、30m、300m以及t=90min代入可得
u=0.00025、0.025、2.5
再代入前面的井函数程序中可得
W(u)=7.7171、3.1365、0.0249
解：由题知导水系数T=100m2/d，导压系数a=100m2/min，流量Q=1256 m3/d
Theis降深公式
其中
（1）将t=10min、100min、1000min以及r=10m代入可得
u=0.025、0.0025、0.00025
再代入前面的井函数程序中可得
W(u)=3.1365、5.4167、7.7171
fori=1:160
u(i)=10^(-15+i/10);%生成等比数列，便于画双对数坐标图像
end
fori=1:160
w(i)=mfun('Ei',1,u(i));
end
loglog(u,w);
gridon;
井函数数值的验证：
U
10-15
10-14
10-13
10-12
10-11
10-10
W(u)
33.9615607
此处最优解的原则是使目标函数值最小，也即目标函数决定了整体的优化评判标准。常见的方法是使理论值与实际值的方差最小，并认为此时函数曲线最契合实际数值点。所以Theis模型的目标函数如下
其中 表示权数，代表数据的可信度，一般设为1。 与 分别表示某时刻某一点水位值的理论值与实际值。Matlab语句编写如下：
6.331539
4.037930
1.822924
0.219384
0.000004
三、利用非线性规化获得Theis解析模型数值解
Theis公式既可以用于水位预测，也可以用于求含水层参数。而在解决后者这一逆问题时，传统的方法是配线法或Jacob直线图解法，其精度不高且随意性较大。利用Matlab中非线性规化函数fmincon的功能实现统一的评判标准，保证解的唯一性与可靠性。由于我们在上面已经解决了井函数的求值问题，故只需编写非线性规划的目标函数即可。
再代入Theis公式中即可求得降深
S=7.7132m、3.1349m、0.0249m
由此可知同一时刻，观测点与抽水井的距离越大，其水位降深越小。
习题2.在承压含水层中有一完整井，以抽水量Q=0.0058 m3/s进行抽水试验，在距抽水井10 m处有一观测孔，其观测资料如下表所示。试用配线法求该承压含水层的导水系数T和贮水系数μ*。
x=0:0.01:10;
y0=besselk(0,x);
y1=besselk(1,x);
loglog(x,y0,x,y1);
gridon;
二、井函数的实现
地下水向完整井的非稳定运动中需要运用井函数 ，其指数积分式为
在Matlab中利用quad或quad8等积命令可实现求其近似值。但Matlab中内置的Maple函数库中包含Ei函数，但不可直接显示其函数值，可直接利用mfun函数调用Matlab中的Maple函数库，以达到求值的要求。相应的语句及图像如下：
functionE =fune(x)
t=T;%调用时间点的M文件
s=S;%调用降深值的M文件
r=R;%调用观测距离的M文件
Q=q;%调用流量的M文件
E=0;
forn=1:length(r)
form=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;
累加时间/min
水位降深/m
累加时间/min
水位降深/m
1.5
0.130
15
0.78
2
0.171
20
0.90
2.5
0.218
30
1.00
3
0.30
50
1.20
4
0.36
70
1.32
6
0.50
110
1.59
8
0.58
200
1.71
10
0.62
解：将观测时间t=[1.5 2 2.5 3 4 6 8 10 15 20 30 50 70 110 200]/(60*24)
四、验证与计算
利用上述函数可根据所测得的数据计算得含水层的导水系数T与贮水系数 ，而所需要提供的数据包括观测时间，观测点与井的距离以及观测点的降深。由于编写的fune函数中可直接调用上述数据，故使用时只需将对应的数据矩阵导入文件夹即可。现利用课本上的习题进行验证结果的可靠性。
（1）利用《地下动力学》课本中100页中观2与观15的两个孔的数据进行测试，结果如下：
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)
end
end
目标函数实现后即可在非线性规划函数中调用，fmincon中变量很多，但一般Theis模型没有对应的线性不等约束或者线性等式约束，故只需限定上下限即可，而初值可根据经验取得相应的数量级即可。具体实现代码如下：
function[xfval]=funm
[x,fval] = fmincon('fune',[300;0.0001],[],[],[],[],[1;0.00001],[10000;0.01]);
通过上述验证我们可以确定Matlab的数值方法可靠度高，且计算过程简便。现计算所提供的两道习题：
习题1.在某均质、各向同性的承压含水层中，有一完整抽水井，其抽水量为1256 m3/d，已知含水层的导水系数为100 m2/d，导压系数为100 m2/min。试求：（1）抽水后10min、100min、1000min时，距抽水井10m处的水位降深，以及所反映水位降深的分布规律；（2）抽水后90min后距水井3m、30m、300m处的水位降深，以及所反映水位降深的分布规律。
目标函数值E=0.4092
书上利用配线法的计算结果如下：
目标函数值E=0.7426
（2）利用所提供的PDF文件中观测数据，上述Matlab程序计算结果如下：
目标函数值E=0.0394
PDF文件中自身程序计算结果为：
目标函数值E=0.1255
两者的区别源于井函数的求值部分的差异，本报告使用的是Matlab内置的Maple函数库，而PDF中利用的quad8积分函数。但两者计算结果相差不大，且理论上讲利用内置函数库的方法更加精确。
一、越流含水层中贝塞尔函数的实现
越流含水层中地下水向承压水井运动的问题中，贝塞尔函数大量运用，其中精确解中运用了零阶第二类虚宗量Bessel函数 ，一阶第二类虚宗量Bessel函数 。
经简化后的Hantush配线法中使用的是Hantush-Jacob公式，需要在双对数纸上绘制 曲线，而这在Matlab中很容易实现。Matlab中内置大量函数，其中包括五类Bessel函数，即besselj(nu,Z)、bessely(nu,Z)、besselh(nu,Z)、besseli(nu,Z)、besselk(nu,Z)，分别对应第一类贝塞尔函数，诺依曼函数，汉克尔函数，第一类修正贝塞尔函数以及第二类修正贝塞尔函数。而我们利用的即为第二类修正贝塞尔函数，相应的语句及图像如下：
再代入Theis公式中即可求得降深
S=3.1333m、5.4140m、7.7132m
由此可知同一观测点的水位降深随时间增加而增大，但增大的速率逐渐减小。
（2）将r=3m、30m、300m以及t=90min代入可得
u=0.00025、0.025、2.5
再代入前面的井函数程序中可得
W(u)=7.7171、3.1365、0.0249
解：由题知导水系数T=100m2/d，导压系数a=100m2/min，流量Q=1256 m3/d
Theis降深公式
其中
（1）将t=10min、100min、1000min以及r=10m代入可得
u=0.025、0.0025、0.00025
再代入前面的井函数程序中可得
W(u)=3.1365、5.4167、7.7171
fori=1:160
u(i)=10^(-15+i/10);%生成等比数列，便于画双对数坐标图像
end
fori=1:160
w(i)=mfun('Ei',1,u(i));
end
loglog(u,w);
gridon;
井函数数值的验证：
U
10-15
10-14
10-13
10-12
10-11
10-10
W(u)
33.9615607
此处最优解的原则是使目标函数值最小，也即目标函数决定了整体的优化评判标准。常见的方法是使理论值与实际值的方差最小，并认为此时函数曲线最契合实际数值点。所以Theis模型的目标函数如下
其中 表示权数，代表数据的可信度，一般设为1。 与 分别表示某时刻某一点水位值的理论值与实际值。Matlab语句编写如下：
6.331539
4.037930
1.822924
0.219384
0.000004
三、利用非线性规化获得Theis解析模型数值解
Theis公式既可以用于水位预测，也可以用于求含水层参数。而在解决后者这一逆问题时，传统的方法是配线法或Jacob直线图解法，其精度不高且随意性较大。利用Matlab中非线性规化函数fmincon的功能实现统一的评判标准，保证解的唯一性与可靠性。由于我们在上面已经解决了井函数的求值问题，故只需编写非线性规划的目标函数即可。
再代入Theis公式中即可求得降深
S=7.7132m、3.1349m、0.0249m
由此可知同一时刻，观测点与抽水井的距离越大，其水位降深越小。
习题2.在承压含水层中有一完整井，以抽水量Q=0.0058 m3/s进行抽水试验，在距抽水井10 m处有一观测孔，其观测资料如下表所示。试用配线法求该承压含水层的导水系数T和贮水系数μ*。
x=0:0.01:10;
y0=besselk(0,x);
y1=besselk(1,x);
loglog(x,y0,x,y1);
gridon;
二、井函数的实现
地下水向完整井的非稳定运动中需要运用井函数 ，其指数积分式为
在Matlab中利用quad或quad8等积命令可实现求其近似值。但Matlab中内置的Maple函数库中包含Ei函数，但不可直接显示其函数值，可直接利用mfun函数调用Matlab中的Maple函数库，以达到求值的要求。相应的语句及图像如下：
functionE =fune(x)
t=T;%调用时间点的M文件
s=S;%调用降深值的M文件
r=R;%调用观测距离的M文件
Q=q;%调用流量的M文件
E=0;
forn=1:length(r)
form=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;
累加时间/min
水位降深/m
累加时间/min
水位降深/m
1.5
0.130
15
0.78
2
0.171
20
0.90
2.5
0.218
30
1.00
3
0.30
50
1.20
4
0.36
70
1.32
6
0.50
110
1.59
8
0.58
200
1.71
10
0.62
解：将观测时间t=[1.5 2 2.5 3 4 6 8 10 15 20 30 50 70 110 200]/(60*24)
四、验证与计算
利用上述函数可根据所测得的数据计算得含水层的导水系数T与贮水系数 ，而所需要提供的数据包括观测时间，观测点与井的距离以及观测点的降深。由于编写的fune函数中可直接调用上述数据，故使用时只需将对应的数据矩阵导入文件夹即可。现利用课本上的习题进行验证结果的可靠性。
（1）利用《地下动力学》课本中100页中观2与观15的两个孔的数据进行测试，结果如下：
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)
end
end
目标函数实现后即可在非线性规划函数中调用，fmincon中变量很多，但一般Theis模型没有对应的线性不等约束或者线性等式约束，故只需限定上下限即可，而初值可根据经验取得相应的数量级即可。具体实现代码如下：
function[xfval]=funm
[x,fval] = fmincon('fune',[300;0.0001],[],[],[],[],[1;0.00001],[10000;0.01]);
通过上述验证我们可以确定Matlab的数值方法可靠度高，且计算过程简便。现计算所提供的两道习题：
习题1.在某均质、各向同性的承压含水层中，有一完整抽水井，其抽水量为1256 m3/d，已知含水层的导水系数为100 m2/d，导压系数为100 m2/min。试求：（1）抽水后10min、100min、1000min时，距抽水井10m处的水位降深，以及所反映水位降深的分布规律；（2）抽水后90min后距水井3m、30m、300m处的水位降深，以及所反映水位降深的分布规律。
目标函数值E=0.4092
书上利用配线法的计算结果如下：
目标函数值E=0.7426
（2）利用所提供的PDF文件中观测数据，上述Matlab程序计算结果如下：
目标函数值E=0.0394
PDF文件中自身程序计算结果为：
目标函数值E=0.1255
两者的区别源于井函数的求值部分的差异，本报告使用的是Matlab内置的Maple函数库，而PDF中利用的quad8积分函数。但两者计算结果相差不大，且理论上讲利用内置函数库的方法更加精确。