% asg 12 1-D, 1-Phase pressure response in linear system % Constant flux in finite system; Carslaw and Jaeger clear all perm=50. ; por=0.3 ; vis=1.0 ; c=2.e-5 ; Lx=1000.; dy=100. ; dz=10. ; b=1. ; q=-10. ; pic=1000. ; kap=0.00633*perm/(por*c*vis); F=5.615*q/(por*c*Lx*dy*dz); FK=887.2*q*vis/(dy*dz*perm*b); % profiles t=[1. 3. 10. 30.]; x=linspace(0,1,10); x=Lx.*x; p=zeros(length(t),length(x)); for i=1:length(t) pavg=pic+F*t(i); for j=1:length(x) steady=(3*x(j)^2-Lx^2)/(6*Lx^2); sum=0.; for n=1:40 ex=-kap*(pi*n/Lx)^2*t(i); arg=pi*n*x(j)/Lx; sum=sum+(((-1)^n)/n^2)*exp(ex)*cos(arg); end p(i,j)=pavg+FK*Lx*(steady-(2*pi^(-2))*sum); end end figure(1) clf reset plot(x,p) % Analytical solution title('Pressure Profile, t=1, 3, 10, 30 days') xlabel('Distance, feet') ylabel('Pressure, psi') NX=input('NX of F.D.? 0 if none. ') if (NX~=0) hold on clear x p load prof1.dat x=prof1(1:NX,1); for i=1:4 ist=1+NX*(i-1); ilst=NX*i; p=prof1(ist:ilst,2); plot(x,p,'kx') %Numerical solution end end % Pressure history figure(2) clf reset if NX>0 dx=Lx/NX; else dx=0; end clear t p t=linspace(0,50); for i=1:length(t) pavg=pic+F*t(i); steady=1/3; sum=0.; for n=1:40 ex=-kap*(pi*n/Lx)^2*t(i); arg=pi*n; sum=sum+(((-1)^n)/n^2)*exp(ex)*cos(arg); end p(i)=pavg+FK*Lx*(steady-(2*pi^(-2))*sum); end plot(t,p) % Analytical solution title('Pressure History, x=Lx') xlabel('Time, days') ylabel('Pressure, psi') if (NX~=0) hold on load prod1.dat tfd=prod1(:,1); pfd=prod1(:,4); plot(tfd,pfd,'kx') % Numerical solution end figure(3) % Plot square root of time clf reset sqt=zeros(length(t)); sqt=sqrt(kap.*t.*Lx^(-2)); plot(sqt,p) % Analytical solution axis([0 1.5 600 1000]) title('Pressure History, x=Lx') ylabel('Pressure, psi') xlabel('sqrt(kappa*t/Lx**2)') if (NX~=0) hold on sqtfd=sqrt(kap.*tfd.*Lx^(-2)); plot(sqtfd,pfd,'kx') % Numerical solution end