function [tj1,tj2,Ja,M2a,M2b,Wa,Wb,Na,Nb]=clutchengagement(w1,w2)
beta=2.2; %输入离合器储备系数
elta=0.1; %输入循环液粘度系数
rho=1220; %输入循环液密度
m1=6.04; %输入每米钻杆质量
f=0.4; %输入钻头与岩石摩擦系数
G=10; %输入钻压
lamda=2.2; %输入电机过载倍数
PN=55; %输入电机额定功率
nN=1480; %输入电机额定转速
n2=125; %输入钻杆转速
phi3=56.5; %输入钻孔直径
phi2=50; %输入钻杆外径
phi1=44; %输入钻杆内径
L=2000; % 输入钻孔深度
J1=0.66; %55kw电机转子转动惯量
J2=0.625; %55kw电机联轴器转动惯量
Je=J1+J2; %主动部分转动惯量
icr=nN/n2; %离合器从动盘到钻柱传动比icr
J3=pi*rho*L*phi1^4/(32*10^12);%钻井液转动惯量
J4=m1*L*(phi2^2+phi1^2)/(8*10^6);%钻柱转动惯量
Ja=(J3+J4)/icr^2;%等效到离合器从动盘转动惯量
TN=9550*PN/nN; %计算电动机额定输出转矩
h=0.01; %时间步长0.01秒
w1=pi*nN/30; %离合器主动盘初始角速度值pi*nN/30
w2=0; %离合器从动盘角速度初始值为0
a1=0; a2=0; %主、从动盘角加速度初始值为0
Lh=0; %滑磨功初始值为0
t=0; %离合器接合时间初始值为0
while w2
Mm=Mm(t+h);
a1=a1_1;
w1=w1+a1*h; %递推得出主动盘下一个角速度值
a2=a2_1;
w2=w2+a2*h; %递推得出从动盘下一个角速度值
else Mm=beta*TN;
a1=a1_2;
w1=w1+a1*h;
a2=a2_2;
w2=w2+a2*h;
end
wh=w1-w2; %主、从动盘角速度之差
Nh=Mm*wh; %瞬时滑磨功率
Lh=Lh+Nh*h; %滑磨功
disp(w1); disp(w2);disp(Nh); disp(Lh); %输出中间时刻各参数
t=t+h;
end %w1=w2退出循环,并输出此时的时间t
disp(t); %输出接合时间
function M1=M1(w1) %主动部分输入转矩子函数
M1=-13.7422*w1^2+4119.6*w1-307982;
function Mm=Mm(t) %离合器摩擦力矩子函数
if t<0.2
Mm=t*beta*TN/0.2
else t>0.2
Mm=beta*TN
function a1=a1_1(t,w1) %t<0.2秒时主动盘角加速度子函数
a1=(-13.7422*w1^2+4119.6*w1-307982-t*beta*TN/0.2)/Je
function a1=a1_2(t,w1) %t>0.2秒时主动盘角加速度子函数
a1=(-13.7422*w1^2+4119.6*w1-307982-beta*TN)/Je
function a2=a2_1(t,w2) %t<0.2秒时从动盘角加速度子函数
a2=(t*beta*TN/0.2-elta*pi*L*w2*phi2^3/(2*(phi3-phi2)*icr*10^6)+G*f*(phi1^2+phi1*phi2+phi2^2)/(3*icr*(phi1+phi2)))/Ja
function a2=a2_2(t,w2) %t>0.2秒时从动盘角加速度子函数
a2=(beta*TN-elta*pi*L*w2*phi2^3/(2*(phi3-phi2)*icr*10^6)+G*f*(phi1^2+phi1*phi2+phi2^2)/(3*icr*(phi1+phi2)))/Ja