关于ODE45积分器的一些代码实现:
[t,theta] = ode45(@(t,theta) odefun(t,theta,omega_pz,GearSet,slices_num,k_single,J_wz,r_wb,Tin,Transmission_ratio), t_span, [theta_wz_0,omga_wz_0]);%%注意这里的t是系统自动变化的
% out = [t,theta];
% legend('wheel theta in rad','omega in rad per s','pinion theta in rad')
legend('wheel theta in rad','pinion theta in rad')
%
% plot(1:length(t)/381:length(t),t*omega_pz,'-o')
figure(3)
TE=(Pinion_ang-Wheel_ang*Transmission_ratio)*db(1)/2*1e3;%[um]
plot(t,TE)
xlabel('Time in s')
ylabel('TE in um')
%%
function theta_dt = odefun(t,theta,omega_pz,GearSet,slices_num,k_single,J_wz,r_wb,Tin,Transmission_ratio)%%注意这里的t是系统自动变化的
r_pz_all = [];
r_wz_all = [];
delt_L = [];
dir = Tin/abs(Tin)
theta_input = dir*(omega_pz*t)*180/pi%这里输入是degree
Points_data = Cal_Contact_Points_th_Fun(GearSet,slices_num,theta_input);
r_pz_all = Points_data(:,4);%小齿轮 接触点半径
r_wz_all = Points_data(:,5);%大齿轮接触点半径
for i = length(r_pz_all)
delt_L(i,1) = omega_pz*t*r_pz_all(i,1)-theta(1)*r_wz_all(i,1);%%注意这里的t是系统自动变化的
if delt_L(i,1)>0
delt_L(i,1)=delt_L(i,1);
else
delt_L(i,1)=0;
end
end
theta_dt = zeros(2,1);
theta_dt(1) = theta(2);%θ的一阶导数是ω这里的theta(2)是角速度
theta_dt(2) = (-dir*k_single*sum(delt_L)*r_wb/J_wz+Tin*Transmission_ratio/J_wz);%θ的二阶导数是α
figure(1)
hold on
plot(t,(-dir*k_single*sum(delt_L)*r_wb/J_wz+Tin*Transmission_ratio/J_wz),'or','MarkerSize',10,'MarkerFaceColor','r')
xlabel('time in s')
ylabel('T_e_x in Nmm')
end