一维地下水运动MATLAB模拟.docx
- 文档编号:7535773
- 上传时间:2023-01-24
- 格式:DOCX
- 页数:14
- 大小:166.98KB
一维地下水运动MATLAB模拟.docx
《一维地下水运动MATLAB模拟.docx》由会员分享,可在线阅读,更多相关《一维地下水运动MATLAB模拟.docx(14页珍藏版)》请在冰豆网上搜索。
一维地下水运动MATLAB模拟
Problem1
Thegroundflowequationfor1Dheterogeneous,isotropicporousmediumwithaconstantaquiferthicknessisgivenby:
Whereh(x,t)isthehydraulichead,T(x)isthetransmissivity,andSisthestorativity.TheboundaryconditionsimposedareconstantheadhL(15m)andhR(5m)attheleftandrightendsofsoilcolumn,respectively.Theinitialconditionis0orrandomnumberateachnode.Thelengthofaquiferis1m(L=1m).Thestorativityis1(S=1).
DevelopaMATLABprogramwhichcanhandleaheterogeneoustransmissivityfieldusingtheimplicitmethod.
1)Testtheprogramforthecaseofhomogenousparameters,andcomparetheresultswiththeresultsgeneratedfromflowequationwithhomogenoustransmissivityfield(i.e.,
)
2)Testtheprogramforpatch-heterogeneousTfield,whichcomprisesatwo-zoneTfieldwithafivefolddifferenceintransmissivity(T1=5T2,T1=0.01m2s-1,T2=0.002m2s-1;or,T1=0.2T2,T1=0.002m2s-1,T2=0.01m2s-1)andaninterfaceinthemiddleofthedomain.
Problem1解答
1理论推导
非均质一维地下水运动方程为:
如右图,各个点的水头分别为
,水头
和
之间水力传导系数为
。
将方程
(1)进行离散化:
设
则
设
则
所以:
设
则
2编码实现
根据第一部分的理论推导编写代码
2.1均质情况下一维地下水运动
functionplot_Data=oneDimGroudwaterFlowHom(hL,hR,L,S,T,method,dx,dt)
%%Finitedifferencemethodtosolve1-Dflowequation
%%WritedByDongdongKong,2014-01-07
%SunYat-SenUniversity,Guangzhou,China
%emal:
%------------------------------------------------------
%onedimensiongroundwaterflowmodelscript
%whichcanhandleaheterogeneoustransmissivityfield
%------------------------------------------------------
%becauseofTvalue,thisfunctiononlysuitforpatch-heterogeneousT
%fieldcomprisesatwo-zoneTfield,orhomogeneousfield.youcanModify
%inputdiscretedTvaluetosuitotherheterogeneoussituation.
%------------------------------------------------------
%hL%leftconstanthydraulichead(m)
%hR%leftconstanthydraulichead(m)
%L%thelengthofaquifer
%S%storativity
%T%transmissivity
%method%'forward'or'backward'approximation
%------------------------------------------------------
x=0:
dx:
L;nx=length(x);
w=T/S*dt/dx/dx;
h=rand(nx,1);
h
(1)=hL;h(end)=hR;
times=100;
Result=zeros(nx,times+1);
ifstrcmp(method,'forward')
%%forwardapproximation
W=(1-2*w)*diag(ones(nx,1),0)+...
w*diag(ones(nx-1,1),1)+...
w*diag(ones(nx-1,1),-1);
W(1,1)=1;W(1,2)=0;
W(nx,nx)=1;W(nx,nx-1)=0;
fori=1:
times+1
Result(:
i)=h;
h=W*h;
end
elseifstrcmp(method,'backward')
%%backwardapproximation
w=T/S*dt/dx/dx;
W=(1+2*w)*diag(ones(nx,1),0)-...
w*diag(ones(nx-1,1),1)-...
w*diag(ones(nx-1,1),-1);
W(1,1)=1;W(1,2)=0;
W(nx,nx)=1;W(nx,nx-1)=0;
fori=1:
times+1
Result(:
i)=h;
h=W\h;
end
else
warning('errormethodinput,methodcanbeonlyset''forward'',''backward''')
end
figure
timeID=[015102040100]+1;
plot_Data=Result(:
timeID);
plot(x,plot_Data,'linewidth',1.5)
legend({'t=0','t=1','t=5','t=10','t=20','t=40','t=100'})
2.2非均质情况下一维地下水运动
functionplot_Data=oneDimGroudwaterFlowHet(hL,hR,L,S,T)
%%Finitedifferencemethodtosolve1-Dflowequation
%%WritedByDongdongKong,2014-01-07
%SunYat-SenUniversity,Guangzhou,China
%emal:
%------------------------------------------------------
%onedimensiongroundwaterflowmodelscript
%whichcanhandleaheterogeneoustransmissivityfield
%------------------------------------------------------
%becauseofTvalue,thisfunctiononlysuitforpatch-heterogeneousT
%fieldcomprisesatwo-zoneTfield,orhomogeneousfield.youcanModify
%inputdiscretedTvaluetosuitotherheterogeneoussituation.
%------------------------------------------------------
%hL%leftconstanthydraulichead(m)
%hR%leftconstanthydraulichead(m)
%L%thelengthofaquifer
%S%storativity
%T%transmissivity
%------------------------------------------------------
nx=50;
x=linspace(0,L,nx*2+1);
dx=x
(2)-x
(1);
dt=0.1;
N=length(x);
Ti=[repmat(T
(1),1,nx),repmat(T
(2),1,nx)];
W=zeros(N);
w=dt/dx/dx/S;
fori=2:
N-1
W(i,i+1)=-w*Ti(i);
W(i,i)=1+w*(Ti(i)+Ti(i-1));
W(i,i-1)=-w*Ti(i-1);
end
W(1,1)=1;
W(end,end)=1;
h=rand(N,1);
h
(1)=hL;h(end)=hR;
times=10000;%simualtiontimes
Result=zeros(N,times+1);
figure,holdon
fori=1:
times+1
Result(:
i)=h;
h=W\h;
end
timeID=[0151020401001000]*10+1;
plot_Data=Result(:
timeID);
plot(x,plot_Data,'linewidth',1.5)
legend({'t=0','t=1','t=5','t=10','t=20','t=40','t=100','t=1000'})
2.3主函数
clear,clc
%%WritedByDongdongKong,2014-12-30
%onedimensiongroundwaterflowmodelscript
%whichcanhandleaheterogeneoustransmissivityfield
hL=15;%leftconstanthydraulichead(m)
hR=5;%leftconstanthydraulichead(m)
L=1;%thelengthofaquifer
S=1;%storativity
%#question1.1
dt=1;
dx=0.2;%fordelta_w<0.5
T=0.01;method='forward';
Hom_forward=oneDimGroudwaterFlowHom(hL,hR,L,S,T,method,dx,dt);
dx=0.01;
T=0.01;method='backward';
Hom_backward=oneDimGroudwaterFlowHom(hL,hR,L,S,T,method,dx,dt);
T=[0.010.01];
Hom=oneDimGroudwaterFlowHet(hL,hR,L,S,T);
%#question1.2
T=[0.010.002];
Het_situ1=oneDimGroudwaterFlowHet(hL,hR,L,S,T);
T=[0.0020.01];
Het_situ2=oneDimGroudwaterFlowHet(hL,hR,L,S,T);
save('oneDimFlowPlotData.mat','Hom_forward','Hom_backward','Hom','Het_situ1','Het_situ2')
3.结果输出
Question1.1result:
由图一可以看出非均质插分和均质情况下向前插分、向后插分,计算结果和收敛速度大致相同,考虑到向前插分中
,因此向前差分时取
,其他情况取
。
图一.非均质与均质求解下的差异,从左到右分别为非均质T插分、均质情况下向后插分、均质情况下向前插分
Question1.2result:
T1=0.01
T2=0.002
和T1=0.2T2,T1=0.002
T2=0.01
非均质T情况下插分结果如图二。
图二.非均质情况下两个情景下插分结果
Problem2
Userainfallfrom8stationstoestimatetheannualrainfallfor6stationsusingthekrigingordinarymethods.
Thelocationsofthe8stationswhererainfallisavailableareasfollows:
No.
stationname
latitude
logitude
1
boulder
40.033
105.267
2
castle
39.383
104.867
3
green9ne
39.267
104.750
4
lawson
39.767
105.633
5
longmont
40.252
105.150
6
manitou
38.850
104.933
7
parker
39.533
104.650
8
woodland
39.100
105.083
The15-yearrainfallrecordsfor8thestationsareavailableintheExcelfile“Rain_8Stations_Known.xls”
Thelocationsofthe6stationswhererainfallneedtobeestimatedareasfollows:
No.
stationname
latitude
logitude
1
byers
39.750
104.133
2
estes
40.383
105.517
3
golden
39.700
105.217
4
green9se
39.100
104.733
5
lakegeor
38.917
105.483
6
morrison
39.650
105.200
Theannualrainfallfor6thestationsareavailableintheExcelfile“AnnualRain_6Stations_ToBeEstimated.xls”.Thiswillbeusedtocalculatetheestimateerrorbyusingdifferentinterpolationmethods
1)Foreachof6stations(i.e.,byers,estes,golden,green9se,lakegeor,andmorrison)wheretheannualwillbeestimated,estimatetheannualrainfallbasedonrainfallat8stations(i.e.,boulder,castle,green9ne,lawson,longmont,manitou,parkerandwoodland)usingthekrigingmethod.
2)Estimatetheannualrainfallusingtheinverse-distance-squaremethodandthearithmeticmeanmethod,andcomparetheerrorfromthekrigingmethod.
3)Calculatetheerrorvarianceforkrigingmethodateachof6stations
Note:
useanexponentialfunctiontomodelcovariance:
whereaandbareparameterstobedetermined,anddisthedistancebetweentwopoints.
Problem2解答
1.KrigingMATLAB编码如下
KrigineInterpolate.m
clear,clc
%%WritedByDongdongKong,2014-01-08
%SunYat-SenUniversity,Guangzhou,China
%email:
%------------------------------------------------------
dataDir='./data/';
stationInfo_known=xlsread([dataDir,'LatitudeLogitude_8Stations_Known.xls']);
stationInfo_pre=xlsread([dataDir,'LatitudeLogitude_6Stations_ToBeEstimated.xls']);
Rain_known=xlsread([dataDir,'Rain_15year_8Stations_Known.xls']);
Rain_preReal=xlsread([dataDir,'AnnualRain_6Stations_ToBeEstimated.xls']);
stationInfo_known=stationInfo_known(:
[3,4]);
stationInfo_pre=stationInfo_pre(:
[3,4]);
nYear=size(Rain_known,1);
Np_know=8;%stationnumberofknowdata
Np_pre=6;%stationnumbertopredict
RainPredict=zeros(nYear,Np_pre);
fork=1:
Np_pre
stationInfo=[stationInfo_pre(k,:
);stationInfo_known];
nStation=size(stationInfo,1);
dispPoint=zeros(nStation);
%calculatetwostationdistance
fori=1:
nStation
dispPoint(:
i)=distance(stationInfo,stationInfo(i,:
));
end
%covariance
mat_cor=cov(Rain_known);
%simulateCovfunction
x=dispPoint(2:
end,2:
end);
y=mat_cor;
x_tri=tril(x);y_tri=tril(y);
X=x_tri(:
);ind=X~=0;
Y=y_tri(:
);
%deletediagzerosdata
X_nozeros=X(ind);Y_nozeros=Y(ind);
Y_zeros=diag(y);
%GeneralmodelExp1:
%y=a*exp(b*x)
[fitResult,gof]=expFit(X_nozeros,Y_nozeros);
a=fitResult.a;
b=fitResult.b;
C0=mean(Y_zeros)-b;
Cov_fit=reshape(feval(fitResult,x),8,8);
fori=1:
length(Cov_fit)
Cov_fit(i,i)=C0+b;
end
C=[Cov_fit,ones(8,1);ones(1,8),0];
dh=dispPoint(2:
end,1);
D=[feval(fitResult,dh);1];
W=C\D;%lastelementismu
w=W(1:
end-1);
RainPredict(:
k)=Rain_known*w;
end
Rain_kriging=mean(RainPredict)
expFit.m
function[fitresult,gof]=expFit(X1,Y1)
%MATLABAUTOGENERATECODEexpFitFUNCTION
[xData,yData]=prepareCurveData(X1,Y1);
%Setupfittypeandoptions.
ft=fittype('exp1');
opts=fitoptions(ft);
opts.Display='Off';
opts.Lower=[-Inf-Inf];
opts.StartPoint=[871049.657333104-6.876];
opts.Upper=[InfInf];
%Fitmodeltodata.
[fitresult,gof]=fit(xData,yData,ft,opts);
%----------------------------------------------------------------
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 地下水 运动 MATLAB 模拟