常规钻井钻进分析

%% 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=Tb