水侵量
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
附录：
图1 无因次水侵曲线（有限水区）
源程序
（1）Laplace变换函数
% 无因次水侵量的Laplace变换计算公式
function Qs=Laplace(s,reD)
Qs=-(besselk(1,sqrt(s)*reD).*besseli(1,sqrt(s))-besseli(1,sqrt(s)*reD ).*besselk(1,sqrt(s)))./(s.^(3/2).*(besselk(1,sqrt(s)*reD).*besseli(0 ,sqrt(s))+besselk(0,sqrt(s)).*besseli(1,sqrt(s)*reD)));
Return
（2）主程序
%主程序
clear
N=12;
for i=1:1:N
Vi=0;
n1=fix((i+1)/2);
if i>N/2
n2=N/2;
else
n2=i;
end
for k=n1:1:n2
Vi=Vi+(-1)^(N/2+i)*k^(N/2)*factorial(2*k)/(factorial(N/2-k)*factorial (k)*factorial(k-1)*factorial(i-k)*factorial(2*k-i));
end
V(i)=Vi;
end
%绘图
for reD=1.5:0.5:5
QtD=0.0;
switch reD
case 1.5
tD=0.1:0.05:0.5;
case 2
tD=0.1:0.05:3;
case 2.5
tD=0.1:0.05:9;
case 3
tD=0.1:0.05:15;
case 3.5
tD=0.1:0.05:30;
case 4
tD=0.1:0.05:45;
case 4.5
tD=0.1:0.05:70;
case 5
tD=0.1:0.05:110;
end
for i=1:1:N
QtD=QtD+log(2)./tD*V(i).*Laplace(log(2)*i./tD,reD);
end
loglog(tD,QtD)
title('无因次水侵曲线');
xlabel('无因次时间tD');
ylabel('Q(tD)');
hold on
end
for reD=6.0:1.0:10.0
QtD=0.0;
switch reD
case 6.0
tD=0.1:0.05:160;
case 7.0
tD=0.1:0.05:200;
case 8.0
tD=0.1:0.05:400;
case 9.0
tD=0.1:0.05:500;
case 10.0
tD=0.1:0.05:500;
end
for i=1:1:N
QtD=QtD+log(2)./tD*V(i).*Laplace(log(2)*i./tD,reD);
end
loglog(tD,QtD)
title('无因次水侵曲线');
xlabel('无因次时间tD');
ylabel('Q(tD)');
hold on
end