shijing.m

用matlab语言编写的试井分析的程序

M=8;
CDS=[0.01,0.1,1,10,100,1000];
%tD=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700,800,900,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,20000,30000,40000,50000,60000,70000,80000,90000,100000,200000,300000,400000,500000,600000,700000,800000,900000,1000000];
for i=1:110
    tD(i)=0.001*10^(9*i/99);
end
L=length(tD);
%leng=length(tD);
%leng=length(tD);
for j=1:M
    %计算系数v(j)
    v(j)=0;
    for k=fix((j+1)/2):1:min(j,M/2)
        v(j) = v(j) + k ^ (M / 2) * factorial(2 * k) / (factorial(M / 2 - k) * factorial(k) * factorial(k - 1) * factorial(j - k) * factorial(2 * k - j));
    end
    v(j)=(-1)^(M/2+j)*v(j);
end
%反解Laplace空间压力到实空间
for k=1:6
   for i=1:L
      ft(i,k)=0;
      u=log(2)/tD(i);
      for j=1:M
          s=j*u;
          temp=sqrt(s/CDS(k));
          ft(i,k) = ft(i,k) + (v(j) / s) * besselk(0,temp) / (s * besselk(0,temp) + temp * besselk(1,temp));
      end
      ft(i,k)=ft(i,k)*u;
   end
end


%三点法求导数（试井分析意义下的压力导数）
for k=1:6
    D_ft(1,k) = tD(1) * (-3 * ft(1,k) + 4 * ft(2,k) - ft(3,k)) / (4 * tD(2) - 3 * tD(1) - tD(3));
    D_ft(2) = tD(2) * (ft(3,k) - ft(1,k)) / (tD(3) - tD(1));
    for i=3:L
        D_ft(i,k) = tD(i) * (ft(i - 2,k) - 4 * ft(i - 1,k) + 3 * ft(i,k)) / (tD(i - 2) - 4 * tD(i - 1) + 3 * tD(i));
    end
end


%两点法求导数（试井分析意义下的压力导数）
%for k=1:5
    %D_ft(1,k) = tD(1) * (ft(2,k)-ft(1,k))/(tD(2)-tD(1));
    %for i=2:64
        %D_ft(i,k)=tD(i)*(ft(i+1,k)-ft(i-1,k))/(tD(i+1)-tD(i-1));
    %end
    %D_ft(65,k)=tD(65) * (ft(65,k) - ft(64,k) ) / (tD(65) -  tD(64) );
%end
for k=1:6
    MD_ft(k)=max(D_ft(:,k))
 %   for i=1:L
  %      ft(i,k)=log(ft(i,k));
   %     D_ft(i,k)=log(D_ft(i,k));
    %end
end
loglog(tD,D_ft(:,1),tD,ft(:,1),tD,D_ft(:,2),tD,ft(:,2),tD,D_ft(:,3),tD,ft(:,3),tD,D_ft(:,4),tD,ft(:,4),tD,D_ft(:,5),tD,ft(:,5),tD,D_ft(:,6),tD,ft(:,6))
axis([10^(-2), 10^7, 10^(-2), 10^3]) 
xlabel('logtD/cD ','FontSize',16) 
ylabel('logPwD'' logPwD','FontSize',16) 
title('无限大均质油藏试井曲线','FontName','楷体' ,'FontSize',25)