z1 = 35;                %齿数z2 = 25;

m = 4;                 %模数

a = 20;                %齿轮压力角

c = 0.25;              %顶隙系数

ha = 1;                  %齿顶高系数

beta = 20;               %啮合线与x轴的夹角???

E = 2*10^9;              %弹性模量

v = 0.3;                 %泊松比

L = 40;                   %齿宽

r_b1 = m*z1*cosd(a);         %基圆半径

r_b2 = m*z2*cosd(a);

%

inva = tand(a)-deg2rad(a);      %渐开线函数

theta_b1 = pi/(2*z1)+inva;         %代表基齿角的一半

theta_b2 = pi/(2*z2)+inva;

x_beta1 = r_b1*((deg2rad(beta)+theta_b1)*cosd(beta)-sind(beta));

y_beta1 = r_b1*((deg2rad(beta)+theta_b1)*sind(beta)+cosd(beta));

x_beta2 = r_b1*((deg2rad(beta)+theta_b2)*cosd(beta)-sind(beta));

y_beta2 = r_b1*((deg2rad(beta)+theta_b2)*sind(beta)+cosd(beta));

%

r_C1 = sqrt((r_b1*tand(a)-ha*m/sind(a))^2+r_b1^2);

alpha_C1 = acos(r_b1/r_C1 );       %弧度

invalpha_C1 = tan(alpha_C1)-alpha_C1;

theta_C1 = theta_b1-invalpha_C1;

tau_C1 = alpha_C1-theta_C1;

y_c1 = r_b1*((tau_C1+theta_b1)*sin(tau_C1)+cos(tau_C1)); %渐开线起始点到原点之间的

%水平距离

r_C2 = sqrt((r_b2*tand(a)-ha*m/sind(a))^2+r_b1^2);

alpha_C2 = acos(r_b2/r_C2 );       %弧度

invalpha_C2 = tan(alpha_C2)-alpha_C2;

theta_C2 = theta_b2-invalpha_C2;

tau_C2 = alpha_C2-theta_C2;

y_c2 = r_b2*((tau_C2+theta_b2)*sin(tau_C2)+cos(tau_C2)); %渐开线起始点到原点之间的

%水平距离

% 定义变量tau和gamma

syms gamma

r1 = m*z1/2;

r2 = m*z2/2;                  %分度圆半径

r_p = c*m/(1-sind(a));        % 为刀顶圆角半径

a_1 =(ha+c)*m-r_p;            %刀顶圆角圆心距中线的距离

b_1 =pi*m/4+ha*m*tan(a)+r_p*cosd(a);   %为刀顶圆角圆心距刀具齿槽中心线的距离

phi1 =(a_1/tan(gamma)+b_1)/r1;  %

phi2 =(a_1/tan(gamma)+b_1)/r2;  %

%齿轮1

x_1 = r1*sin(phi1)-(a_1/sin(gamma)+r_p)*cos(gamma-phi1);

y_1 = r1*cos(phi1)-(a_1/sin(gamma)+r_p)*sin(gamma-phi1);

x_2 = r_b1*((gamma+theta_b1)*cos(gamma)-sin(gamma));

y_2 = r_b1*((gamma+theta_b1)*sin(gamma)-cos(gamma));

A_y1 = 2*x_1*L; A_y2 = 2*x_2*L; I_y1 = 2/3*x_1^3*L; I_y2 = 2/3*x_2^3*L;

%齿轮2

x_3 = r2*sin(phi2)-(a_1/sin(gamma)+r_p)*cos(gamma-phi2);

y_3 = r2*cos(phi2)-(a_1/sin(gamma)+r_p)*sin(gamma-phi2);

x_4 = r_b2*((gamma+theta_b2)*cos(gamma)-sin(gamma));

y_4 = r_b2*((gamma+theta_b2)*sin(gamma)-cos(gamma));

