一维地下水运动MATLAB模拟.doc

(word完整版)一维地下水运动MATLAB模拟
Error! Bookmark not defined.Problem 1
The ground flow equation for 1D heterogeneous, isotropic porous medium with a constant aquifer thickness is given by:

Where h(x，t） is the hydraulic head， T(x） is the transmissivity， and S is the storativity. The boundary conditions imposed are constant head hL （15m) and hR (5m) at the left and right ends of soil column， respectively. The initial condition is 0 or random number at each node. The length of aquifer is 1m （L=1m）. The storativity is 1 （S = 1）。

Develop a MATLAB program which can handle a heterogeneous transmissivity field using the implicit method。
1) Test the program for the case of homogenous parameters, and compare the results with the results generated from flow equation with homogenous transmissivity field (i.e., )
2) Test the program for patch- heterogeneous T field， which comprises a two—zone T field with a fivefold difference in transmissivity (T1=5T2, T1=0.01m2s—1， T2=0.002m2s-1 ; or, T1=0。2T2, T1=0。002m2s—1, T2=0。01m2s-1) and an interface in the middle of the domain。

Problem 1解答
1 理论推导
非均质一维地下水运动方程为： \＊ MERGEFORMAT (Error! Bookmark not defined.)
如右图，各个点的水头分别为，水头和之间水力传导系数为.
将方程(1)进行离散化：设，则
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.)
（Error! Bookmark not defined.)
\＊ MERGEFORMAT (3）
Error! Bookmark not defined.(Error! Bookmark not defined.）
设 则
(Error! Bookmark not defined.）
Error! Bookmark not defined.(Error! Bookmark not defined.）
所以:
(5)
设， 则

， ,
\＊ MERGEFORMAT Error! Bookmark not defined.（Error! Bookmark not defined.）
2 编码实现
根据第一部分的理论推导编写代码
2.1 均质情况下一维地下水运动
function plot_Data = oneDimGroudwaterFlowHom(hL，hR，L，S,T,method,dx，dt）
%％ Finite difference method to solve 1—D flow equation
％% Writed By Dongdong Kong, 2014—01—07
% Sun Yat-Sen University， Guangzhou， China
％ emal： [e-mail]。
％ -—--—-————---——-—--—--—-——-—-——---—----————-—----—-——-
% one dimension groundwater flow model script
% which can handle a heterogeneous transmissivity field
% ——-—--————----—-—--———-————-----——--—-——-———-——-—--—-—
% because of T value, this function only suit for patch-heterogeneous T
% field comprises a two—zone T field， or homogeneous field。 you can Modify
% input discreted T value to suit other heterogeneous situation.
% ——-—-—--————--—-————-————--—-———---——----——----—--——--
% hL ％ left constant hydraulic head (m)
％ hR % left constant hydraulic head （m）
％ L ％ the length of aquifer
% 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）;

if strcmp（method,’forward'）
％% forward approximation
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;

for i=1：times+1
Result（:，i）=h；
h=W＊h；
end
elseif strcmp(method,'backward’)
％％ backward approximation
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；

for i=1:times+1
Result（:,i)=h;
h=W\h；
end
else
warning('error method input, method can be only set ’’forward’’,’'backward'’’)
end
figure
timeID=[0 1 5 10 20 40 100]+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 非均质情况下一维地下水运动
function plot_Data = oneDimGroudwaterFlowHet（hL， hR, L， S， T）
％% Finite difference method to solve 1-D flow equation
%％ Writed By Dongdong Kong， 2014—01-07
％ Sun Yat-Sen University， Guangzhou, China
％ emal: kongdd＠
% ——————---—-—-------———-——-—-—-----———---————-—————----
% one dimension groundwater flow model script
％ which can handle a heterogeneous transmissivity field
% -—-——————-—----—-———--—--—------——--——--——-—---—-————-
% because of T value, this function only suit for patch—heterogeneous T
％ field comprises a two—zone T field, or homogeneous field. you can Modify
％ input discreted T value to suit other heterogeneous situation.
％ —-—----—-—-------—-—----———-----——-----—————--——-—————
％ hL ％ left constant hydraulic head （m)
％ hR ％ left constant hydraulic head (m)
％ L % the length of aquifer
% 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；
for i=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； %simualtion times
Result=zeros（N,times+1）；

figure，hold on
for i=1：times+1
Result(:,i）=h;
h=W\h;
end
timeID=[0 1 5 10 20 40 100 1000］*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
％% Writed By Dongdong Kong， 2014-12-30
％ one dimension groundwater flow model script
% which can handle a heterogeneous transmissivity field
hL = 15； % left constant hydraulic head (m）
hR = 5; % left constant hydraulic head (m）
L = 1； % the length of aquifer
S = 1; % storativity
% ＃ question1.1
dt=1；
dx=0.2; %for delta_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。01 0.01］;
Hom = oneDimGroudwaterFlowHet(hL，hR，L,S，T）;
% ＃ question1。2
T=[0。01 0.002］；
Het_situ1 = oneDimGroudwaterFlowHet（hL,hR，L，S,T）;
T=[0.002 0.01]；
Het_situ2 = oneDimGroudwaterFlowHet(hL，hR，L，S，T）；
save（’oneDimFlowPlotData.mat’,’Hom_forward'，’Hom_backward’，’Hom',’Het_situ1','Het_situ2')

3。 结果输出
Question1.1 result：
由图一可以看出非均质插分和均质情况下向前插分、向后插分,计算结果和收敛速度大致相同，考虑到向前插分中，因此向前差分时取,其他情况取。
图一。 非均质与均质求解下的差异,从左到右分别为非均质T插分、均质情况下向后插分、均质情况下向前插分

Question1。2 result：
T1=0.01, T2=0。002 和 T1=0.2T2， T1=0。002， T2=0。01非均质T情况下插分结果如图二.
图二. 非均质情况下两个情景下插分结果

Problem 2
Use rainfall from 8 stations to estimate the annual rainfall for 6 stations using the kriging ordinary methods.
The locations of the 8 stations where rainfall is available are as follows：
No。
station name
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

The 15-year rainfall records for 8 the stations are available in the Excel file “Rain_8Stations_Known。xls”
The locations of the 6 stations where rainfall need to be estimated are as follows：
No。
station name
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

The annual rainfall for 6 the stations are available in the Excel file “AnnualRain_6Stations_ToBeEstimated。xls”. This will be used to calculate the estimate error by using different interpolation methods

1) For each of 6 stations (i。e.， byers, estes， golden， green9se, lakegeor， and morrison） where the annual will be estimated， estimate the annual rainfall based on rainfall at 8 stations （i。e。, boulder， castle， green9ne, lawson, longmont, manitou， parker and woodland) using the kriging method.
2) Estimate the annual rainfall using the inverse-distance-square method and the arithmetic mean method, and compare the error from the kriging method.
3) Calculate the error variance for kriging method at each of 6 stations

Note: use an exponential function to model covariance: ， where a and b are parameters to be determined, and d is the distance between two points。

Problem 2解答
1. Kriging MATLAB 编码如下
KrigineInterpolate.m
clear，clc
％% Writed By Dongdong Kong， 2014-01—08
% Sun Yat—Sen University， Guangzhou， China
％ email: kongdd＠
% -——--——--———---———--—-—-——-———-——----———-—-—-—-————-——
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; % station number of know data
Np_pre = 6; % station number to predict
RainPredict = zeros（nYear, Np_pre)；
for k = 1:Np_pre
stationInfo = [stationInfo_pre（k，：)； stationInfo_known］；
nStation = size(stationInfo，1）；
dispPoint=zeros(nStation）;
％ calculate two station distance
for i=1:nStation
dispPoint（:，i） = distance(stationInfo, stationInfo(i,:))；
end
% covariance
mat_cor = cov(Rain_known)；
% simulate Cov function
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（:)；
％ delete diag zeros data
X_nozeros = X（ind）； Y_nozeros = Y(ind)；
Y_zeros = diag(y);
％ General model Exp1：
% 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)；
for i=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; %last element is mu
w = W（1:end—1）;
RainPredict（：，k） = Rain_known＊w;
end
Rain_kriging = mean(RainPredict）
expFit。m
function ［fitresult， gof] = expFit（X1， Y1)
% MATLAB AUTO GENERATE CODE expFit FUNCTION
[xData， yData］ = prepareCurveData( X1， Y1 ）；
％ Set up fittype and options.
ft = fittype（ ’exp1' );
opts = fitoptions( ft );
opts。Display = ’Off’；
opts。Lower = [-Inf -Inf］;
opts。StartPoint = ［871049.657333104 —6.29140824382876］;
opts.Upper = [Inf Inf];
% Fit model to data。
［fitresult， gof］ = fit( xData， yData, ft， opts )；
％ ——-————-——--—-——--—-—-——--———-——-—--——---—-——--——-———-----—-—--—
计算结果保留在第二问的表格里。
2. 反距离权重插值、算术平均值和kriging插值的误差比较
反距离权重插值采用R语言下的工具包gstat的idw函数进行，对于每个站点采用8个雨量站进行插值。
IDW Interpolation R 代码如下：
library(gstat)
library(maptools)
#＃ Loading required package： sp
＃# Checking rgeos availability: TRUE
rm(list=ls（））
setwd（"F:/Users/kongdd/Documents/MATLAB/finalWork/krigineInterpolation"）

stationInfo_known=read.table('。/data/LatitudeLogitude_8Stations_Known.txt'，sep=”\t”,head=T)
stationInfo_pre = read。table（'./data/LatitudeLogitude_6Stations_ToBeEstimated。txt',sep=”\t”,head=T)
Rain_known = read。table(’。/data/Rain_15year_8Stations_Known。txt’,sep="\t”,head=T)
Rain_preReal = read.table('./data/AnnualRain_6Stations_ToBeEstimated.txt'，head=T)
Rain_known <- t(Rain_known［-1］)
colnames（Rain_known）〈—paste(”year",1:15,sep=””)
stationInfo_known = stationInfo_known[，c（3,4)]
stationInfo_pre = stationInfo_pre[，c(3,4）]

Rain_known<- data。frame(stationInfo_known,Rain_known)
coordinates（stationInfo_pre)<-～long+lat
coordinates(Rain_known)<—～long+lat

nYear <— 15
Result〈—list（）
RainPredict<—matrix(0,nYear，6)
for(i in 1:nYear）{
varName <- sprintf(”year%d”，i)
evalStr<- sprintf（"Result［［％d]] 〈— idw(％s～1,Rain_known,stationInfo_pre,nmin=8)",i，varName）
eval（parse（text=evalStr）)
RainPredict[i，］ 〈— Result[［i］]$var1。pred
}
AnnualRainPredict<—apply（RainPredict,2,mean)
station
byers
estes
golden
green9se
lakegeor
morrison
Real
8916.3
7949.7
9054.9
8859。6
8276.5
7743。4
kriging
8416。6
8259。8
8262。8
8464.5
8383.9
8283.5
idw
8531.8
8191。1
8233。6
8515.3
8460。3
8277。9
Mean
8363。9
8363。9
8363.9
8363.9
8363。9
8363.9
表一. 不同插值方法的年均降雨插值结果

station
byers
estes
golden
green9se
lakegeor
morrison
kriging
—499.69
310.14
—792.04
—395。11
107。42
540。11
idw
115.23
-68。78
-29。20
50。76
76.35
—5.57
Mean
—167。94
172.83
130。24
-151.38
—96.41
85。99
表一. 年均降雨不同插值方法的相对误差
3. 计算kriging插值的误差的方差
， 其中代表真实值、代表预测值，n为样本长度。Var= 268791.14