% gndw3.m Last modified 8/4/2008 % Solves the Theis model for a single layer confined aquifer with % radial symmetry. % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. % Copyright Raymond G. Jacquot amd Thomas V. Edgar clear all % working in meters and days here are the parameters K=8.7;b=5.30; T=K*b;S=2e-5; r=linspace(1,100,41); framedata=zeros(400,41); t=linspace(.005,2,400); Sover4T=S/(4*T); u=zeros(800,41); a=[-.57721566 .99999193 -.24991055 .05519968 -.00976004 .00107857]; for j=1:400% time loop for i=1:41 % space loop u(j,i)=(r(1,i)*r(1,i)*Sover4T)/(t(1,j)); u1=u(j,i); E1=a(1,1)+a(1,2)*u1+a(1,3)*(u1^2)+a(1,4)*(u1^3)+a(1,5)*(u1^4)+a(1,6)*(u1^5)-log(u1); I=(1/(4*pi))*(E1); framedata(j,i)=I; end end figure(1);clf;%polt s vs. r for various values of t h1=axes('Ydir','reverse'); hold on index=[1 121 241 400 ]; for j=1:40:361 plot(r,framedata(j,:)); hold on; end box on xlabel('Radial Distance from Well, r [m]') ylabel('Dimensionless Drawdown, s(r,t)T/Q') text(60,1.2,['b = ' num2str(b) ' m']); text(60,1.1,['K = ' num2str(K) ' m/day']); text(60,1.,['S = ' num2str(S)]); t1=5; text(10,1.3,'Press Enter to Continue') pause hold off figure(2);clf;%plot s vs. t for various r h2=axes('Ydir','reverse'); axis([0 2 0 1.4]); hold on index2=[1 3 5 10 15]; for i=1:4:41 plot(t,framedata(:,i)) hold on end box on text(1.4,.3,['b = ' num2str(b) ' m']); text(1.4,.2,['K = ' num2str(K) ' m/day']); text(1.4,.1,['S = ' num2str(S)]); xlabel('Time, t [days]') ylabel('Dimensionless Drawdown, s(r,t)T/Q') text(.5,1.3,'Press Enter to Continue') hold off pause t=[0 t]; framedata=[zeros(1,41);framedata]; figure(3);clf;% plot s vs r and t axis([0 100 0 4 0 0.6]); [R,T1]=meshgrid(r,t); mesh(R,T1,framedata); view(45,30) xlabel('Distance from Well, r [m]','Rotation',-20) ylabel('Time, t [days]', 'Rotation',20) zlabel('Dimensionless Drawdown, s(r,t)T/Q') text(20,2,1.4,'Press Enter to Continue','Rotation',-20) pause lim=max(max(framedata)); figure(4);clf;%plot animation of well drawdown as time evolves h3=axes('Ydir','reverse'); r1=[-fliplr(r) r]; framedata2=[fliplr(framedata) framedata]; axis([r1(1,1) r1(1,82) -1.1*lim 1.5*lim]) hold on patch([r1(1,1) r1(1,82) r1(1,82) r1(1,1) r1(1,1)],[-.8*lim -.8*lim 0 0 -.8*lim],'m'); hold on xpatch=[r1 r1(1,82) r1(1,1) r1(1,1)]; ypatch=[framedata2(1,:) 1.5*lim 1.5*lim framedata2(1,1)]; P1=patch(xpatch,ypatch,'b','EraseMode','Background'); hold on plot([0 0],[1.4*lim -lim],'k','LineWidth',5);%well pipe plot([-2 20],[-lim -lim],'k','LineWidth',5)%horiz. well pipe plot([20 25 30],[-lim -.95*lim -.85*lim],'b','Linewidth',5)% water spout hold on text(-80,-.4*lim,'Impermeable layer') text(-80,lim,'Aquifer') T1=text(-80,-.7*lim,'Press Enter to Animate'); pause set(T1,'String',' ') for i=1:401 ypatch=[framedata2(i,:) 2*lim 2*lim framedata2(i,1)]; set(P1,'Ydata',ypatch); plot([0 0],[1.9*lim -lim],'k','LineWidth',5); hold on pause(.02) end text(-80,.3*lim,['t = ' num2str(t(1,401)) 'days']); xlabel('Radial Distance from Well, r [m]') ylabel('Dimensionless Drawdown, s(r,t)T/Q') text(40,lim,['b = ' num2str(b) ' m']); text(40,1.15*lim,['K = ' num2str(K) ' m/day']); text(40,1.3*lim,['S = ' num2str(S)]); set(T1,'String','Press Enter to Continue'); set(gca,'Box','on') pause figure(5);clf; for j=1:40:401 semilogx(r,framedata(j,:)); hold on; end xlabel('Radial Distance from Well, r [m]') ylabel('Dimensionless Drawdown, s(r,t)T/Q') text(2,.2,'Press Enter to Continue') pause