A_y3 = 2*x_3*L; A_y4 = 2*x_4*L; I_y3 = 2/3*x_3^3*L; I_y4 = 2/3*x_4^3*L;

%

dy_1_dgamma = (a_1*sin((a_1/tan(gamma)+b_1)/r1)*(1+(tan(gamma))^2))/(tan(gamma))^2 ...

+(a_1*cos(gamma))/((sin(gamma))^2)*sin(gamma-(a_1/tan(gamma)+b_1)/r1) ...

-(a_1/(sin(gamma))^2 +r_p)*cos(gamma-(a_1/tan(gamma)+b_1)/r1) ...

*(1+(a_1*(1+(tan(gamma))^2)/(r1*(tan(gamma))^2)));

dy_2_dtau = r_b1*(gamma+theta_b1)*cos(gamma);

dy_3_dgamma = (a_1*sin((a_1/tan(gamma)+b_1)/r2)*(1+(tan(gamma))^2))/(tan(gamma))^2 ...

+(a_1*cos(gamma))/((sin(gamma))^2)*sin(gamma-(a_1/tan(gamma)+b_1)/r2) ...

-(a_1/(sin(gamma))^2 +r_p)*cos(gamma-(a_1/tan(gamma)+b_1)/r2) ...

*(1+(a_1*(1+(tan(gamma))^2)/(r2*(tan(gamma))^2)));

dy_4_dtau = r_b2*(gamma+theta_b2)*cos(gamma);

%弯曲刚度

a1 = deg2rad(a);

beta1 = deg2rad(beta);

f_11 = (cosd(beta)*(y_beta1-y_1)-x_beta1*sind(beta))^2/(E*I_y1)*dy_1_dgamma;

f_12 = (cosd(beta)*(y_beta1-y_2)-x_beta1*sind(beta))^2/(E*I_y2)*dy_2_dtau;

k_b1 = 1/(int(f_11,pi/2,a1)+int(f_12,tau_C1,beta1));

f_13 = (cosd(beta)*(y_beta2-y_3)-x_beta2*sind(beta))^2/(E*I_y3)*dy_3_dgamma;

f_14 = (cosd(beta)*(y_beta2-y_4)-x_beta2*sind(beta))^2/(E*I_y4)*dy_4_dtau;

k_b2 = 1/(int(f_13,pi/2,a1)+int(f_14,tau_C2,beta1));

%剪切刚度

G = E/(2*(1+v));               %为剪切模量

f_21 = (cosd(beta))^2/(G*A_y1)*(dy_1_dgamma);

f_22 = (sind(beta))^2/(G*A_y2)*(dy_2_dtau);

k_s1 = 1/(int(f_21,pi/2,a1)+int(f_22,tau_C1,beta1));

f_23 = (cosd(beta))^2/(G*A_y3)*(dy_3_dgamma);

f_24 = (sind(beta))^2/(G*A_y4)*(dy_4_dtau);

k_s2 = 1/(int(f_23,pi/2,a1)+int(f_24,tau_C2,beta1));

%压缩刚度

f_31 = (sind(beta))^2/(E*A_y1)*dy_1_dgamma;

f_32 = (sind(beta))^2/(E*A_y2)*dy_2_dtau;

k_a1 = 1/(int(f_31,pi/2,a1)+int(f_32,tau_C1,beta1));

f_33 = (sind(beta))^2/(E*A_y3)*dy_3_dgamma;

f_34 = (sind(beta))^2/(E*A_y4)*dy_4_dtau;

k_a2 = 1/(int(f_33,pi/2,a1)+int(f_34,tau_C2,beta1));

%接触刚度与齿轮轮体刚度

k_h = 6.9011e+07;

k_f1 = 6.5893e+09;

k_f2 = 8.1926e+09;

%一对啮合轮齿的总刚度

K = 1/(1/k_h+1/k_b1+1/k_s1+1/k_a1+1/k_f1+1/k_b2+1/k_s2+1/k_a2+1/k_f2);

plot(gamma,K)