求助？关于直齿圆柱齿强度计算的问题

[yylyshyy]

程序如下：
function[sigmab sigmaf]=GearStress
global m BH fs YZ b Ft
JK=GearJofK;
Kv=GearKofV;
KH=GearKofH;
sigmab=Kv*KH*Ft/m/b/JK;
sigmafp=0.533*BH+88.3;
YN=GearYofN;
sigmaf=sigmafp*YN/fs/YZ;

function Kv=GearKofV
global Rp n Qv
vt=2*pi*Rp*n/60000;
B=0.25*(12-Qv).^(2/3);
A=50+56*(1-B);
vtmax=(A+Qv-3).^2/200;
if vt>vtmax
error('Maximum tangential velocity exceeded for given Qv')
end
Kv=((A+sqrt(200*vt))/A).^B;

function Kh=GearKofH
global b typeg Rp
class=[0.247 0.127 0.0675 0.0380;...
0.657e-3 0.622e-3 0.504e-3 0.402e-3;...
-1.186e-7 -1.69e-7 -1.44e-7 -1.27e-7];
Khma=class(1,typeg)+class(2,typeg)*b+class(3,typeg)*b.^2;
size(Khma);
ko=0.05*b/Rp;
if ko<0.05
ko=0.05;
end
if b<=25
Khpf=ko-0.025;
elseif b<=432
Khpf=ko-0.0375+0.000492*b;
else
Khpf=ko-0.1109+0.000815*b-0.353e-6*b.^2;
end

Kh=1+Khpf+Khma;

function GearParameters
global m phis ar rrT C N1 N2 rT1 rT2 ts Pb1 Pb2
global urmin urmax e xD rf yrp gammaw phi
xrp=-ar+rrT;
h=ar-rrT+rrT*sin(phis);
yrp=-pi*m/4+h*tan(phis)+rrT*cos(phis);
e=(ts-pi*m/2)/tan(phis)/2;
urmin=(e+xrp)/tan(phis)-yrp;
urmax=-yrp;
Pb=m*pi*cos(phis);
Rb1=N1*m/2*cos(phis);
Rb2=N2*m/2*cos(phis);
phi=acos((Pb1+Pb2)/C);
mc=(sqrt(rT1.^2-Rb1.^2)+sqrt(rT2.^2-Rb2.^2)-(Rb1+Rb2)*tan(phi))/Pb;
Rw=sqrt(Rb1.^2+(sqrt(rT1.^2-Rb1.^2)-(mc-1)*Pb).^2);
phiw=acos(Rb1/Rw);
thetaw=ts/m*N1+involute(phis)-involute(phiw);
gammaw=phiw-thetaw;
xD=Rw*cos(thetaw)-Rw*sin(thetaw)*tan(gammaw);
rf=rrT+(ar-e-rrT).^2/(N1*m/2+ar-e-rrT);


function JK=GearJofK

global urmin urmax phi
GearParameters
options=optimset('display','off');
ur=fminbnd(@GearKofJ,urmin,urmax,options);
JK=-cos(phi)/GearKofJ(ur);

function KJ=GearKofJ(ur)
global m rrT N1
global e xD xrp yrp gammaw
xip=e+xrp;
etap=ur+yrp;
s=1+rrT/sqrt(xip.^2+etap.^2);
xi=s*xip;
eta=s*etap;
thetaR=atan(eta/(N1*m/2+xi)-(ur-pi*m/2)/(N1*m/2));
R=sqrt((N1*m/2+xi).^2+eta.^2);
x=R*cos(thetaR);
y=R*sin(thetaR);
Kf=GearKofF(x,y);
KJ=-m*cos(gammaw)*Kf*(1.5*(xD-x)/y.^2-0.5*tan(gammaw)/y);


function KF=GearKofF(x,y)
global phis rf xD
d=phis*180/pi;
k1=0.3054-0.00489*d-0.000069*d.^2;
k2=0.362-0.01268*d+0.000104*d.^2;
k3=0.2934+0.00609*d+0.0000087*d.^2;
KF=k1+(2*y/rf).^k2*(2*y/(xD-x)).^k3;


function inv=involute(angle)
inv=tan(angle)-angle;


function YN=GearKofN(YN)
global n L BH
nL=60*n*L;
fBH=-9.2595e-6-6*(BH-160).^2+0.009722*(BH-160)+1.6;
if nL<=1e3
YN=fBH;
elseif nL<=3e6
D=0.8628*(log10(fBH)-0.0169);
E=-0.2876*(log10(fBH)-0.0169);
YN=(fBH*10.^D)*nL.^E;
else
YN=1.683*nL.^(-0.0323);
end




global m phis ar rrT ts C N1 N2 rT1 rT2 torque
global fs YZ b typeg Rp Ft
m=10;phis=20*pi/180;ar=12.5;ts=18.51; C=525; N1=28; N2=75;
rT1=153.9; rT2=391.1; torque=2400; BH=260;fs=1.2; YZ=1.0; b=45;
n=1800; typeg=2; L=4000;rrT=3.8;
Rp=C/(1+N2/N1);
Ft=1000*torque/Rp;
[sigmab,sigmaf]=GearStress;
disp(['The bending stress is' num2str(sigmab,5)' N/mm^2'])
disp(['The permissiable level is' num2str(sigmaf,5)' N/mm^2']);
执行后出现错误，
Error in ==> GearStress>GearKofH at 41
Kh=1+Khpf+Khma;

Error in ==> GearStress at 5
KH=GearKofH;


THANKS