% gndw2.m Last mofified 8/4/2008 % This m-file displays the solution to the two layer aquifer problem % for radial symmetry. Data for the animation is stored in the mat-file % framedata.mat. % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); % working in meters and days here are the parameters bprime=5;Kprime=1.728;S=5e-5; K=60.48;b=10; T=K*b; B=sqrt(T*bprime/Kprime); Bsq4=B*B*4; % looks like the upper limit for t should be 200 days SoverT=S/T; r=(B/10)*linspace(1, 40,40); r=[B/40 r]; roverB=r/B; framedata=zeros(400,41); t=linspace(.5,200,400); load framedata; figure(1);%polt s vs. r/B for various values of t h1=axes('Ydir','reverse'); hold on index=[1 11 41 81 161]; for j=1:5 plot(roverB(1,1:21),framedata(index(1,j),1:21)); hold on end set(gca,'Box','on'); text(.5,.5,['B = ' num2str(B) ' meters']); text(.5,.45,'o o o Denotes Infinite Time Solution'); text(.2,.025,'t = 0.5'); text(.27,.05,'5'); text(.4,.07,'20'); text(.6,.07,'40'); xlabel('Dimensionless Distance from Well, r/B'); ylabel('Dimensionless Drawdown, s(r,t)T/Q'); assym=(1/(2*pi))*besselk(0,r/B); plot(roverB(1,1:21),assym(1,1:21),'o'); hold off text(.2,.65,'Press Enter to Continue'); pause figure(2);clf%plot s vs. t for various r h2=axes('Ydir','reverse'); hold on; index2=[1 2 3 5 10]; for i=1:5 plot(t,framedata(:,index2(1,i))); hold on end set(gca,'Box','on'); xlabel('Time, t [days]') ylabel('Dimensionless Drawdown, s(r,t)T/Q') text(120,.65,['B = ' num2str(B) ' meters']); text(20,.6,'r/B = 0.025') text(30,.4,'r/B = 0.1') text(35,.28,'r/B = 0.3') text(40,.17,'r/B = 0.5') text(50,.07,'r/B = 1') hold off text(15,.65,'Press Enter to Continue') pause t=[0 t]; framedata=[zeros(1,41);framedata]; figure(3);clf;% plot s vs r and t axis([0 4 0 200 0 0.6]); [R,T1]=meshgrid(roverB,t); mesh(R,T1,framedata); view(45,30) text(.5, 0,-.2,'Distance, r/B','Rotation',-20) text(4,50,-.2,'Time, t [days]', 'Rotation',20) zlabel('Dimensionless Drawdown, s(r,t)T/Q') hold off text(.5, 200, .7,'Press Enter to Continue','Rotation', -20) pause lim=max(max(framedata)); figure(4);clf;%plot animation of well drawdown as time evolves r1=[-fliplr(r) r]/B; framedata2=[fliplr(framedata) framedata]; h3=axes('Ydir','reverse'); axis([r1(1,1) r1(1,82) -.6*lim 3*lim]); hold on patch([r1(1,1) r1(1,82) r1(1,82) r1(1,1) r1(1,1)],[-.4*lim -.4*lim 0 0 -.4*lim],'m'); hold on text(1,-.3*lim,['B = ' num2str(B) ' m']); xpatch=[r1 r1(1,82) r1(1,1) r1(1,1)]; ypatch=[framedata2(1,:) lim lim framedata2(1,1)]; P1=patch(xpatch,ypatch,'c','EraseMode','background'); hold on patch([r1(1,1) r1(1,82) r1(1,82) r1(1,1) r1(1,1)],[lim lim 3*lim 3*lim lim],'b'); hold on plot([0 0 .5],[2.9*lim -.5*lim -.5*lim],'k','LineWidth',5); hold on; P2=plot([.5 .9 1.3],[-.5*lim -.45*lim -.42*lim],'b','LineWidth',5, 'EraseMode','background'); hold on text(-3.5,-.1*lim,'Permeable Layer') text(-3.5,.4*lim,'Clay-Aquitard') text(-3.5,1.4*lim,'Aquifer') T1=text(-3.5,-.3*lim,'Press Enter to Continue'); pause set(T1,'String',' '); for i=1:401 ypatch=[framedata2(i,:) 1*lim 1*lim framedata2(i,1)]; set(P1,'Ydata',ypatch); plot([0 0],[2.8*lim -.5*lim],'k','LineWidth',5); hold on pause(.02) end index2=[1 6 20 40]; for i=1:4 plot(r1,framedata2(index(1,i),:)); hold on end xlabel('Radial Distance from Well, r/B') ylabel('Dimensionless Drawdown, s(r,t)T/Q') text(-3.5,-.1*lim,'Permeable Layer') text(1,.3,['Kprime = ' num2str(Kprime) ' m/day']) text(1,.2,['bprime = ' num2str(bprime) ' m']) text(1,2*lim,['S = ' num2str(S)]) text(1,1.6*lim,['b = ' num2str(b),' m']) text(1,1.8*lim,['K = ' num2str(K) ' m/day']) hold off set(T1,'String','Press Enter to Continue'); set(gca,'Box','on'); pause figure(5):clf;% plot 3-d plot of steady-state s(t) for i=1:40 for j=1:40 r=sqrt(i^2+j^2)*.1; w(i,j)=(1/(2*pi))*besselk(0,r); end end y1=linspace(.1,4,40); ax=(1/(2*pi))*besselk(0,y1); w1=flipud(w); w2=fliplr(w); w3=fliplr(w1); ax1=fliplr(ax); ax2=ax'; ax3=flipud(ax2); W=[w3 ax3 w1;ax1 ax(1,1) ax;w2,ax2 w]; x=linspace(-4,4,81); y=x; [X, Y]=meshgrid(x,y); h5=axes('Zdir','reverse'); hold on mesh(X,Y,W); text(2,4,.35,'Distance, x/B','Rotation',33); text(4,-2,.35,'Distance, y/B', 'Rotation',-11); view(120,30); zlabel('Drawdown, sT/Q') grid on text(-4,-1,.35,'Press Enter to Continue','Rotation',-11) text(4,-3, .3,'Steady-State Cone of Depression','Rotation',-11) hold off pause