main_helical_dynamic.rar_时变啮合刚度_石川_石川法_齿轮_齿轮 刚度

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc;clear;close all; format long; %------------------------------------------------------------1、基本参数输入 Z_p=31;%齿数 Z_g=40; beita=31*pi/180;%螺旋角 alpha_n=17.5*pi/180;%法面分度圆压力角 m_n=1.75e-3;%法面模数 B=18e-3;%齿宽 ha_x_n=1.48;%法面齿顶高系数 c_x_n=0.3;%法面顶隙系数 x_n=0.2422; %法向变位系数 T_p=150;%输入扭矩Nm（T_in=9550*P(kW)/Nz1(r/min)） T_g=T_p*Z_g/Z_p; zeta=0.05;%阻尼比 b=40*1e-6;% %% %-------------------------------------------------------------2、参数计算 m=m_n/cos(beita);%端面模数 x_t=x_n*cos(beita); alpha=atan(tan(alpha_n)/cos(beita));%端面压力角 beita_b=atan(tan(beita)*cos(alpha));%基圆柱螺旋角 ha_x=ha_x_n*cos(beita);%端面齿顶高系数 c_x=c_x_n*cos(beita);%端面顶隙系数 h_a=ha_x_n*m_n;%齿顶高 h_f=(ha_x_n+c_x_n)*m_n;%齿根高 Pn=pi*m_n;%法向齿距 Pt=pi*m;%端面齿距 Pb_n=Pn*cos(alpha_n); Pb_t=Pt*cos(alpha);%基圆齿距 r_p=m*Z_p/2;%分度圆半径 r_g=m*Z_g/2; r_ap=r_p+h_a;%齿顶圆半径 r_fp=r_p-h_f;%齿根圆半径 r_ag=r_g+h_a; r_fg=r_g-h_f; r_bp=r_p*cos(alpha);%基圆半径 r_bg=r_g*cos(alpha); % r_op = 100e-3;%轴孔半径 % r_og = 90e-3; % rho = 7850; % I_p = rho*B*pi*((2*r_p)^4-(2*r_op)^4)/32; % I_g = rho*B*pi*((2*r_g)^4-(2*r_og)^4)/32; % m_p=pi*(r_p^2-r_op^2)*B*rho; % m_g=pi*(r_g^2-r_og^2)*B*rho; I_p = 265.2e-6;%小齿轮转动惯量 I_g = 540.4e-6;%大齿轮转动惯量 m_p = 0.513;%小齿轮质量 m_g = 0.582;%大齿轮质量 sp0 = xlsread('stiffness_wj.xlsx');%读取刚度值 k_m=mean(sp0(:,2));%求刚度的平均值（sp0中第二列所有数） wn=sqrt((r_bp^2/I_p+r_bg^2/I_g)*k_m);%固有频率，为根号下km/me,其中km为啮合刚度平均值，me为等效质量，即基圆半径二次方与转动惯量的比值。 C0=2*zeta*sqrt(k_m/(r_bp^2/I_p+r_bg^2/I_g));%阻尼 parameter1=[Z_p;beita;Z_p;T_p;T_g;b;zeta]; parameter2=[m_p;m_g;I_p;I_g;wn;C0;r_bp;r_bg]; N_p=linspace(400,18000,200);%主动轮的转速r/min,%linspace是Matlab中的一个指令，用于产生指定范围内的指定数量点数，相邻数据跨度相同，并返回一个行向量. %调用方法：linspace(x1,x2,N)功能：用于产生x1，x2之间的N点行矢量，相邻数据跨度相同。其中x1、x2、N分别为起始值、终止值、元素个数。若缺省N，默认点数为100。 NN=100;%求解的周期数 N_d=100;%每个周期内取的点个数 x0 = zeros(10,1);%创建10行1列的零矩阵 opt=odeset('abstol',1e-8,'reltol',1e-9);%解微分方程时的选项设置，abstol是绝对误差，reltol是相对误差 [t,x]=ode45('helicaldynamics_eqs',[0,1000],x0,opt,parameter1,parameter2,sp0,N_p(1)); x0 = x(end,:); x1rms = zeros(length(N_p),1); x2rms = zeros(length(N_p),1); %下面求解的是一个啮合周期内，如需多个周期，将tspan延长即可 XX1=zeros(50,200);XX2=zeros(50,200);XX3=zeros(50,200); tic for h = 1:length(N_p) TT1=60*wn/(Z_p*N_p(h));%h为对每个元素求绝对值，TT1为基频与啮合频率的比值，即无量纲啮合频率w=TT1=wk/wn tf=[0:TT1/500:NN*TT1]; [t,x]=ode45('helicaldynamics_eqs',tf,x0,opt,parameter1,parameter2,sp0,N_p(h)); x0 = x(end,:); DTE_y = x(:,1)-x(:,3)+x(:,5);D_DTE_y = x(6)-x(8)+x(10);%y方向动态传递误差 DTE_z = x(:,2)-x(:,4);D_DTE_z = x(7)-x(9);%z方向动态传递误差 DTE = DTE_y*cos(beita)+DTE_z*sin(beita);%齿轮副动态传递误差 D_DTE = D_DTE_y*cos(beita)+D_DTE_z*sin(beita); XX1(:,h)=x(50*500+230:500:end,5); XX2(:,h)=x(50*500+230:500:end,1); XX3(:,h)=DTE(50*500+230:500:end,1); % figure(121);plot(N_p(h),x(50*500+200:500:end,5),'k.');hold on; x1rms(h) = sqrt(sum(DTE(50*500:1:end,1).^2)/length(DTE(50*500:1:end,1)));%动态误差的均方根值 x2rms(h) = sqrt(sum(x(50*500:1:end,5).^2)/length(x(50*500:1:end,5))); h end toc save rms.mat N_p x1rms x2rms XX1 XX2 XX3 t x DTE D_DTE figure(11) %subplot(221) plot(N_p,x1rms,'-ko') %subplot(222) %plot(N_p,x2rms,'-ko') %subplot(223) %plot(N_p,XX3,'.k') %subplot(224) %plot(x(5000:500:end,5),x(5000:500:end,10),'*k') figure hold on plot(N_p/60*31,XX2) % kt = helical_stiffness(t,seta,ny,N_p(end),Z_p,wn);ct = C0; % et0 = 10*1e-6; et1 = 30*1e-6;et = (et0+et1*sin(w1*t)); % deta_y = x(:,1)+x(:,3)-x(:,4)+x(:,6);d_deta_y = x(:,7)+x(:,9)-x(:,10)+x(:,12); % deta_z = x(:,2)-x(:,5);d_deta_z = x(:,8)-x(:,11); % deta = deta_y.*cos(beita)+deta_z.*sin(beita); % d_deta = d_deta_y.*cos(beita)+d_deta_z.*sin(beita); % for i=1:length(deta) % if deta(i)-et(i)/b>1 % Fn(i) = kt(i)*(deta(i)-et(i)/b-1)+ct*wn*(d_deta(i)-et1*w1*cos(w1*t(i))/b); % elseif deta(i)-et(i)/b<-1 % Fn(i) = kt(i)*(deta(i)-et(i)/b+1)+ct*wn*(d_deta(i)-et1*w1*cos(w1*t(i))/b); % else % Fn(i) = 0; % end % end % Fm = Fn*b; % plot(t,Fm) % figure(2) % plot(t,b*x(:,1:6))