% asg 6.4 & 6.5: 2-D, 1-Phase pressure response in finite system % Well in finite system, Matthews and Russell % Radial symmetry clear all clf reset perm=50. ; por=0.3 ; vis=1.0 ; c=2.e-5 ; Lx=1000.; Ly=1000.; Lz=1.0 ; rw=0.4 ; re=1128.38; red=re/rw; redf=(3*red^4-4*red^4*log(red)-2*red^2-1)... /(4*(red^2-1)^2); b=1. ; q=-10. ; q=4*q ; % If rate is for one quadrant, then multiply by factor of 4 pic=1000. ; kap=0.00633*perm/(por*c*vis); FK=887.2*q*vis/(perm*b*Lz)*(1/(2*pi)); % profiles t=[1. 3. 10. 30. 50.]; load alf.dat [m,n]=size(alf); rd=alf(1,2:n); r=zeros(size(rd)); r(:)=rw*rd(:); alpha=alf(2:m,1); bess=alf(2:m,2:n); p=zeros(length(t),length(r)); for i=1:length(t) td=kap*t(i)/rw^2; for j=1:length(r) sum=0.; nsum=0; term=1.0; while (abs(term) > 1.e-7*abs(sum) & nsum < m) nsum=nsum+1; term=exp(-alpha(nsum)^2*td)*bess(nsum,j); sum=sum+term; end p(i,j)=pic+FK*((2/(red^2-1))*(rd(j)^2/4+td)... -(red^2*log(rd(j)))/(red^2-1)-redf... +pi*sum); end end figure(1) clf reset semilogx(r,p) axis([10 2000. 0. 1000.]) title('Pressure Profile, t=1, 3, 10, 30, 50 days') xlabel('Radial Distance, feet') ylabel('Pressure, psi') NX=input('NX of F.D.? 0 if none. ') if (NX>0) hold on clear r p load prof1.dat r=sqrt(2)*prof1(1:NX,1); % Radial distance along diagonal for i=1:5 ist=1+NX*(i-1); ilst=NX*i; p=diag(prof1(ist:ilst,2:(NX+1))); % Pressure along diagonal semilogx(r,p,'kx') end end if(NX>0) % Contour plots of simulation figure(2) clf reset clear x p v=0:50:1000 ; tim=[1 3 10 30] ; load prof1.dat x=prof1(1:NX,1); y=x ; p=zeros(length(x),length(y)); for i=1:4 ist=1+NX*(i-1); ilst=NX*i; for j=1:length(y) p(:,j)=prof1(ist:ilst,j+1); end subplot(2,2,i) contour(x,y,p,v) title(['50 psi intv., t = '... ,int2str(tim(i)),' days']) % xlabel('Distance, feet') % ylabel('Distance, feet') axis('square') axis('equal') axis([0 1000 0 1000]) end % Pressure history % expression for quarter well in corner of grid block figure(3) clear t p clf reset t=logspace(-2,log10(50),100); pinf=zeros(size(t)); for i=1:length(t) r=exp(-0.5)*(2*Lx/NX)/sqrt(pi); arg=r^2/(4*kap*t(i)); pinf(i)=pic-FK*0.5*ei(-arg); % Infinite reservoir end plot(t,pinf,'b-') hold on % Use image wells to represent a confined well in a square psq=zeros(size(t)); for i=1:length(t) arg=(2*Lx)^2/(4*kap*t(i)); arg2=(2*sqrt(2)*Lx)^2/(4*kap*t(i)); arg3=(4*Lx)^2/(4*kap*t(i)); arg4=(2*sqrt(5)*Lx)^2/(4*kap*t(i)); arg5=(2*sqrt(5)*Lx)^2/(4*kap*t(i)); arg6=(4*sqrt(2)*Lx)^2/(4*kap*t(i)); psq(i)=pinf(i)-4*FK*0.5*ei(-arg)-4*FK*0.5*ei(-arg2)... -4*FK*0.5*ei(-arg3)-4*FK*0.5*ei(-arg4) ... -4*FK*0.5*ei(-arg5)-4*FK*0.5*ei(-arg6) ... ; % Well in square end plot(t,psq,'g-') % Finite radial reservoir clear t p rd prof1 r=exp(-0.5)*(2*Lx/NX)/sqrt(pi); rd=r/rw; t=logspace(-2,log10(50),100); p=zeros(size(t)); j=14; % j should be between 13 and 14, to be exact % i.e., 151 < rd = 68/0.4 < 229 for i=1:length(t) td=kap*t(i)/rw^2; sum=0.; nsum=0; term=1.0; while (abs(term) > 1.e-7*abs(sum) & nsum < m) nsum=nsum+1; term=exp(-alpha(nsum)^2*td)*bess(nsum,j); sum=sum+term; end if (nsum < m) p(i)=pic+FK*((2/(red^2-1))*(rd^2/4+td)... -(red^2*log(rd))/(red^2-1)-redf... +pi*sum); else p(i)=pic end end title('Pressure History, x=x(1,1)') xlabel('Time, days') ylabel('Pressure, psi') hold on load prod1.dat tfd=prod1(:,1); pfd=prod1(:,4); plot(tfd,pfd,'kx') plot(t,p,'r') figure(4) % Plot with log time clf reset semilogx(tfd,pfd,'kx') title('Pressure History, x=x(1,1)') xlabel('Time, days') ylabel('Pressure, psi') hold on semilogx(t,pinf,'b-') semilogx(t,psq,'g-') semilogx(t,p,'r-') end