常规钻井钻进分析

 %% time frame and time increament

    
 tau=2*pi/n/omiga0;
    
 t=-80*tau:h:100;
    
 %% matrix construction
    
 y1=zeros(size(t));y2=zeros(size(t));y3=zeros(size(t));y4=zeros(size(t));
    
 tao=zeros(size(t));
    
 %% initial value
    
 m=tau/h;
    
 for i=1:80*m+1
    
     y1(i)=v0;y2(i)=v0*t(i);y3(i)=omiga0;y4(i)=omiga0*t(i);
    
     tao(i)=m;
    
 end
    
 %% main program
    
 for j=80*m+1:length(t)-1
    
     %% equation 1
    
     A=[];
    
     kk=1:80*m;
    
     N=[A abs(y4(j)-y4(j-kk)-2*pi/n)];
    
     index=find(N==min(N));
    
     tao(j)=kk(:,index);
    
     d=n*(y2(j)-y2(j-tao(j)));
    
     Wbc=KFc*d*H(d)*H(y3(j));
    
     Wbf=KFf0*H(d)*H(y1(j));
    
     Wb=Wbc+Wbf;
    
     Fd=F0*cos(200*pi*t(j));
    
     y1(j+1)=y1(j)+h*(W0+Fd-Wb-Ca*(y1(j)-v0)-Ka*(y2(j)-v0*t(j)))/m1;
    
     %% equation 2
    
     y2(j+1)=y2(j)+h*y1(j);
    
     %% equation 3
    
     Tbc=KTc*d*H(d)*H(y3(j));
    
     Tbf=KTf0*sign(y3(j))*H(d)*H(y1(j));
    
     Tb=Tbc+Tbf;
    
     y3(j+1)=y3(j)+h*(-Tb-Ct*(y3(j)-omiga0)-Kt*(y4(j)-omiga0*t(j)))/I1;
    
     %% equation 4
    
     y4(j+1)=y4(j)+h*y3(j);
    
 end